[moneylancy]

关于循环谱
大家好！！

我最近在做调制信号识别，其中需要用到循环谱。

我用了一篇论文里的程序，但是运行却有问题，还请大家指教。

这篇论文的名字是《Detection and identification of cyclostationary signals》。
程序如下：
function [Sx,alphao,fo]=autofam(x,fs,df,dalpha)

if nargin > 4 | nargin < 4
error('Wrong number of arguments.');
end

% Definition of Parameters %
Np = pow2(nextpow2(fs/df));
L = Np/4;
P = pow2(nextpow2(fs/dalpha/L));
N = P*L;

% Input Channelization %
if length(x) < N
x(N) = 0;
elseif length(x) > N
x = x(1:N);
end

NN = (P-1)*L+Np;
xx = x;
xx(NN) = 0;
xx = xx(;
X = zeros(Np,P);

for k=0-1
X(:,k+1) = xx(k*L+1:k*L+Np); end

% Windowing %
a = hamming(Np);
XW = diag(a)*X; %每段加窗
% XW = X; %I think this sentence should be kicked out.

% First FFT %
XF1 = fft(XW);
XF1 = fftshift(XF1);
XF1 = [XF1(:,P/2+1) XF1(:,1/2)];

% Downconversion %
E = zeros(Np,P);

for k = -Np/2:Np/2-1
for m = 0-1
E(k+Np/2+1,m+1) = exp(-i*2*pi*k*m*L/Np);
end
end

XD = XF1.*E;
XD = conj(XD'); %XD转置的复共轭

% Multiplication %
XM = zeros(P,Np^2);

for k = 1:Np
for q = 1:Np
XM(:,(k-1)*Np+q) = (XD(:,k).*conj(XD(:,q)));
end
end

% Second FFT %
XF2 = fft(XM);
XF2 = fftshift(XF2);
XF2 = [XF2(:,Np^2/2+1:Np^2) XF2(:,1:Np^2/2)];
XF2 = XF2(P/4:3*P/4,;
M = abs(XF2); %Absolute value and complex magnitude
alphao = (-1:1/N:1)*fs;
fo = (-.5:1/Np:.5)*fs;
Sx = zeros(Np+1,2*N+1);

for k1 = 1/2+1
for k2 = 1:Np^2
if rem(k2,Np) == 0
q = Np/2-1; %increase 1
else
q = rem(k2,Np)-Np/2-1; %increase 1
end
k = ceil(k2/Np)-Np/2-1; %increase 1
p = k1-P/4-1;
alpha = (k-q)/Np+(p-1)/L/P;
f = (k+q)/2/Np;
if alpha < -1 | alpha > 1
k2 = k2+1;
elseif f < -.5 | f > .5
k2 = k2+1;
elseif rem(k+q,2)==0
kk = 1+Np*(f+.5);
ll = 1+N*(alpha+1);
Sx(kk,ll) = M(k1,k2); %出现问题处，这里计算出的kk有时为非整数。为什么？？M到Sx的对应是怎样的？？
end
end
end