close all

clear all

clc

x=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 ];

x=x';

y=[4.02725 3.3011 2.8372 2.4234 2.06155 1.8402 1.7439 1.645 1.546 1.44 1.3349 1.2369 1.1397 1.0046];

y=y';

z=[3.75 9.75 13.75 17.25 21.5 25.25 27.825 30.475 33.125 35.95 38.775 41.4 44 46.05];

z=z';

h=14;%节点个数

tqd=1;

v10=29.665;%10米处风速

k=1:h;

v(k)=v10*(z(k)/10).^0.16;  %任意高度处的风速指数风速廓线

Tstep=6000;%步长的个数

dt=0.1;%步长的大小

f=0.01:0.01:10;%脉动风频率

u=0.005*v10^2;

P=4;

t=(1:6000)*dt;

nt=length(t);

%以下是求RO~R4

R0=zeros(h);

for i=1:h

for j=i:h

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

dy=y(i)-y(j);

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

h1=z(i)*f/v10/(z(i)/10)^0.16;

h2=z(j)*f/v10/(z(j)/10)^0.16;

w1=f.*(1+50*h1).^(5/3);

w2=f.*(1+50*h2).^(5/3);

s1=200*h1*u./w1;  %求kaimal谱

s2=200*h2*u./w2;

xiangguanhanshu=exp(-sqrt(dx.^2./50^2+dy.^2./50^2+dz.^2./60^2)).*cos(2*pi.*f*0.*dt);   %求相关函数

y=sqrt(s1.*s2).*xiangguanhanshu;  %求互功率谱水平顺风向和竖向互谱

R0(i,j)=trapz(f,y);

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

end

end

R1=zeros(h);

for i=1:h

for j=i:h

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

dy=y(i)-y(j);

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

h1=z(i)*f/v10/(z(i)/10)^0.16;

h2=z(j)*f/v10/(z(j)/10)^0.16;

w1=f.*(1+50*h1).^(5/3);

w2=f.*(1+50*h2).^(5/3);

s1=200*h1*u./w1;

s2=200*h2*u./w2;

xiangguanhanshu=exp(-sqrt(dx.^2./50^2+dy.^2./50^2+dz.^2./60^2)).*cos(2*pi.*f*1.*dt);

y=sqrt(s1.*s2).*xiangguanhanshu;

R1(i,j)=trapz(f,y);

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

end

end

R2=zeros(h);

for i=1:h

for j=i:h

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

dy=y(i)-y(j);

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

h1=z(i)*f/v10/(z(i)/10)^0.16;

h2=z(j)*f/v10/(z(j)/10)^0.16;

w1=f.*(1+50*h1).^(5/3);

w2=f.*(1+50*h2).^(5/3);

s1=200*h1*u./w1;

s2=200*h2*u./w2;

xiangguanhanshu=exp(-sqrt(dx.^2./50^2+dy.^2./50^2+dz.^2./60^2)).*cos(2*pi.*f*2.*dt);%Shiotani谱

y=sqrt(s1.*s2).*xiangguanhanshu;

R2(i,j)=trapz(f,y);

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

end

end

R3=zeros(h);

for i=1:h

for j=i:h

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

dy=y(i)-y(j);

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

h1=z(i)*f/v10/(z(i)/10)^0.16;

h2=z(j)*f/v10/(z(j)/10)^0.16;

w1=f.*(1+50*h1).^(5/3);

w2=f.*(1+50*h2).^(5/3);

s1=200*h1*u./w1;

s2=200*h2*u./w2;

xiangguanhanshu=exp(-sqrt(dx.^2./50^2+dy.^2./50^2+dz.^2./60^2)).*cos(2*pi.*f*3.*dt);

y=sqrt(s1.*s2).*xiangguanhanshu;

R3(i,j)=trapz(f,y);

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

end

end

R4=zeros(h);

for i=1:h

for j=i:h

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

dy=y(i)-y(j);

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

h1=z(i)*f/v10/(z(i)/10)^0.16;

h2=z(j)*f/v10/(z(j)/10)^0.16;

w1=f.*(1+50*h1).^(5/3);

w2=f.*(1+50*h2).^(5/3);

s1=200*h1*u./w1;

s2=200*h2*u./w2;

xiangguanhanshu=exp(-sqrt(dx.^2./50^2+dy.^2./50^2+dz.^2./60^2)).*cos(2*pi.*f*4.*dt);

y=sqrt(s1.*s2).*xiangguanhanshu;

R4(i,j)=trapz(f,y);

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

end

end

Q1=zeros(h); %14*14阶零阵

Q2=zeros(h);

Q3=zeros(h);

Q4=zeros(h);

A=[-R0 -R1 -R2 -R3;-R1 -R0 -R1 -R2;-R2 -R1 -R0 -R1;-R3 -R2 -R1 -R0];

B=[R1;R2;R3;R4];

X=inv(A)*B;

q1=X(1:h,:);

q2=X(h+1:2*h,:);

q3=X(2*h+1:3*h,:);

q4=X(3*h+1:4*h,:);

Q1=q1';

Q2=q2';

Q3=q3';

Q4=q4';

RN=R0+Q1*R1+Q2*R2+Q3*R3+Q4*R4;  %求RN，协方差矩阵

L=chol(RN); %对RN柯西分解，求L

L=L';

a=zeros(h,nt);

a=normrnd(0,1,h,nt);  %构造 以0为均值，以1为方差，h*nt阶彼此相互独立的 的正态分布的随机函数

v(1:h,1)=L*a(:,1);

v(1:h,2)=-Q1*v(1:h,1)+L*a(:,2);

v(1:h,3)=-(Q1*v(1:h,2)+Q2*v(1:h,1))+L*a(:,3);

v(1:h,4)=-(Q1*v(1:h,3)+Q2*v(1:h,2)+Q3*v(1:h,1))+L*a(:,4);

for it=5:nt  %14个空间点相关脉动风速时程列向量的AR模型

v(1:h,it)=-(Q1*v(1:h,it-1)+Q2*v(1:h,it-2)+Q3*v(1:h,it-3)+Q4*v(1:h,it-4))+L*a(:,it);

end

v;

%以下是平均风加脉动风,得总风速

fs1=v(1,:)+v10*(z(1)/10).^0.16;

fs2=v(2,:)+v10*(z(2)/10).^0.16;

fs3=v(3,:)+v10*(z(3)/10).^0.16;

fs4=v(4,:)+v10*(z(4)/10).^0.16;

fs5=v(5,:)+v10*(z(5)/10).^0.16;

fs6=v(6,:)+v10*(z(6)/10).^0.16;

fs7=v(7,:)+v10*(z(7)/10).^0.16;

fs8=v(8,:)+v10*(z(8)/10).^0.16;

fs9=v(9,:)+v10*(z(9)/10).^0.16;

fs10=v(10,:)+v10*(z(10)/10).^0.16;

fs11=v(11,:)+v10*(z(11)/10).^0.16;

fs12=v(12,:)+v10*(z(12)/10).^0.16;

fs13=v(13,:)+v10*(z(13)/10).^0.16;

fs14=v(14,:)+v10*(z(14)/10).^0.16;

fs=[fs1' fs2' fs3' fs4' fs5' fs6' fs7' fs8' fs9' fs10' fs11' fs12' fs13' fs14'];

v1=v(1,:);

v2=v(2,:);%1,3,5点脉动风时程曲线对比

v3=v(3,:);

v5=v(5,:);

v14=v(14,:);

figure(1);

plot(t,v2,'b-');%画出v1时程

xlabel('t/s');

ylabel('v/(m/s)');

axis([0 600 -30 30]);  %x轴的范围0~600，y轴的范围-30~30

figure(2);

N=3;

wpt1=wpdec(v2,N,'db1');

plot(wpt1);           %分解树结构

%以下是基于小波包方法的分析

delta=dt;

b_w=1;

for i=1:2^N%wpcoef(wpt1,[n,i-1])是求第n层第i个节点的系数

E(i,:)=(wpcoef(wpt1,[N,i-1]).^2);%求第i个节点的范数平方，其实也就是平方和

delta_f=b_w*1/delta;%J尺度带宽/scale_a^(j-1)

end

figure(3)

[E_total]=sum(E); %求总能量

P=E_total./8;

plot(P);

legend('分信号平均功率谱');

xlabel('时间/(t/s)');

ylabel('功率谱/(m^2/s)');

%功率谱估计

N2=N;

M=600;

wpt2=wpdec(v2,N2,'db1');

for i=1:2^N2; %wpcoef(wpt1,[n,i-1])是求第n层第i个节点的系数

E2(i,:)=(wpcoef(wpt2,[N2,i-1]).^2);%求第i个节点的范数平方，其实也就是平方和

for j=1:M

stf(i,j)=E2(i,j)^2;

end

end

Gxx=sum(stf)./delta_f;

%显示处理

figure(4);

subplot(111);

tt=0:dt:(M-1)*dt;

ff=10/8:10/8:10;

[X,Y]=meshgrid(tt,ff);%全部

[Z]=stf;%全部

mesh(X,Y,Z);

%plot3(X,Y,Z);

xlabel('time / sec');

ylabel('frequency / Hz');

这是我写的程序，我知道程序一定不对，很乱，希望有大神可以指导一下