cyclic_cross_periodogram.m

matlab编写基于循环平稳谱估计的认知无线电频谱感程序。

function S=cyclic_cross_periodogram(x,y,N,a,g,L)
% 
% CYCLIC_CROSS_PERIODOGRAM
% 
% 	       calculate the spectral correlation density function
%              of signals x and y using the cyclic cross periodogram
%              algorithm
%              
%              Input parameters are:
%              x,y signals
%              N   length of time window used for estimating frequency 
%                  segments
%              a   window used for smoothing segments
%              g   window for smoothing correlation
%              L   decimation factor
%              
% USAGE        S=cyclic_cross_periodogram(x,y,N,a,g,L)
%
%              e.g s=cyclic_cross_periodogram(s1,s1,64,'hamming','hamming',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        %这边的x，y的频谱中频率范围是【-fs，fs】
  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         %对x，y的谱分别进行平滑,i相当于时间，在一个周期内对其积分
    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,floor(N/2)+1);
for alpha=0:floor(N/2)        %alpha的范围是【0 fs】，那么alpha/2的范围是【0，fs/2】,这边的alpha可以理解成实际中的alpha/2，因此范围是到N/2
for f=-floor(N/2):floor(N/2)  %f的范围是【-fs/2,fs/2】
    f1=f+alpha;
    f2=f-alpha;
    if ((abs(f1)<N/2)&(abs(f2)<N/2))
       S(f+floor(N/2)+1,alpha+1)=sum(g.*X(f1+N+1,:).*Y(f2+N+1,:));
    end
end
end
S=abs(S);