用四阶累积量估计信号载频

EccoC

这是我产生的16QAM信号
function [s] = singnal()
%% 参数定义
NUM = 4*1000;  
SYM_LEN = NUM/4;                                                              % 符号个数
fd = 151200;                                                                  % 码元速率                                                            % 采样率
fs_send = 1e8;                                                                % 采样率
fc = 1.49e7;                                                                  % 载波频率
fc_deviation = 1.49e7-100;                                                    % 频偏
fc_deviation1=1.49e7-200;
SNR = -30;                                                                     % 信噪比
alpha = 0.25;                                                                 % 滚降系数
INSERT_TIMES = 661;                                                           % 内插长度
delay = 3;                                                                    % 滤波器延迟
N = 2097152*2;                                                                % 傅里叶变换点数

%% 产生信源和I,Q分量,第一个信号
message = randi([0 1],1,NUM);
message(message<=0.5)=0;
message(message>0.5)=1;
signal_base_band1I=message(1:4:end);
signal_base_band2I=message(2:4:end);      
signal_base_band3Q=message(3:4:end);
signal_base_band4Q=message(4:4:end);
signal_base_bandI=zeros(1,SYM_LEN);
signal_base_bandQ=zeros(1,SYM_LEN);
for i=1:1:SYM_LEN
    tempI=[signal_base_band1I(i),signal_base_band2I(i)];
    if     tempI==[1,0]
            signal_base_bandI(i)=1;
    elseif tempI==[1,1]
            signal_base_bandI(i)=3;
    elseif tempI==[0,1]
            signal_base_bandI(i)=-3;
    else
            signal_base_bandI(i)=-1;
    end;
    tempQ=[signal_base_band3Q(i),signal_base_band4Q(i)];
    if     tempQ==[1,0]
            signal_base_bandQ(i)=1;
    elseif tempQ==[1,1]
            signal_base_bandQ(i)=3;
    elseif tempQ==[0,1]
            signal_base_bandQ(i)=-3;
    else
            signal_base_bandQ(i)=-1;
         
    end;
    signal_base_band(i)=signal_base_bandI(i)+signal_base_bandQ(i)*1i;
end;
%% 匹配滤波
% 插值
signal_after_insertI=zeros(1,SYM_LEN*INSERT_TIMES);
signal_after_insertI(1:INSERT_TIMES:end)=signal_base_bandI;
signal_after_insertQ=zeros(1,SYM_LEN*INSERT_TIMES);
signal_after_insertQ(1:INSERT_TIMES:end)=signal_base_bandQ;
% % 成型滤波
signal_after_fiterI2 = rcosflt(signal_after_insertI,1,661,'fir/fs/sqrt',alpha,delay)';
signal_after_fiterQ2 = rcosflt(signal_after_insertQ,1,661,'fir/fs/sqrt',alpha,delay)';
signal_sendI = signal_after_fiterI2(INSERT_TIMES*delay+1:end-INSERT_TIMES*delay);
signal_sendQ = signal_after_fiterQ2(INSERT_TIMES*delay+1:end-INSERT_TIMES*delay);
%% 加载波
t = 0:1/fs_send:(661000-1)/fs_send;
%% 加噪声
s = signal_sendI.*cos(2*pi*fc_deviation1*t)+signal_sendQ.*sin(2*pi*fc_deviation1*t); %变为复数信号
s = awgn(s,SNR,'measured'); %加噪声

end


下面是我自己编写的四阶累积量的程序
%利用四阶累积量来估算载频%
y = singnal();
%四阶累积量的求法应该要与载波频率挂上钩，因此不能用求均值的方法来计算%
N = 2097152*2;%进行傅里叶变换的点数
fs = 1e8;%信号采样频率
fc = 1.49e7;%信号的载波频率
NUM = length(y);
M4_4 = [];
M2_2 = [];
k = 1;%用于存放所有的循环矩
%求循环矩
for i = fc-300:1:fc-100
    t = 0:1/fs:(NUM-1)/fs;
    M4 = y.^4.*exp(-1i*8*pi*i*t);
    M4_4(k) = mean(M4);
    M2 = y.^2.*exp(-1i*4*pi*i*t);
    M2_2(k) = mean(M2);
    k = k+1;   
end

%求四阶累积量%
Cy = M4_4 - 3*M2_2.*M2_2;
% Cy = abs(Cy);
f = fc-300:1:fc-100;
plot(f,abs(Cy));
但是结果不对，请大家指点指点，谢谢了