test.m

一些关于谱分析的matlab程序

R=reshape(ones(5,1)*(round(rand(1,1000))*2-1),1,5*1000);
fs=100;t=1/fs:1/fs:1/fs*5*1000;y=sin(2*pi*20*t);
Q=y.*R;
x=Q;y=Q;N=64;a='hamming';g='hamming';L=1;

if ~exist('L')
  L=1;
end

lx=length(x);
if (length(y)~=lx)
  error('Time series x and y must be same length')
end


n=0:floor((lx-N)/L);
ln=length(n);
a=feval(a,N)';
g=feval(g,ln)';
g=g/sum(g);
a=a/sum(a);

X=zeros(2*N+1,ln);
Y=zeros(2*N+1,ln);

Ts=1/N;

for f=-N:N
  xf=x.*exp(-j*2*pi*f*(0:lx-1)*Ts);
  yf=y.*exp(-j*2*pi*f*(0:lx-1)*Ts);
  for i=1:ln
    n_r=n(i)*L+(1:N);
    X(f+N+1,i)=sum(a.*xf(n_r));
    Y(f+N+1,i)=conj(sum(a.*yf(n_r)));	
  end
end

S=zeros(N+1,N+1);
for alpha=-floor(N/2):floor(N/2)
  for f=-floor(N/2):floor(N/2)
    f1=f+alpha;
    f2=f-alpha;
    if ((abs(f1)<N/2)&(abs(f2)<N/2))
      S(f+floor(N/2)+1,floor(N/2)+alpha+1)=sum(g.*X(f1+N+1,:).*Y(f2+N+1,:));
    end
  end
end
S=abs(S);

%谱相关系数
rho=zeros(N,N);
for f=-floor(N/2):floor(N/2)-1;
  for alpha=-floor(N/2):floor(N/2)-1;
    fplus=f+alpha;
    fminus=f-alpha;
    if ((abs(fplus)<N/2)&(abs(fminus)<N/2))
    if (S(fplus+N/2+1,N/2+1)*S(fminus+N/2+1,N/2+1)~=0);
    rho(f+floor(N/2)+1,alpha+floor(N/2)+1)=abs(S(f+N/2+1,floor(N/2)+alpha+1))/sqrt(S(fplus+N/2+1,floor(N/2)+1)*S(fminus+N/2+1,floor(N/2)+1));
    end
    end
  end
end