function [FF,CS,FF_S, CS_S] = Cal_CSPD(UU,VV,Fa)
%% matlab 计算互谱子程序
%% 输入量：UU,VV 分别为输入的两个待计算的随机序列，Fa 为采用频率
%% 输出量：FF,CS 为未平滑的频率和互谱，FF_S,CS_S 为结果平滑的频率和相应的互谱
%% 可以验证 UU 和 VV 的方差等于 CS 积分
CS=[];
CS_S=[];
Row=length(UU);
if(Row~=length(VV))
return;
end
X=fft(UU);
Y=fft(VV);
CS=real(X.*conj(Y))/Row*2/Fa;%%
FF=[0:Row-1]/Row*Fa;
FF_S=[];
CS_S=[];
NUM_SS=0;
ii=1;
kk=0;
nn=0;
jj=0;
while ii<Row/2&jj<Row/2&kk<Row/2 &nn<Row/2
nn=floor(ii+ii/20);
CS_S(NUM_SS+1)=0.0;
FF_S(NUM_SS+1)=0.0;
kk=0;
for jj=ii:nn
if(jj>=Row/2)
break;
end
CS_S(NUM_SS+1)=CS_S(NUM_SS+1)*kk/(kk+1.0)+CS(jj)/(kk+1.0);
FF_S(NUM_SS+1)=FF_S(NUM_SS+1)*kk/(kk+1.0)+FF(jj)/(kk+1.0);
kk=kk+1;
end
NUM_SS=NUM_SS+1;
ii=floor(ii+1+ii/20);
end
评论6