clear

tic

v10=24;

z=116;

n=0.001:0.001:10;

xn=z*n./v10;

k=0.003;
%地面粗糙长度

ti=0.1;%时间步长

s1=200*k*v10^2*xn./n./(1+50*xn).^(5/3);%K谱

%产生空间点坐标

for i=1:20

x(i)=5+i;

z(i)=8+i;

end

%求R矩阵

syms f

R0=zeros(20);

for i=1:20

for j=i:20

H0=inline('200*k* v10^2*116*f./v10./f./(1+(50*116*f./v10)).^(5/3).*exp(-sqrt(dx^2/50^2+dz^2/60^2))','f','k','dx','dz','v10');

dx=x(i)-x(j);

dz=z(i)-z(j);

R0(i,j)=quadl(H0,0.001,10,0.001,0,k,dx,dz,v10);

R0(j,i)=R0(i,j);

end

end

R1=zeros(20);

for i=1:20

for j=i:20

H1=inline('200*k* v10^2*116*f./v10./f./(1+(50*116*f./v10)).^(5/3).*exp(-sqrt(dx^2/50^2+dz^2/60^2)).*cos(2*pi*f*1*ti)','f','k','dx','dz','ti','v10');

dx=x(i)-x(j);

dz=z(i)-z(j);

R1(i,j)=quadl(H1,0.001,10,0.001,0,k,dx,dz,ti,v10);

R1(j,i)=R1(i,j);

end

end

R2=zeros(20);

for i=1:20

for j=i:20

H2=inline('200*k* v10^2*116*f./v10./f./(1+(50*116*f./v10)).^(5/3).*exp(-sqrt(dx^2/50^2+dz^2/60^2)).*cos(2*pi*f*1*ti)','f','k','dx','dz','ti','v10');

dx=x(i)-x(j);

dz=z(i)-z(j);

R2(i,j)=quadl(H2,0.001,10,0.001,0,k,dx,dz,ti,v10);

R2(j,i)=R2(i,j);

end

end

R3=zeros(20);

for i=1:20

for j=i:20

H3=inline('200*k* v10^2*116*f./v10./f./(1+(50*116*f./v10)).^(5/3).*exp(-sqrt(dx^2/50^2+dz^2/60^2)).*cos(2*pi*f*1*ti)','f','k','dx','dz','ti','v10');

dx=x(i)-x(j);

dz=z(i)-z(j);

R3(i,j)=quadl(H3,0.001,10,0.001,0,k,dx,dz,ti,v10);

R3(j,i)=R3(i,j);

end

end

R4=zeros(20);

for i=1:20

for j=i:20

H4=inline('200*k* v10^2*116*f./v10./f./(1+(50*116*f./v10)).^(5/3).*exp(-sqrt(dx^2/50^2+dz^2/60^2)).*cos(2*pi*f*1*ti)','f','k','dx','dz','ti','v10');

dx=x(i)-x(j);

dz=z(i)-z(j);

R4(i,j)=quadl(H4,0.001,10,0.001,0,k,dx,dz,ti,v10);

R4(j,i)=R4(i,j);

end

end

Q1=zeros(20);

Q2=zeros(20);

Q3=zeros(20);

Q4=zeros(20);

A=[R0 R1 R2 R3;R1 R0 R1 R2;R2 R1 R0 R1;R3 R2 R1 R0];%80X80矩阵

B=[R1;R2;R3;R4];%80X20矩阵

X=A\B;
%此式相当于A*X=B，X(80×20矩阵)为自回归系数Ψ

q1=X(1:20,:);%取X的第一个20×20矩阵

q2=X(20+1:2*20,:);%取X的第二个20×20矩阵

q3=X(2*20+1:3*20,:);%取X的第三个20×20矩阵

q4=X(3*20+1:4*20,:);%取X的第四个20×20矩阵

Q1=q1';

Q2=q2';

Q3=q3';

Q4=q4';

RN=R0+Q1*R1+Q2*R2+Q3*R3+Q4*R4;%计算时出现

L=zeros(20);

L=chol(RN);

L=L';

a=zeros(20,2048);

for i=1:20

for j=1:2048

a(i,j)=normrnd(0,1,1,1);
%产生均值0，方差1的正态随机数矩阵

end

end

V(1:20,1)=L*a(:,1);

V(1:20,2)=Q1*V(1:20,1)+L*a(:,2);

V(1:20,3)=(Q1*V(1:20,2)+Q2*V(1:20,1))+L*a(:,3);

V(1:20,4)=(Q1*V(1:20,3)+Q2*V(1:20,2)+Q3*V(1:20,1))+L*a(:,4);

for t=5:2048

V(1:20,t)=(Q1*V(1:20,t-1)+Q2*V(1:20,t-2)+Q3*V(1:20,t-3)+Q4*V(1:20,t-4))+L*a(:,t);

end

V1=V(2,:);%取第一点的风速

t=(1:2048)*ti;

figure

subplot(2,1,1);

plot(t,V1,'b-');

xlabel('t(s)');

ylabel('v(t)');

axis([0 200 -60 60]);

%与目标谱进行比较

%[power,freq]=psd(V1,2048,10,boxcar(1024),0,'mean');

[power,freq]=cpsd(V1,V1,blackman(1024),0.5,2048,10);

power=power.*2./10;%规一化修正

subplot(2,1,2);

loglog(freq,power,'r-',n,s1,'g-');

xlabel('freq');

ylabel('s1/power');

toc