我在MATLAB中编写了这段代码来存储矩阵。我使用的是cells数组,但在Python中我不知道如何做到这一点。
有人知道我该怎么做吗?
MATLAB代码为:
S =[0.5 0.7 0.9 1.1]; % distância entre tx e rx[m]
d = 0.2*ones(1,10);
h = [ 0 0.1 0.2 0.3 0.4]
n_cam = length(d); % numero de camadas
n_alt = length(h); % numero de alturas
n_S = length(S); % numero de receptores
z = zeros(n_cam,n_alt); % profundidade
Rv_h = zeros(n_S,n_alt);
Rv_v = zeros(n_S,n_alt);
Rv = zeros(n_cam,n_alt);
Rh = zeros(n_cam,n_alt);
S_Rv = cell(1,n_S);
S_Rh = cell(1,n_S);
sigma = 0.3*ones(1,n_cam);
sigmaah = zeros(n_S,n_alt);
for i = 1:n_S
for j = 1:n_alt
for k = 1:n_cam
z(k,j)= (sum(d(1:k))+h(j))/S(i);
Rv(k,j) = 1/((4*z(k,j)^2+1)^0.5);
Rh(k,j) = ((4*z(k,j)^2+1)^0.5)-2*z(k,j);
end
Rv_h(i,j) = 1/((4*(h(j)/S(i))^2+1)^0.5);
Rh_h(i,j)=((4*(h(j)/S(i))^2+1)^0.5)-2*(h(j)/S(i));
end
S_Rv(:,i) = {Rv}; % z para cada camada em cada altura, para cada S
S_Rh(:,i) = {Rh};
end
for i = 1:n_S
for j = 1:n_alt
Rv = cell2mat(S_Rv(1,i));
Rh = cell2mat(S_Rh(1,i));
sigma_ah = sigma(1)*(Rh_h(i,j)-Rh(1,j));
sigma_av = sigma(1)*(Rv_h(i,j)-Rv(1,j));
for k = 2:(n_cam-1)
sigma_ah_ant = sigma_ah;
sigma_av_ant = sigma_av;
sigma_ah = sigma_ah_ant + sigma(k)*(Rh(k-1,j)-Rh(k,j));
sigma_av = sigma_av_ant + sigma(k)*(Rv(k-1,j)-Rv(k,j));
end
sigmaah (i,j) = sigma_ah + sigma(end)*Rh(n_cam-1,j)
sigmaav (i,j) = sigma_av + sigma(end)*Rv(n_cam-1,j)
end
end
我在想,在Python中,我可以做一些类似的事情:
n_S = 4
n_alt = 9
n_cam = 6
Rv =[]
for i in range(1,n_S):
for j in range(1,n_alt):
for k in range(1,n_cam):
z[k][j]= (sum(d[0:k])+h[j])/S[i]
Rv[i][j][k] = 1/((4*z[k,j]**2+1)**0.5)
但它不起作用,我收到的错误消息是
列表索引超出范围。
我建议使用广播功能
然后可以更换环路以澄清代码:
from numpy import array, ones
S =array([0.5, 0.7, 0.9, 1.1])
d = 0.2*ones((10));
h = array([ 0, 0.1, 0.2, 0.3, 0.4])
z = ((d.cumsum()[:, None] + h).ravel() / S[:, None]).reshape((S.size, d.size, h.size))
Rv = 1 / (4 * z ** 2 + 1)** .5
# ...etc
print(Rv[-1])
哪个输出:
[[0.93979342 0.87789557 0.80873608 0.73994007 0.67572463]
[0.80873608 0.73994007 0.67572463 0.61782155 0.56652882]
[0.67572463 0.61782155 0.56652882 0.52145001 0.48191875]
[0.56652882 0.52145001 0.48191875 0.4472136 0.41665471]
[0.48191875 0.4472136 0.41665471 0.38963999 0.36565237]
[0.41665471 0.38963999 0.36565237 0.34425465 0.32507977]
[0.36565237 0.34425465 0.32507977 0.30782029 0.29221854]
[0.32507977 0.30782029 0.29221854 0.27805808 0.26515648]
[0.29221854 0.27805808 0.26515648 0.25335939 0.24253563]
[0.26515648 0.25335939 0.24253563 0.23257321 0.22337616]]
与倍频程/matlab中的计算重叠:
Rv =
0.93979 0.87790 0.80874 0.73994 0.67572
0.80874 0.73994 0.67572 0.61782 0.56653
0.67572 0.61782 0.56653 0.52145 0.48192
0.56653 0.52145 0.48192 0.44721 0.41665
0.48192 0.44721 0.41665 0.38964 0.36565
0.41665 0.38964 0.36565 0.34425 0.32508
0.36565 0.34425 0.32508 0.30782 0.29222
0.32508 0.30782 0.29222 0.27806 0.26516
0.29222 0.27806 0.26516 0.25336 0.24254
0.26516 0.25336 0.24254 0.23257 0.22338
这减少了厄运的金字塔,并且由于愚蠢的魔法,可能比for循环更快。编辑:格式化第2版:给检查
在执行此操作之前,需要将z
和Rv
定义为已知大小的2D和3D阵列。
还要注意,python(和numpy(数组是基于零的,所以只需直接使用range
。
未测试:
import numpy as np
d = 0.2*np.arange(10) ### not sure if this is what you meant
h = np.array([ 0 0.1 0.2 0.3 0.4])
n_S = 4
n_alt = 9
n_cam = 6
Rv = np.zeros((n_S, n_alt, n_cam))
z = np.zeros((n_cam, n_alt))
for i in range(n_S):
for j in range(n_alt):
for k in range(n_cam):
z[k][j]= (sum(d[0:k])+h[j])/S[i] ## not sure about 0:10
Rv[i][j][k] = 1/((4*z[k,j]**2+1)**0.5)
然而,正如@GlobalTraveler所指出的,做到这一点的最愚蠢/最愚蠢的方法是利用广播,而不是使用循环:
Rv = 1/np.sqrt(4 * z**2 + 1)