1. 先把MSK的原理和项目难点拆开MSKMinimum Shift Keying最小频移键控是通信原理里绕不开的一个调制方式。很多人第一眼看到它觉得简单——不就是频移键控FSK把两个频率拉近一点吗但真正用MATLAB去把调制、解调、误码率整套跑通的时候就会发现连续相位、I/Q两路延迟、抽样判决的极性翻转这些问题一个接一个冒出来。这篇文章就把我从零开始搭MSK仿真链路的过程完整写一遍包含可直接复现的MATLAB代码、关键波形怎么看、误码率怎么数以及几个只有踩过坑才写得出来的注意事项。适合正在做通信课程设计、毕业设计或者想自己动手把理论变成代码的同学。1.1 信号模型它和FSK、QPSK到底是什么关系MSK本质上是一种连续相位FSK两个发送频率之间的间隔被压到最小刚好是1/(2Tb)这里的Tb是一个数据比特的持续时间。这个最小间隔让调制指数等于0.5所以它叫“最小频移键控”。它的信号可以写成s(t) cos(2π·fc·t θ(t))其中θ(t)在每一个比特周期内线性增加或者减少π/2。也就是说发送比特为1时相位在Tb时间内正向走π/2发送比特为0时相位反向走π/2。频率上正好对应两个频率点fc 1/(4Tb)和fc - 1/(4Tb)。既然本质是FSK为什么又要跟QPSK扯上关系因为MSK还有一种非常经典的等效视角它就是一个在Q路数据上延迟了一个比特周期Tb的OQPSK只是每个符号的成形脉冲不是矩形而是持续2Tb的半正弦脉冲。半正弦成形加上I/Q两路错开Tb才保证了相位路径连续、包络恒定。理解这一层后面写MATLAB调制代码就有了明确方向。1.2 为什么调制难点全都集中在I/Q两路“错位”上MSK的正交展开形式可以写成s(t) I(t)·cos(πt / 2Tb)·cos(2π·fc·t) Q(t)·sin(πt / 2Tb)·sin(2π·fc·t)其中I(t)和Q(t)是两路双极性数据流。I路每个符号持续2TbQ路每个符号也持续2Tb但Q路整体比I路滞后一个Tb。这一点不是人为选择的而是MSK连续相位条件本身推导出来的结果。我在帮人调试代码时发现绝大多数“为什么我的MSK相位不连续”的问题最后都出在这个“Q路延迟Tb”上。有人把Q路延迟写成了2Tb有人干脆没做延迟还有人把I路和Q路的加权函数写反了。一个很短的MATLAB仿真前前后后可能因为这一行代码折磨一个晚上。所以这里先把它当重点记住串并转换后I路取奇数位数据Q路取偶数位数据Q路数据序列在送入调制前必须额外错开一个Tb。2. 调制端MATLAB实现从串并转换到相位路径验证2.1 串并转换和I/Q延迟一份可以直接跑的代码我先给出调制端完整代码。代码里用了两个关键思路一是用repmat把每个2Tb符号展开成采样点序列二是Q路通过索引平移实现延迟Tb。用循环也可以但MATLAB里这种矢量化写法更接近工程习惯也更快。clear; clc; close all; % 基础参数 Rb 1000; % 比特速率 1000 bps Tb 1 / Rb; % 比特周期 fs 200e3; % 采样率 200 kHz sps fs / Rb; % 每个比特的采样点数 fc 20e3; % 载波频率 20 kHz numBits 400; % 仿真比特数取偶数 rng(10); % 随机数据转成双极性 data randi([0 1], numBits, 1); bipolar 2 * data - 1; % 串并转换 dI bipolar(1:2:end); % I路奇数索引 dQ bipolar(2:2:end); % Q路偶数索引 M length(dI); % 每路符号数 sps2 2 * sps; % 每个I/Q符号持续2Tb % I路零阶保持每个符号持续2Tb I_wave reshape(repmat(dI, sps2, 1), [], 1); % Q路先做同样的零阶保持 Q_wave_raw reshape(repmat(dQ, sps2, 1), [], 1); % Q路整体延迟Tb头部用第一个Q符号值补 Q_wave [Q_wave_raw(1:sps); Q_wave_raw(1:end-sps)]; % 总时间轴总时长为 numBits * Tb t (0:length(I_wave)-1) / fs; % 注意这里 I_wave 长度 numBits*sps % 但 Q_wave_raw 长度也为 numBits*sps所以 t 可以直接使用 t t(:); % MSK正交调制 msk I_wave .* cos(pi * t / (2 * Tb)) .* cos(2 * pi * fc * t) ... Q_wave .* sin(pi * t / (2 * Tb)) .* sin(2 * pi * fc * t);这段代码生成的msk就是标准的MSK通带信号。有一个地方值得展开说Q_wave [Q_wave_raw(1:sps); Q_wave_raw(1:end-sps)]并不是随便写的。Q_wave_raw中原本第一个Q符号占前2Tb第二Q符号占接下来的2Tb我们要让整个Q序列延迟Tb也就是头部空出来Tb。头部用第一个Q符号补齐是因为Q路第一个符号本身在Tb到3Tb之间会出现前面Tb时间内如果填0会使初始段相位突变。用第一个Q符号值填充能让相位路径从头到尾都连续。为了让读者少走弯路我再补一句如果只是做理论学习很多人会把Q路延迟直接写成Q_wave [zeros(sps,1); Q_wave_raw(1:end-sps)]这样在大多数情况下也能解调出数据但画相位路径时开头会有一小段异常误码率曲线在低信噪比区域也可能出现莫名的抖动。所以填第一个Q符号而不是填零是更讲究的做法。2.2 相位累加视角另一种验证调制正确性的方法除了上面基于I/Q正交调制的做法MSK还有一个更贴近“连续相位FSK”本源的生成方式直接把相位累加出来每个比特周期内相位线性走±π/2。等效复基带可以这样写base_cp zeros(numBits * sps, 1); phase 0; for k 1:numBits phase_inc (pi / 2) * bipolar(k); % 每个比特相位增量 idx (k-1)*sps 1 : k*sps; base_cp(idx) exp(1j * (phase phase_inc * (0:sps-1) / sps)); phase phase phase_inc; end这个base_cp就是MSK的等效复基带幅度恒为1相位在每比特区间内线性变化±π/2。你可以用这个结果去和上一节正交调制的结果做对比把正交调制信号去掉载波、再分别取I/Q分量得到的复基带与base_cp只差一个固定的初始相位偏移。二者描述的是同一个东西只是坐标系不同。我觉得这个“双视角”特别适合写进报告或者用来自查。当你怀疑正交I/Q代码写错时直接跑一下相位累加版本把unwrap后的相位画出来如果斜率不是±π/2每Tb说明调制代码有问题。如果两个版本画出的相位路径趋势一致基本可以确认调制端没毛病。2.3 用图形检查调制结果相位路径、频谱、包络代码跑通之后强烈建议先画三张图再进入解调环节。第一张是I/Q正交调制的相位路径第二张是等效复基带的频谱第三张是通带波形局部放大。% 从通带信号得到复基带解调端也会用到类似的思路 base_rx I_wave .* cos(pi * t / (2 * Tb)) 1j * Q_wave .* sin(pi * t / (2 * Tb)); figure; subplot(3,1,1); plot(t / Tb, unwrap(angle(base_rx)) / pi); xlabel(t / Tb); ylabel(相位 / π); title(MSK相位路径); grid on; xlim([0 20]); subplot(3,1,2); [psd, f] pwelch(base_rx, hann(1024), 512, 1024, fs, centerdc); plot(f / Rb, db(psd / max(psd))); xlabel(f / Rb); ylabel(归一化功率谱/dB); title(MSK等效基带频谱); xlim([-3 3]); grid on; subplot(3,1,3); plot(t / Tb, msk); xlabel(t / Tb); ylabel(幅度); title(MSK通带波形前10个比特); xlim([0 10]); grid on;相位路径画出来后你会看到一条连续的折线每一段斜率都是±0.5π/Tb没有跳变。这个“连续”是MSK最核心的卖点也是它和普通FSK最大的区别。频谱图画出来后主瓣宽度大约在±0.75Rb附近旁瓣衰减很快。频率占用比QPSK的主瓣略宽但因为相位连续带外泄漏更小。这些结论直接写进实验报告或课程设计说明里都是很扎实的素材。通带波形在时间轴上能看到明显的恒定包络效果。如果你画出来的波形包络在某个位置出现凹陷或者毛刺多半是I/Q延迟没做好或者加权函数没写对。不要急着往下做解调先把调制波形调干净。3. 接收端解调实现相干解调的正确打开方式3.1 为什么我不推荐直接在码元中心抽样很多刚接触MSK的人做完调制之后会下意识地按照QPSK的思维去解调正交下变频、低通滤波、在每个符号中心抽样、判决。这个流程放在QPSK上没问题放在MSK上就很容易翻车。原因在于MSK的I/Q支路经过半正弦成形之后每个支路的“符号中心”并不一定是信号幅度的最大值点。I路加权函数是cos(πt / 2Tb)在t0、2Tb、4Tb这些时刻取到±1Q路加权函数是sin(πt / 2Tb)在tTb、3Tb、5Tb这些时刻取到±1。如果你在码元中心抽样抽到的反而可能是零附近的值。再加上半正弦加权会在相邻符号间导致极性翻转抽出来的序列不做处理误码率会稳定在0.5附近。更可靠的做法是“半正弦匹配滤波”也就是把接收到的基带I/Q信号在每个2Tb窗口内与本地半正弦参考波形做相关。这样才能把能量完整积累起来也天然消除了极性翻转问题。这个思路和最优接收机理论是一致的。3.2 匹配滤波解调代码完整可跑的接收端第2章的调制代码生成的是通带信号这里为了加噪声和统计误码率方便我在接收端先用等效基带信号做解调。等效基带直接由I_wave、Q_wave和半正弦加权构成省去载波之后相干接收的逻辑更清楚。% 发射端等效复基带 base_tx I_wave .* cos(pi * t / (2 * Tb)) 1j * Q_wave .* sin(pi * t / (2 * Tb)); % 加入复高斯噪声第4章会详细讲噪声功率换算 snr_dB 10; Eb mean(abs(base_tx).^2) * Tb; % 每比特平均能量 N0 Eb / 10^(snr_dB / 10); sigma sqrt(N0 / 2); noise sigma * (randn(size(base_tx)) 1j * randn(size(base_tx))); rx base_tx noise; % 相干解调取实部、虚部得到I/Q基带 I_soft real(rx); Q_soft imag(rx); % 逐符号匹配滤波 dI_hat zeros(M, 1); dQ_hat zeros(M, 1); for k 1:M % I路第k个符号区间: [2(k-1)Tb, 2kTb] segI ( (2*k-2)*sps 1 ) : ( 2*k*sps ); wI cos(pi * t(segI) / (2 * Tb)); dI_hat(k) I_soft(segI) * wI / length(segI); % Q路第k个符号区间: [(2k-1)Tb, (2k1)Tb] segQ ( (2*k-1)*sps 1 ) : ( (2*k1)*sps ); if segQ(end) length(Q_soft) wQ sin(pi * t(segQ) / (2 * Tb)); dQ_hat(k) Q_soft(segQ) * wQ / length(segQ); else % 最后一个Q符号窗口超出信号长度时丢弃该符号 dQ_hat(k) NaN; end end % 判决 dI_hat dI_hat 0; dQ_hat dQ_hat 0; % 合并成串行比特最后一个Q符号不计入误码 valid_Q ~isnan(dQ_hat); data_hat zeros(numBits, 1); data_hat(1:2:end) dI_hat; data_hat(2:2:end) dQ_hat; % 统计误码 valid_idx find(~( ... (2:2:numBits) numBits isnan(dQ_hat) ... ));这段代码有几处细节需要说明。Q路第k个符号的窗口是从(2k-1)Tb到(2k1)Tb最后一个Q符号会跨出信号总长度所以我在代码里做了越界判断越界时直接把该符号置NaN统计误码时不参与。实际系统里会在帧尾补一个尾符号这里为了演示简单直接丢掉这个比特即可对整体误码率影响可以忽略。还有一点判决门限是0但相关输出值的正负完全可以代表数据极性不需要再额外做差分编码或者极性翻转处理这正是匹配滤波相比直接抽样的优势。3.3 通带版接收和等效基带版接收的关系可能有读者会问上面一直在操作等效复基带没有真正出现载波fc这样算不算“通带信号解调”从数学本质上讲MSK等效基带已经包含了所有调制信息载波只负责把频谱搬到高频。只要本地载波频率和相位都正确通带信号经正交下变频、低通滤波之后得到的恰好就是等效基带的实部和虚部。也就是说如果非要在MATLAB里走一遍完整通带流程代码也不复杂% 先乘以本地载波 r_I_tmp 2 * msk .* cos(2 * pi * fc * t); r_Q_tmp 2 * msk .* sin(2 * pi * fc * t); % 低通滤波 b fir1(128, 1.2 * Rb / (fs / 2)); I_lp filter(b, 1, r_I_tmp); Q_lp filter(b, 1, r_Q_tmp); % 补偿滤波群延迟后再做匹配滤波 delay (length(b) - 1) / 2; I_soft I_lp(delay1:end); Q_soft Q_lp(delay1:end);这一步在演示完整通信链路时很有意义能让学生看到载波恢复、下变频、滤波这些真实模块。但在做误码率蒙特卡洛仿真时直接用等效基带更快也不会因为滤波器设计不当引入额外误差。实际项目里如果发现通带版本的BER曲线比理论差很多首先要检查的就是低通滤波器的延迟补偿。这个坑我在后面“常见问题”里还会再讲。4. 误码率仿真与结果验证4.1 加噪声的换算方式为什么直接给复基带加噪声做误码率仿真最麻烦的往往不是解调算法本身而是噪声功率怎么设置。如果对通带实信号直接加高斯白噪声噪声带宽和信号带宽不一致换算Eb/N0很容易出错。我的建议是在做BER蒙特卡洛仿真时直接用等效复基带加复高斯噪声。对于复基带信号噪声方差和Eb/N0的关系是sigma² N0 / 2其中N0 Eb / 10^(EbN0dB / 10)Eb是每个比特平均能量sigma²是复噪声每维的方差。所以生成噪声时用noise sqrt(N0/2) * (randn(size(base_tx)) 1j * randn(size(base_tx)));这个公式在BPSK、QPSK、MSK这类正交调制里都通用因为它们的理论误码率最后都能统一到BPSK等价模型上。用等效基带加噪声既保留了调制解调的核心过程又避开了带通噪声带宽换算这个最容易翻车的环节。4.2 BER统计代码框架与理论曲线对比下面是完整的误码率仿真核心代码。我没有把整个文件贴出来只保留最核心的循环方便你移植到自己的项目里。EbN0_dB -2:2:10; ber zeros(size(EbN0_dB)); for idx 1:length(EbN0_dB) EbN0 10^(EbN0_dB(idx) / 10); N0 Eb / EbN0; sigma sqrt(N0 / 2); numErr 0; numBitsStat 0; for trial 1:20 % 每次重新生成随机数据 data randi([0 1], numBits, 1); bipolar 2 * data - 1; dI bipolar(1:2:end); dQ bipolar(2:2:end); M length(dI); I_wave reshape(repmat(dI, sps2, 1), [], 1); Q_wave_raw reshape(repmat(dQ, sps2, 1), [], 1); Q_wave [Q_wave_raw(1:sps); Q_wave_raw(1:end-sps)]; base_tx I_wave .* cos(pi * t / (2 * Tb)) ... 1j * Q_wave .* sin(pi * t / (2 * Tb)); rx base_tx sigma * (randn(size(base_tx)) 1j * randn(size(base_tx))); % 解调匹配滤波参考第3.2节代码 % 这里简写实际运行时替换为完整匹配滤波循环 % 统计误码 numErr numErr sum(data_valid ~ data_hat_valid); numBitsStat numBitsStat length(data_valid); end ber(idx) numErr / numBitsStat; end % 理论BPSK/MSK相干误码率 EbN0_lin 10.^(EbN0_dB / 10); ber_theory qfunc(sqrt(2 * EbN0_lin)); figure; semilogy(EbN0_dB, ber, o-); hold on; semilogy(EbN0_dB, ber_theory, x-); grid on; xlabel(Eb/N0 (dB)); ylabel(BER); legend(MSK 仿真, MSK 相干理论, Location, southwest);MSK相干解调的误码率理论值和BPSK完全一样都是Pb Q(sqrt(2Eb/N0))。原因就是MSK可以拆成两个正交的BPSK支路每支路数据率是总比特率的一半合并之后总误码率还是BPSK的理论公式。第一次看到这个结论可能会觉得意外但仿真结果会证明它是对的。正常情况下仿真曲线和理论曲线的差距应该在0.3dB以内。如果差距明显偏大先检查是不是噪声方差算错再看匹配滤波窗口是否正确。很多人在这一步会不小心把Eb算成“符号能量”而不是“比特能量”导致整条曲线右移3dB这是高频问题。5. 常见问题与排查实录5.1 常见问题速查表下面这张表是我在帮别人排查MSK仿真问题时最常用的一份对照表基本覆盖了从调制到解调最容易踩的坑。问题现象可能原因排查方向相位路径出现跳变Q路延迟没有做或延迟为2Tb检查Q_wave的索引平移长度是否等于sps误码率稳定在0.5附近判决前没有做半正弦相关用匹配滤波代替码元中心直接抽样误码率BER曲线整体右移3dBEb和N0换算时把符号能量当成了比特能量确认Eb mean(abs(base_tx).^2) * Tb通带版本BER比等效基带差很多低通滤波器群延迟未补偿滤波后统一扣除delay个采样点信号包络不是恒定包络I/Q延迟或加权函数不匹配分别画I/Q基带波形并与本地参考对比仿真速度慢每符号都用for循环numBits不超过几千即可课设足够最后一个符号判决报错Q路窗口超出信号尾部丢弃最后一个Q符号或帧尾补尾符号这张表里最常出问题的还是前三条。尤其“误码率0.5”十个报这个问题的人里八个都是因为用了中心点直接抽样没有做半正弦相关。我在一次帮别人调试时对方拿着星座图看了半天始终不明白为什么抽出来的点全在零点附近直到我把MSK的I/Q支路波形画出来他才意识到中心点对应的是半正弦过零点。5.2 几个只有实战才会踩的细节第一个细节是滤波器的边界处理。如果你选择做完整通带下变频fir1低通滤波器会产生一个固定的群延迟也就是filter输出信号相对于输入信号整体延后了(length(b)-1)/2个采样点。这个延迟不补偿后面的抽样点和相关窗口全部会错位。最简单的做法是滤波后直接从delay1开始取数据让后续所有处理都建立在“对齐后”的信号上。第二个细节是随机种子。蒙特卡洛仿真里每次重新生成随机数据如果不固定随机种子不同EbN0点之间会有随机波动导致BER曲线毛刺很多。建议在仿真外层用rng设定全局种子比如rng(2025)这样结果可复现。如果追求更平滑的曲线可以在每个信噪比点增加蒙特卡洛次数而不是增大单次数据长度因为短帧更容易出现统计波动。第三个细节是时域波形和相位路径只适合用小数据量展示。我建议展示波形和相位路径时用numBits400以下画出来清晰做BER仿真时再用更长数据比如每帧400bit乘20次。不要用一个超长序列同时干两件事否则波形图细节会被压缩到看不清。最后再说一个我自己觉得特别值得养成的习惯写完调制端代码后第一件事不是急着写解调而是先画相位路径。如果相位路径不是一条连续折线后面做再多解调都是白搭。我在第一次做MSK时就是在这里栽了跟头——Q路延迟少了一个Tb相位路径每隔两个符号就出现一次跳变当时还以为公式写错了后来对着教材逐行核对才发现是索引平移出了问题。这个习惯救了我很多次也希望你能用上。