马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
本帖最后由 benyou 于 2012-6-15 17:30 编辑
由于matlab不熟,现在写论文编了一个程序不知道为什么,我改变参数10到48(也就是矩阵的维数)就会提示Index exceeds matrix dimensions。谁能帮忙指点下错误的地方,程序如下:
function main
clear all
rand('state',0);
randn('state',0);
N_dot=48;%改成10就行了可是我就得用到48乘以48的啊,怎么办?
p=1000;
ax=[106.382 119.036 131.365 143.158 154.211 164.338 173.363 181.133 187.515 192.399 195.703 197.369 197.369 195.703 192.399 187.515 181.133 173.363 164.338 154.211 143.158 131.365 119.036 106.382 93.618 80.964 68.635 56.842 45.789 35.662 26.637 18.867 12.485 7.601 4.297 2.631 2.631 4.297 7.601 12.485 18.867 26.637 35.662 45.789 56.842 68.635 80.964 93.618];
az=[10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 130 135 140 145 150 155 160 165 170 175 180 185 190 195 200 205 210 215 220 225 230 235 240 245];
v10=30;
alpha=0.22;
k=0.00464;
delta_t=0.2;
Uz=v10*(az/10).^alpha; %产生功率谱
vvd=0.4*v10/2.302585;
N=p;
posi=zeros(N_dot,N_dot,N);
ch=zeros(N_dot,N_dot,N);
S=zeros(N_dot,N_dot,N);
H=zeros(N_dot,N_dot,N);
dw=10*pi;
dltw=dw/N;
for m=1:1:N
mid1_S_H=zeros(N_dot,N_dot);
mid2_S_H=zeros(N_dot,N_dot);
n(m)=dw*m/N;%频率
f(m)=1200*n(m)/v10/2/pi;
Pn(m)=4*k*v10^2*f(m).^2/(1+f(m).^2).^(4/3);% 风速谱,如计算竖向Panofsky 谱将其替换为相应函数Uf=1;%摩擦速度
Su(m)=Pn(m)/n(m);
for g=1:N_dot
for h=1:N_dot
if g==h
ch(g,h,m)=1;
else
delta_z=sqrt(256*(ax(g)-ax(h)).^2+100*(az(g)-az(h)).^2);
ch(g,h,m)=exp(-2*n(m)/2/pi*delta_z/(Uz(g)+Uz(h)));
f_star=2*n(m)/2/pi*abs(az(g)-az(h))/(Uz(g)+Uz(h));
if f_star<=0.1
posi(g,h,m)=pi*f_star/4;
else if f_star<=0.125
posi(g,h,m)=-10*pi*f_star+1.25;
else posi(g,h,m)=pi-2*pi*rand(1);
end
end
end
end
end
for g=1:N_dot
for h=1:N_dot
S(g,h,m)=Su(m)*ch(g,h,m)*exp(i*posi(g,h,m));
mid1_S_H(g,h)=S(g,h,m);
end
end
[mid2_S_H,a]=chol(mid1_S_H);
mid2_S_H=mid2_S_H';
clear g;
clear h;
for g=1:N_dot
for h=1:N_dot
H(g,h,m)=mid2_S_H(g,h);
end
end
end |