风谱到时程的模拟.txt

采用自回归滑动模型进行风速时程的模拟

%-----------------------------------------
% Simultate a multivariate ststionary vector
% process using the spectral representatio method
% based on the cross-spectral density matrix.
% The vector process of wind speeds over a tall building
% is used as a specific example.
%-----------------------------------------
%Note:
% [Pxy,f] = csd(x,y,nfft,Fs) : Estimate the cross spectral density (CSD) of two signals.
% [Pxx,freq] = pburg(x,p,nfft,Fs) : Power spectrum estimate using the Burg method.
% or other methods 
% c = xcorr(x,y,maxlags,'option')
% c = xcorr(x,maxlags,'option'), 'option' = 'biased'
% Cross-correlation function estimate  
%
function [H,Uz,st,vt,dt,M,f]=simu2(filename)
tic;
fr = fopen('windzbeam.txt', 'rt');
% Property of the wind
ustar = fscanf(fr,'%f',[1,1]);%1.76
z0 = fscanf(fr,'%f',[1,1]);%0.001266
NV = fscanf(fr,'%i',[1,1]);%3
% NV = number of variables, or number of heights
for i=1:NV
   H(i) = fscanf(fr,'%f',[1,1]);
end
H
% heights for each variable, ordered from low to high.
% Data controling the simulation
wu = fscanf(fr,'%f',[1,1]); % the upper cutoff frequency  %6.28
N = fscanf(fr,'%i',[1,1]); % number of intervals of [0,wu] to be divided. %2048
M = fscanf(fr,'%i',[1,1]); % values at NV*M time instants will be simulated;%4096
% Condition: M>=2N
% then by formulation
dw = wu/N; % length of frequency interval%6.28/2048=0.003066
dt = 2*pi/(M*dw); % length of time intervals%0.5
Fs=1/dt;%2
%
nfft = ceil(log2((4*pi*Fs)/dw));
nfft = 2^(nfft-1);
%nfft = 256;
%
T0 = NV*(2*pi/dw);
% period of the simulated process, sample function will usually be simulated for one period
%
NT = fscanf(fr,'%i',[1,1]); % Number of sample points to be calculated and output%8192
VNT = fscanf(fr,'%i',[1,1]); % Number of sample points to be used for verification%8192
if VNT < nfft
   fprintf('VNT < fft\n');
   quit;
end
% Calculate the mean wind speed
H=H(:);
Uz=Uz_f(z0,ustar,H);
Uz=Uz(:);
rc=RC(Uz,H);
HUz=H./Uz;
SC1 = sqrt(15.91549430918953*ustar*ustar)*sqrt(HUz);
SC2 = 7.957747154594*HUz;
%
B = B_f(SC1,SC2,rc,NV,H,wu,N,M,dw,dt);
% FFT
for m=1:NV
   fprintf('Computing g, Step No. %i of %i \n ',m,NV);
   for j=m:NV
      temp = B((j-1)*j/2+m,:);
      B((j-1)*j/2+m,:)=fft(temp);
   end
end
% Now the sample function
for j=1:NV
   fprintf('Computing f, Step No. %i of %i \n ',j,NV);
   for p=1:NT
      sum = 0.0;
      q=mod(p-1,M);
      for m=1:j
         angle3 = m*dw*p*dt/NV;
         whole = B((j-1)*j/2+m,q+1)*(cos(angle3)+1.0i*sin(angle3));
         sum = sum + real(whole);
      end
      f(j,p)=sum;
   end
end
toc;
st=toc;
tic;
% Now compare to verify
NP = fscanf(fr,'%i',[1,1]); % Number of variable pairs to be verified
nplot=0;
for i=1:NP
   fprintf('Verifying, Step No. %i of %i \n ',i,NP);
   NOV = fscanf(fr,'%i',[2,1]); % No. of first and second variable
   % Generate the sample functions
   j=0;
   for j1 = 1:2
      if j == NOV(j1)
         break;
      end
      j=NOV(j1);
      for p=1:VNT
         sum = 0.0;
         q=mod(p-1,M);
         for m=1:j
            angle2 = m*dw*p*dt/NV;
            whole = B((j-1)*j/2+m,q+1)*(cos(angle2)+1.0i*sin(angle2));
            sum = sum + real(whole);
         end
         fx(j1,p)=sum;
      end
   end
   % The spectrum and the auto/cross-correlation function
   % Estimated from the simulated sample functions
   maxlags = 500;
   if NOV(1) == NOV(2)
      [Pxx,freq] = pburg(fx(1,:),4,nfft,Fs);
      Pxx=Pxx*0.5;
      [Pxx1,freq1]=psd(fx(1,:),2048*2,1/0.5,boxcar(1024),0,'mean');
      Pxx1 = Pxx1*0.5;%规一化修正
      %[Pxx,freq] = pwelch(fx(1,:),nfft,Fs,[],nfft/2);
      %Pxx=Pxx/2;
      cx = xcorr(fx(1,:),fx(1,:),maxlags,'biased');
   else
      %[Pxx,freq] = csd(fx(1,:),fx(2,:),256,Fs);
      [Pxx,freq] = csd(fx(1,:),fx(2,:),nfft,Fs,[],nfft/2);
      cx = xcorr(fx(1,:),fx(2,:),maxlags,'biased');
   end
   for loop = 1:2*maxlags+1
      cx_t(loop) = -maxlags*dt+(loop-1)*dt;
   end
   %The target psd and auto/cross-correlation function
   for loop = 1:N
      wx = (loop-1)*dw;
      psd_w(loop) = wx;
      %S0w = S0w_f(ustar,z0,NV,H,wx);
      S0w = S0w_f_1(SC1,SC2,rc,wx);
      psd(loop) = S0w(NOV(1),NOV(2));
   end
   Rx = ifft(psd);
   Rx = real(Rx);
   for loop = 1:maxlags+1
      R_t(loop) = -(maxlags/2)*2*dt+(loop-1)*2*dt;
      po = abs(loop-maxlags/2-1)+1;
      R(loop) = 2*wu*Rx(po);
   end
   % Now plot
   nplot=nplot+1;
  % figure(nplot);
   %loglog(freq*2*pi,Pxx/2/pi,'b',psd_w,psd,'k',freq*2*pi,Pxx1/2/pi,'r');
%   loglog(freq*2*pi,Pxx/2/pi,'b',psd_w,psd,'k');
 %  nplot=nplot+1;
 %  figure(nplot);
 %  plot(cx_t,cx,R_t,R);
   fname = strcat('SRM_',int2str(nplot));
   save(fname,'f','fx','freq','Pxx','Pxx1','psd_w','psd','cx_t','cx','R_t','R','Uz');
   %
end
%
nfft
fclose(fr);
toc;
vt=toc;
%-----------------------------------------
% Compute the B matrix, B(NV,NV,M) 
% But B(j,m,l) with m<=j and 1<= l <=M will be saved as B((j-1)*j/2+m,l)
%-----------------------------------------
function B = B_f(SC1,SC2,rc,NV,H,wu,N,M,dw,dt)
B=zeros(NV*(NV+1)/2,M);
sdw=sqrt(dw);
twopi=2*pi;
% The random phase angle
FI = twopi*rand(NV,N);
%
for l=1:N
   %tic
   fprintf('Computing B, Step No. %i of %i \n ',l,N);
   for m=1:NV
      wlm = (l-1)*dw+m*dw/NV;
      S0wlm = S0w_f_1(SC1,SC2,rc,wlm); 

      Hwlm = chol(S0wlm); % Hwlm'*Hwlm = S0wlm
      lHl = abs(Hwlm);        % magnitude
      theta = angle(Hwlm);  % phase angle
      for j=m:NV
         lBl=2*lHl(m,j)*sdw;
         angle1 = -theta(m,j)+FI(m,l);
         realx = lBl*cos(angle1);
         imaginex = lBl*sin(angle1);
         B((j-1)*j/2+m,l) = realx+imaginex*1.0i;
      end
   end
   %toc
end

%-----------------------------------------
% The Davenport coherence function
% between the velocity fluctuations
% at different heights
%-----------------------------------------
function rdzw = rdzw_f(dz,w,Uz1,Uz2)
% dz = abs(z1-z2), distance between the two heights z1 and z2
% w = vector of cicular frequencies
% Uz1 = mean velocity at z1
% Uz2 = mean velocity at z2
Cz = 10; % a constant
ndz = length(dz);
nw = length(w);
for i=1:ndz
   for j=1:nw
      rdzw(i,j) = exp(-0.31830988618379*w(j)*Cz*dz(i)/(Uz1(i)+Uz2(i)));
   end
end
%
%-----------------------------------------
% The Kaimal (two-sided) power spectral density function
% of the longitudinal wind velocity fluctuations at some heights
%-----------------------------------------
function Szw=Szw_f(ustar,z,w,Uz)
% z = vector of heights
% w = vector of cicular frequencies
% Uz = vector of mean velocities at z
nz = length(z);
nw = length(w);
const = 15.91549430918953*ustar*ustar;
% = (1/2)*(200/(2*pi))*ustar^2;
for i = 1:nz
   z_Uz(i) = z(i)/Uz(i);
end
%
for i=1:nz
   for j=1:nw
      Szw(i,j) = const*z_Uz(i)/(1+7.957747154594*w(j)*z_Uz(i))^1.6667;
      % 7.957747154594 = 50/(2*pi)
   end
end
%
%-----------------------------------------
% The logarithm law for mean velocity profile
%-----------------------------------------
function Uz=Uz_f(z0,ustar,z)
%Uz = mean velocity at height z;
%z0 = the surface roughness length; 
%ustar = shear velocity of the flow
%z = vector of height;
k = 0.4; % Kalman constant;
n = length(z);
for i=1:n
   Uz(i) = ustar*log(z(i)/z0)/k;
end
%
%-----------------------------------------
% Calculate the constant matrix for the coherence function
%-----------------------------------------
function rc=RC(Uz,H)
Cz = 10; % a constant
NV=length(H);
rc=zeros(NV,NV);
for i=1:NV
   for j=i:NV
      dz=abs(H(i)-H(j));
      rc(i,j)=exp(-0.31830988618379*Cz*dz/(Uz(i)+Uz(j)));
      rc(j,i)=rc(i,j);
   end
end
%-----------------------------------------
% The Kaimal (two-sided) power spectral density function
% of the longitudinal wind velocity fluctuations at some heights
%-----------------------------------------
function Szw=Szw_f_1(SC1,SC2,w)
TEMP=1+w*SC2;
TEMP=TEMP.^0.83333333333333;
Szw=SC1./TEMP;
%
%-----------------------------------------
% Compute the cross-spectral density matrix S0w(NV,NV)
% for some w
%-----------------------------------------
function S0w = S0w_f_1(SC1,SC2,rc,w)
% Diagonal Components, auto-spectral
Szw=Szw_f_1(SC1,SC2,w);
TEMP=rc.^w;
S0w=(Szw*Szw').*TEMP;