clear
close('all')
U=load('111.txt');   %读入实测的风速数据
u=U(1:2048,1)-mean(U(1:2048,1));

Fs=10;                %采样频率
N=2^nextpow2(length(u));   %变换长度

z=5;                  %测点高度
v0=27.3;              %该高度处的风速大小
rms=std(u);           %风速均方根
sdz=rms^2;            %风速标准差

[Sv,f]=pmtm(u,4,N,Fs); %利用多窗口法实现频谱估计，本例直接调用pmtm函数进行估计,格式[Pxx[,f]]=pmtm(x[,nw,Nfft,Fs])x为信号序列；nw为时间带宽积，缺省值为4。通常可取2，5/2，3，7/2；Nfft为FFT长度；Fs为采样频率。
fn=z*f/v0;             %归一化功率      
Sv=Sv.*f/sdz;          %归一化脉动风速功率谱

v10=v0*(10/z)^0.16;
x=1200*f/v10;
Sv_Davenport=4*x.^2./(6*(1+x.^2).^(4/3));    �venport风速谱

x1=z*f/(v10*(z/10)^2);
Sv_Simiu=200*x1./(6*(1+50*x1).^(5/3));       %Simiu风速谱

Sv_Kaimal=200*fn./(6*(1+50*fn).^(5/3));       %Kaimal风速谱

x2=1800*f/v10;
Sv_Harris=4*x2./(6.677*(2+x2.^2).^(5/6));     %Harris风速谱

Lu=100*(z/30)^(1/2);
x3=Lu*f/v0;
Sv_Karman=4*x3./((1+70.8*x3.^2).^(5/6));       %Karman风速谱 

figure(2)
loglog(fn,Sv,'k',fn,Sv_Davenport,'r',fn,Sv_Simiu,'b',fn,Sv_Kaimal,'c',fn,Sv_Harris,'m',fn,Sv_Karman,'g');
xlabel('nz/u')
ylabel('nSv(n)/{\sigma}^{2}')
legend('Wind Tunnel','Davenport','Simiu','Kaimal','Harris','Karman')
http://s15/mw690/002XKryigy70b8geZjwfe&690

111.txt