111.rar_moutht2t_www.pu111.com_脉动风速_风速_风速谱

resultt=fopen(['G:\毕设数据\眼镜蛇\B类风场\测风速\空模型\91.44mm\时程\0002 (Ve)','.ap']); fgetl(resultt); fgetl(resultt); fgetl(resultt); fgetl(resultt); fgetl(resultt); fgetl(resultt); fgetl(resultt); fgetl(resultt); fgetl(resultt); fgetl(resultt); fgetl(resultt); result=(fscanf(resultt,'%f %f',[4,15104]))'; fclose(resultt); sp=result(:,1); usp=mean(sp); rusp=sp-usp; var_u=(std(rusp)^2);%方差 c=xcorr(rusp,'unbiased'); r=abs(c); rm=mean(r); Lux=rm*usp/var_u; z=0.25; %%%%%%%%%%%试验%%%%%%%%%%%% [pxx,f]=pmtm(rusp,4,2^14,500); pxx=pxx.*f/var_u; %%%%%%%%%%karman%%%%%%%%% fL1=f*Lux/usp; pxx1=4.*(fL1)./(1+70.8*fL1.^2).^(5/6); %%%%%%%%kaimal%%%%%%% fL2=f*z/usp; pxx2=200.*(fL2)./(6*(1+50*fL2).^(5/3)); %%%%%%%%%%davenport%%%%%%%%%%% fL3=(1200*f/usp).^2; pxx3=(2.*fL3)./(3*(1+fL3).^(4/3)); %%%%%%归一化频率%%%%%%%% fx=f*z/usp; figure loglog(fx,pxx,'k'); hold on loglog(fx,pxx1,'r--'); hold on loglog(fx,pxx2,'b:'); hold on loglog(fx,pxx3,'m-.'); xlabel('nz/u'); ylabel('nS_v(n)/σ^2'); legend('风洞试验谱','Karman谱','Kaimal谱','Davenport谱');%,'Davenport谱' hold off