循环谱密度函数MATLAB程序

%% Initialization
% 参考文献
% E. L. Da Costa “Detection and Identification of
% Cyclostationary Signals“. MS Thesis. 1996. 
% 公式（9）使用对称形式时，本程序参考公式（8），所以不再引入频移

clc;clear;close all;
%% 系统参数配置
span = 6;               % Filter span
fs = 5e2;
sps = 8;                % Samples per symbol
rs = fs/sps;
fc = 1*rs;
m_ord = 2;              % Size of the signal constellation
m_bit = log2（m_ord）;    % Number of bits per symbol

EsNo = 100;
burst_length = 100;
rolloff = 0.35;         % Filter rolloff
% Generate a burst data symbols using the randi function.
%% 生成测试数据
rrcFilter = rcosdesign（rolloff span sps ‘sqrt‘）; 
data = randi（[0 m_ord-1] burst_length 1）;
% PSK modulation
modData = pskmod（data m_ord pi/4）;
% modData = qammod（data m_ord）;
% Using the upfirdn function upsample and filter the input data.
txSig = upfirdn（modData rrcFilter sps）;
% Convert the Eb/N0 to SNR and then pass the signal through an AWGN channel.
EbNo = EsNo - 10*log10（m_bit）;
SNR = EbNo + 10*log10（m_bit） - 10*log10（sps）;
AwgnSig = txSig（1:end） + awgn（txSig（1:end） SNR ‘measured‘）; 
% AwgnSig = txSig（1:end-3） + awgn（txSig（1:end-3） SNR ‘measured‘）; %截取非整数周期，让DFT泄露
t = [0:1/fs:（length（AwgnSig）-1）/fs]‘;
% rxSig = real（AwgnSig.*exp（j*2*pi*fc*t））; %加残留载波
rxSig =