谐波叠加法-求风速.rar_kaimal谱_谐波_随机风速_风力发电谐波_风速 matlab

%风速时程参数设置 m=10; N=2^8; dt=0.5; omegaup=2*pi; % 风速谱参数设定 L=1000; z=50; z0=0.03; %地面粗糙度 Uz=40; %50m处的平均风速 delta=100; %模拟点间距 lambda=10; %空间相关函数中的系数 K=0.4; M=2*N; v=zeros(m,M*m); t=0.5*(0:1:(M*m-1)); domega=omegaup/N; D=zeros(m,m,N); U=K*Uz/log(z/z0); %U为摩擦速度 %形成目标谱 omega1=omegaup/N:omegaup/N:omegaup; Sw1=200*U^2.*z/Uz./(1+50.*omega1.*z./(2*pi*Uz)).^(5./3); for j=1:m rand('state',0); thet=2*pi*rand(j,N); for l=1:N %%??? omega(l)=(l-1)*domega+j/m*domega; end Sw=200*U^2.*z/Uz./(1+50.*omega.*z./(2*pi*Uz)).^(5/3); %计算谱数据库矩阵，功率谱计算，kaimal谱 for j1=1:m for l=1:m for k=1:N Coh(j1,l,k)=(exp(-lambda*omega(k)*delta/(2*pi*Uz)))^(abs(j1-l)); S(j1,l,k)=Sw(k)*Coh(j1,l,k); end end end %Cholesky分解 for i=1:1:N H(:,:,i)=chol(S(:,:,i)); H(:,:,i)=H(:,:,i)'; end %填充谱数据矩阵D D(:,j,:)=H(:,j,:); i=sqrt(-1); B1=sqrt(2*domega).*D(j,:,:); for ii=1:j for jj=1:N B2(ii,jj)=B1(1,ii,jj); end end B2=B2.*exp(i.*thet); for jj=1:j G(jj,1:M)=fft(B2(jj,:),M); for jjj=2:m G(jj,((jjj-1)*M+1):(jjj*M))=G(jj,1:M); end end %谐波叠加生成模拟点的风速时程 for p=1:M*m for k=1:j v(j,p)=v(j,p)+real(G(k,p)*exp(i.*k./m.*domega.*(p-1).*dt)); end end end % 显示风速时程 figure(1) plot(t,v(1,:)) figure(2) plot(t,v(5,:)) [power1,freq1]=psd(v(1,:),M*m,2,boxcar(512),0,'mean'); [power5,freq5]=psd(v(5,:),M*m,2,boxcar(512),0,'mean'); [power15,freq15]=cpsd(v(1,:),v(5,:),[],512,512,2); %计算第一、第五点的模拟互功率谱 Sw15=[]; for i=1:N s15=S(1,5,i); Sw15=[Sw15,s15] end %功率谱检验 figure(3) loglog(freq1,power1,'r',omega,Sw1,'b') figure(4) loglog(freq5,power5,'r',omega,Sw1,'b') figure(5) loglog(freq15,abs(power15),'r',omega,Sw15,'b')