一种简易的频偏估计补偿方法，使用QAM等信号。估计精度受FFT长度限制，可以作为粗频偏估计。

Nfft = 1024;                        % FFT长度
N = 10*Nfft;                        % 仿真符号数
M = 16;                             % 调制QAM16
freq = 1e3;                         % 频偏
fs = 1e5;                           % 采样率

msg = randi([0 M-1],N,1);
iq = qammod(msg,M);
h = rcosdesign(0.2,8,4,'sqrt');     % 成型滤波器
tx = conv(upsample(iq,4),h,'same'); % 滤波

ch = awgn(tx,20,'measured');        % 过信道

t = (0:4*N-1)/fs;
car_offset = exp(1i*2*pi*freq*t');  % 产生频偏用于载波调制
rx0 = ch.*car_offset;                % 频偏
rx = conv(rx0,h,'same');

tstart = 0;
out = zeros(4*N,1);
freq_est = zeros(N/Nfft,1);
for i=1:4*N/Nfft
    spec = fftshift(abs(fft(rx(i*Nfft-Nfft+1:i*Nfft).^4)));  % 求4阶频谱
    [~,maxId] = max(spec);                                   % 求最大频谱分量位置
    offsetIdx = maxId - Nfft/2-1;                            % 以中心频点为参考的位置
    estFreq = offsetIdx*fs/Nfft/4;                           % 频偏计算
    freq_est(i) = estFreq;

    tnew = (0:Nfft-1)/fs+tstart;
    car_compl = exp(-1i*2*pi*estFreq*tnew');                  % 补偿载波
    out(i*Nfft-Nfft+1:i*Nfft) = rx(i*Nfft-Nfft+1:i*Nfft).*car_compl; % 补偿
    tstart = tstart+N/fs;
end

scatterplot(rx0(1:4:end));
title('Channel Constellation');
scatterplot(out(1:4:end));
title('FreqOffsetComp Constellation');
periodogram([rx0,out],[],1024,1e5,'center');