cyclostationary的matlab代码

以下是一个简单的MATLAB代码实现，可以计算信号的循环自相关函数和循环谱密度。假设信号为x，采样频率为fs。

循环自相关函数：

function [R] = cyclo_corr(x,tau,M)
% x: 输入信号
% tau: 循环时延
% M: 循环块大小
L = length(x);
R = 0;
for n = 1:M
    for m = 1:M
        R = R + x(m)*conj(x(mod(n+m-2,L)+1))*...
            exp(-1j*2*pi*(m-1)*(n-1)/M);
    end
end
R = R/M;

循环谱密度:

function [S,f] = cyclo_spec(x,M,fs)
% x: 输入信号
% M: 循环块大小
% fs: 采样频率
L = length(x);
f = linspace(-fs/2,fs/2,2*M+1);
S = zeros(size(f));
for n = 1:M
    R = cyclo_corr(x,n,M);
    S = S + abs(fftshift(fft(R,2*M+1))).^2;
end
S = S/M^2/L/fs;

此处的循环自相关函数和循环谱密度的实现是基于块循环的方法。循环谱密度的计算实际上是对循环自相关函数的傅里叶变换，然后进行平均。