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');

這是我寫的程序，我知道程序一定不對，很亂，希望有大神可以指導一下

總結

以上是生活随笔為你收集整理的matlab小波三维图,matlab小波包变换估计时变功率谱三维图出图和理想不一样的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。