脉动风功率谱源码

 for i=1:63
 Y=fft（udm（:bb（i））8192）/6000;
 F=5*linspace（018192/2+1）;
 Ap（i:）=2*abs（Y（1:8192/2+1））;
 %plot（F2*abs（Y（1:8192/2+1）））;
 %title（‘No.14样本Fourier幅值谱（NFFT=8192）‘）;
 %xlabel（‘频率/Hz‘）;
 %ylabel（‘幅值/（m\cdots^{-1}）‘）;
 Sn（i:）=2*（abs（fft（udm（:bb（i））8192）））.^2/6000/10;
 %plot（FSn（1:8192/2+1）‘r‘）;
 %title（‘No.14样本Fourier功率谱（NFFT=8192）‘）;
 %xlabel（‘频率/Hz‘）;
 %ylabel（‘幅值/（m^2\cdots^{-1}）‘）;
 %di=diff（F）;
 intSn（1i）=sum（Sn（i2:8192/2+1）.*diff（F））;
 sigmafang（1i）=mean（udm（:bb（i））.^2）;
 Luu（1i）=mm（bb（i））./sigmafang（1i）.*1/2*sum（xcorr（udm（:bb（i））‘biased‘））;%求湍流积分尺度用无偏估计子。
 end
 %% 归一化相关函数
  AC=zeros（600163）;
 for i=1:63
 AC（:i）=xcorr（udm（:bb（i））3000‘coeff‘）;
 end
 %% 确定积分区间（0.05）
  for i=1:63
     for k=3001:6001
         if AC（ki）<0.05
             sp（i1）=k;
             break
         end
     end
  end
 %% 相关函数积分
  for i=1:63
      intlag（i1）=sum（AC（3001:sp（i1）））;
  end
 %% 湍流积分尺度（有偏估计）
  for i=1:63
     Luu（1i）=mm（bb（i））./sigmafang（1i）.*intlag（i1）‘;
  end
%% 画图脉动风速谱+Von Karman谱
 for k=51:63
    Psd=Sn（k1:8192/2+1）.*F./sigmafang（1k）;
    figure;
    loglog（F.*Luu（1k）./mm（:bb（k））Psd）;
    title（[‘No.‘num2str（bb（k））‘样本顺风向脉动风速谱（NFFT=8192）‘]）;
    xlabel（‘频率/Hz‘）;
    ylabel（‘\itnS_u（n）/\sigma^2‘）;
    hold on
    loglog（F.*Luu（1k）./mm（:bb（k））4*F*Luu（1k）/mm（:bb（k））./（1+70.8*（F*Luu（1k）/mm（:bb（k）））.^2）.^（5/6）‘r‘）;
    legend（‘实测风速谱‘‘Von Karman谱‘）;
    hold off;
 end
a*（x*1398.38/17.8866）.^（c*d-2/3）./（1+b*（x*1398.38/17.8866）.^c）.^d
hold on
loglog（F14003*（F*Luu（11）/mm（:bb（1）））.^21.0487./（1+102605*（F*Luu（11）/mm（:bb（1）））.^4.9146）.^4.4186