MSK信号频偏估计Matlab代码

以下是一段基于FFT的MSK信号频偏估计的Matlab代码：

% 生成MSK信号
fs = 1000; % 采样率
f0 = 50; % 载波频率
Tb = 1/10; % 每个比特的时间
fc = 2/Tb; % 调制频率
N = 1024; % 信号长度
t = (0:N-1)/fs; % 时间序列
bits = randi([0 1],1,N/2); % 随机生成比特序列
s = zeros(1,N); % 初始化信号
for i = 1:N/2
    if bits(i) == 0
        s((2*i-2)*fs*fc*Tb+1:(2*i-1)*fs*fc*Tb) = cos(2*pi*f0*t((2*i-2)*fs*fc*Tb+1:(2*i-1)*fs*fc*Tb));
        s((2*i-1)*fs*fc*Tb+1:2*i*fs*fc*Tb) = sin(2*pi*f0*t((2*i-1)*fs*fc*Tb+1:2*i*fs*fc*Tb));
    else
        s((2*i-2)*fs*fc*Tb+1:(2*i-1)*fs*fc*Tb) = sin(2*pi*f0*t((2*i-2)*fs*fc*Tb+1:(2*i-1)*fs*fc*Tb));
        s((2*i-1)*fs*fc*Tb+1:2*i*fs*fc*Tb) = cos(2*pi*f0*t((2*i-1)*fs*fc*Tb+1:2*i*fs*fc*Tb));
    end
end

% 加入频偏
f_offset = 10; % 频偏
s = s .* exp(1j*2*pi*f_offset*t);

% FFT求频率响应
S = fft(s);
f = linspace(-fs/2,fs/2,N);
f_step = fs/N;
f_idx = round((f0+f_offset)/f_step);
phase_diff = angle(S(f_idx+N/2)) - angle(S(N/2+1));
if phase_diff > pi
    phase_diff = phase_diff - 2*pi;
elseif phase_diff < -pi
    phase_diff = phase_diff + 2*pi;
end
freq_offset = phase_diff / (2*pi*Tb); % 频偏估计值
disp(['真实频偏：',num2str(f_offset),'Hz，估计频偏：',num2str(freq_offset),'Hz']);

代码中首先生成了一个长度为1024的MSK信号，其中每个比特的时间为1/10秒，采样率为1000Hz，载波频率为50Hz，调制频率为20Hz。然后在信号中加入了10Hz的频偏。接着使用FFT求出信号的频率响应，并计算出位于载波频率和频偏之和处的频率分量的相位差。最后根据相位差计算出频偏的估计值。