咱们今天不整那些虚的直接拿MATLAB把2FSK调制解调从头到尾捋一遍。标题里写了重点对比包络检波和相干解调那这篇就把两套解调方案都做出来从原理到代码再到误码率曲线一步都不落下。不管你是通信课设要用还是面试前临时抱佛脚或者单纯想搞懂这两个解调方式到底差在哪这篇文章都能让你照着重现一遍。先交代一下代码环境MATLAB R2021a及以上版本都能跑不需要额外工具箱所有代码纯手写复制粘贴就能出结果。1. 2FSK为什么值得手动仿真一遍2FSK全称是二进制频移键控用两个不同频率的载波代表二进制里的0和1。它和2ASK最大的区别在于信息承载在频率上而不是幅度上所以抗幅度衰落的能力更强这也是它在无线对讲、低速数传里依旧占有一席之地的原因。但深入做仿真之前有几个概念必须掰扯清楚不然代码写出来也只是一堆数组拼接出了问题根本不知道错在哪。1.1 FSK信号的本质两个频率的开关切换2FSK的时域表达式可以写成s(t) A * cos(2π * f1 * t φ1) 表示发送“0” s(t) A * cos(2π * f2 * t φ2) 表示发送“1”这里最关键的不是两个频率本身而是切换是否连续。相位连续FSKCPFSK的频率切换过程平滑频谱旁瓣衰减得快相位不连续FSK在码元切换瞬间会有相位跳变频谱扩展明显实际工程中一般会做额外的滤波处理。用MATLAB仿真时有两种建模思路一种是根据码元状态分别生成两段正弦波然后拼接另一种是先构造频率控制序列再用频率调制函数。我推荐新手用第一种因为它把FSK的物理含义直接体现在代码里每个码元就是一个固定频率的余弦片段拼接起来就是完整的FSK波形。1.2 仿真参数怎么定才能说明问题参数设计直接决定仿真结果的可靠程度。对比包络检波和相干解调核心指标是误码率随着信噪比的变化所以参数设计必须保证采样率足够高、码元周期内有足够的采样点、两个载波频率不互相干扰。我用的参数如下码元速率Rb 1000 bps即每个码元持续1ms载波频率f1 3000 Hz代表“0”f2 5000 Hz代表“1”采样率fs 50000 Hz每个码元采样点数 fs / Rb 50个点为什么f1取3000Hz、f2取5000Hz这两个频率间距2000Hz正好是码元速率的两倍保证了两个频率的频谱不会重叠太多。采样率取50kHz意味着最高信号频率5kHz还有十倍的过采样余量还原正弦波形完全够用。实际上工程里有更讲究的频率间隔选取方式最小频移键控MSK把频偏压缩到码元速率的0.5倍但那是另一个话题咱们这次就不展开了。2. 调制端代码从比特序列到FSK波形模拟通信系统的第一步是生成发送端的信号。这里看起来简单但“生成比特序列”这个环节就有一个隐蔽的坑不少新手在这里翻车。2.1 随机比特序列生成与两个小坑用randi生成随机序列我得提醒两件事。第一randi([0 1], 1, N)返回的是double类型的0/1如果直接拿去当索引用会出问题第二第一次跑仿真之前务必固定随机种子否则每次都得到不同曲线无法判断代码改动的影响。固定种子的方式很简单rng(42); % 固定随机种子保证仿真可复现 data randi([0 1], 1, 1000); % 生成1000个随机比特这里固定随机种子不是偷懒而是工程仿真的基本素养。调试阶段不锁随机源基本等于闭着眼睛改代码。2.2 码元映射与载波拼接实现接下来把0映射到f11映射到f2为每个码元生成对应的正弦片段。这段代码是调制端的核心也是整个系统的地基fs 50000; % 采样率50kHz Rb 1000; % 码元速率1000bps f1 3000; % 代表0的载波频率3kHz f2 5000; % 代表1的载波频率5kHz samplesPerBit fs / Rb; % 每个码元的采样点数 N length(data); t_bit (0:samplesPerBit-1) / fs; % 单个码元的时间轴 % 预分配调制信号数组 modSignal zeros(1, N * samplesPerBit); % 逐码元生成FSK信号 for k 1:N idx (k-1)*samplesPerBit 1 : k*samplesPerBit; if data(k) 0 modSignal(idx) cos(2*pi*f1*t_bit); else modSignal(idx) cos(2*pi*f2*t_bit); end end这里有个MATLAB性能小技巧预分配modSignal数组避免循环中动态增长。1000个码元还算小事但如果仿真几百万比特动态增长的数组会让运行时间从秒级变成分钟级这不是危言耸听。2.3 观察一下时域波形和功率谱确认调制正确信号生成之后先别急着加噪声看一眼波形再往下走。取前50个码元画时域波形figure; plot((0:length(modSignal)-1)/fs * 1000, modSignal); xlabel(时间 (ms)); ylabel(幅度); title(2FSK调制信号时域波形); xlim([0 5]); % 只看前5ms如果调制正确你会看到波形在3kHz和5kHz之间切换密集程度有明显差异。再用pwelch看功率谱两个频率处会出现明显的谱峰。这一步能快速确认f1和f2是否真的落在预期位置也能看出相位不连续导致的频谱扩展——旁瓣掉得不够快那就是相位不连续的直接体现。3. 高斯白噪声信道与信噪比换算调制信号不经过信道直接解调得到的一定是“完美结果”毫无参考价值。真正的通信系统必然经历噪声干扰所以必须在接收端加上高斯白噪声并且对于每个信噪比都要独立做一次蒙特卡洛仿真。3.1 AWGN信道的MATLAB实现给信号加噪声MATLAB里最直接的方式是awgn函数SNR_dB 10; % 信噪比10dB rxSignal awgn(modSignal, SNR_dB, measured);第三个参数写成measured表示函数会根据输入信号的实际功率自动计算噪声方差这样能保证信噪比的准确性比自己手动算噪声功率靠谱得多。这一段最简单的代码背后藏着一个容易忽略的问题——awgn默认假设输入信号是实信号如果输入复数信号会产生错误的噪声功率计算这次仿真是实信号所以没问题但如果你以后仿真QPSK一定要记得改用复数噪声添加方式。3.2 不同信噪比下波形长什么样拿SNR0dB和SNR15dB两种情况做个对比0dB时噪声幅度和信号幅度几乎一样肉眼已经很难从时域波形直接分辨频率切换15dB时波形仍然干净频率切换一目了然。这就是后面误码率曲线会呈现出来的基本趋势低信噪比下解调恢复正确信息的难度陡增而在高信噪比下两种解调方式的差距在缩窄。理解这一层再去看误码率曲线就会有更直观的感知。4. 包络检波解调实现简单但别轻视理论包络检波的本质是把频率差异转化成功率差异。2FSK信号通过两个中心频率分别为f1和f2的窄带带通滤波器再各自做包络提取比较两个包络大小就可以判断当前码元是0还是1。整个过程中频率信息被“丢弃”了提取的是幅度包络所以这个方法也叫非相干解调。4.1 带通滤波器设计与代码首先要为两条支路分别设计带通滤波器。MATLAB里最方便的是designfilt函数fpass1 [2500 3500]; % 支路1通带中心频率3000Hz fpass2 [4500 5500]; % 支路2通带中心频率5000Hz bpFilt1 designfilt(bandpassiir, ... FilterOrder, 8, ... HalfPowerFrequency1, fpass1(1), ... HalfPowerFrequency2, fpass1(2), ... SampleRate, fs); bpFilt2 designfilt(bandpassiir, ... FilterOrder, 8, ... HalfPowerFrequency1, fpass2(1), ... HalfPowerFrequency2, fpass2(2), ... SampleRate, fs);滤波器设计里的取舍值得多说两句。滤波器阶数越高通带越平坦、过渡带越窄但相位延迟也会变大导致波形失真。8阶对于这个仿真场景足够再高就会出现滤波后码元之间相互拖尾干扰的情况。滤波器的通带宽度也不能太窄否则码元切换瞬间的高频分量被削掉包络会变得圆润边缘模糊判决点不好选。4.2 包络提取与判决包络提取最简单的办法是对滤波后的信号求绝对值再经过一个低通滤波器或者移动平均得到平滑的包络曲线。这里用movmean移动平均来做env1 abs(filter(bpFilt1, rxSignal)); env2 abs(filter(bpFilt2, rxSignal)); windowSize 10; % 滑动平均窗口大小 env1 movmean(env1, windowSize); env2 movmean(env2, windowSize);包络提取后在每个码元的判决时刻比较两支路的包络大小。抽样点取在码元周期的中点附近最合理因为码元切换瞬间的暂态已经结束包络相对稳定% 抽样判决 samplingPoints round(samplesPerBit/2 : samplesPerBit : length(env1)); decoded_data_coherent double(env1(samplingPoints) env2(samplingPoints));这里用而不是是防止包络相等时误判实际信号中两个包络相等的概率极低但代码逻辑上要先定义清楚。4.3 包络检波的误码率统计误码率的计算直接统计判决结果和原始数据的差异errors sum(decoded_data_coherent ~ data); ber_coherent errors / N;到这里包络检波就完整跑通了。但别高兴太早跑完第一次仿真之后大概率会发现一个匪夷所思的现象——不论SNR多高、误码率都降不下去。原因往往出在两个地方一个是滤波器相位延迟导致抽样点偏移另一个是移动平均窗口太大导致包络变化被平滑过头。我建议你在包络检波仿真时画一条env1和env2的曲线和原始信号叠加在一起看亲眼确认抽样时刻包络是不是能正确区分0和1。5. 相干解调两路乘法器加低通滤波器相干解调走的是另一条路线接收端需要产生与发送端同频同相的本地载波将接收信号分别与两路本地载波相乘再做低通滤波恢复基带信号最后比较两路输出。MATLAB仿真时我们可以简化一个操作——直接用发送端的载波作为本地载波这就等于给了接收端“上帝视角”但现实系统中载波同步是要专门解决的难题。5.1 本地载波相乘并用低通滤波器提取基带信号相干解调的MATLAB实现如下% 生成本地载波 t_total (0:length(rxSignal)-1) / fs; localCarrier1 cos(2*pi*f1*t_total); localCarrier2 cos(2*pi*f2*t_total); % 乘法器输出 mixSignal1 rxSignal .* localCarrier1; mixSignal2 rxSignal .* localCarrier2; % 低通滤波器设计 lpFilt designfilt(lowpassiir, ... FilterOrder, 8, ... HalfPowerFrequency, 1500, ... SampleRate, fs); basebandSignal1 filter(lpFilt, mixSignal1); basebandSignal2 filter(lpFilt, mixSignal2);这个低通滤波器的截止频率选1500Hz道理在于乘法器输出中包含直流分量信息所在和位于2f1、2f2附近的二倍频分量。低通滤波器需要把这些二倍频分量滤除掉只保留基带分量截止频率过高就滤不干净过低又会把码元波形压扁。1500Hz对1kbps的码元速率是合理的折中。注意这里用的是filter而不是conv因为我们要保持输出数组长度与输入一致后续抽样才不需要额外处理索引。5.2 抽样判决的实现细节低通滤波之后两路输出分别反映了“这个码元与f1的相关程度”和“这个码元与f2的相关程度”。对每个码元周期取中点抽样比较两路大小samplingPoints round(samplesPerBit/2 : samplesPerBit : length(basebandSignal1)); decoded_data_coherent double(basebandSignal1(samplingPoints) basebandSignal2(samplingPoints));这里有一个MATLAB新手经常踩的坑不要提前对basebandSignal做归一化或标准化除非你有明确理由。相干解调两路信号在同一信道条件下受到同一个噪声影响它们之间的相对大小才是判决依据任何独立的幅度归一化操作都会破坏这种相对比较关系。5.3 为什么相干解调理论误码率更低相干解调的误码率理论上优于包络检波约1.5dB原因在于相干的乘法器实际上是一个相关器把信号能量集中到了基带直流分量上而噪声经过低通滤波后统计特性是均匀分布在宽带内的两者相乘后信噪比获得了处理增益。包络检波本质上是能量检测没有利用信号的相位信息在低信噪比下更容易被噪声把包络顶起来导致误判。6. 蒙特卡洛仿真与误码率曲线对比单次仿真只能得到一个点的误码率要画出完整的曲线必须在多个信噪比下重复仿真。这里用蒙特卡洛法对每个SNR值做若干次独立仿真求平均误码率。6.1 完整蒙特卡洛仿真循环这段代码会遍历从-2dB到15dB的信噪比每一步运行50次独立仿真SNR_dB_list -2:1:15; numTrials 50; numBits 2000; ber_env zeros(size(SNR_dB_list)); ber_coh zeros(size(SNR_dB_list)); for snrIdx 1:length(SNR_dB_list) SNR_dB SNR_dB_list(snrIdx); errEnvTotal 0; errCohTotal 0; for trial 1:numTrials rng(trial * snrIdx); % 每次试验独立随机种子 data randi([0 1], 1, numBits); modSignal zeros(1, numBits * samplesPerBit); for k 1:numBits idx (k-1)*samplesPerBit 1 : k*samplesPerBit; if data(k) 0 modSignal(idx) cos(2*pi*f1*t_bit); else modSignal(idx) cos(2*pi*f2*t_bit); end end rxSignal awgn(modSignal, SNR_dB, measured); % 包络检波 env1 abs(filter(bpFilt1, rxSignal)); env2 abs(filter(bpFilt2, rxSignal)); env1 movmean(env1, 10); env2 movmean(env2, 10); samplingPoints round(samplesPerBit/2 : samplesPerBit : length(env1)); decodedEnv double(env1(samplingPoints) env2(samplingPoints)); errEnvTotal errEnvTotal sum(decodedEnv ~ data); % 相干解调 mix1 rxSignal .* localCarrier1; mix2 rxSignal .* localCarrier2; bb1 filter(lpFilt, mix1); bb2 filter(lpFilt, mix2); samplingPoints round(samplesPerBit/2 : samplesPerBit : length(bb1)); decodedCoh double(bb1(samplingPoints) bb2(samplingPoints)); errCohTotal errCohTotal sum(decodedCoh ~ data); end ber_env(snrIdx) errEnvTotal / (numBits * numTrials); ber_coh(snrIdx) errCohTotal / (numBits * numTrials); end每个SNR点做50次独立仿真2000个比特总统计量为50*2000100000个比特。这个量级算误码率在10^-3左右还能保证统计稳定性再低的误码率就需要更多的仿真次数否则曲线会剧烈抖动。6.2 理论误码率曲线对比把仿真结果和理论公式放到同一张图上。2FSK相干解调的理论误码率为Pe_coherent 0.5 * erfc(sqrt(Eb/N0/2))包络检波的理论误码率为Pe_envelope 0.5 * exp(-Eb/(2*N0))这里要注意理论公式中的Eb/N0和信噪比SNR之间的换算Eb_N0_dB SNR_dB_list - 10*log10(samplesPerBit/2);为什么会有这个换算因为awgn函数的SNR定义是信号功率与噪声功率的比值而误码率公式中的Eb/N0是每比特能量与噪声功率谱密度之比。对于2FSK信号在采样率fs下每个码元有samplesPerBit个采样点信号能量分布在频率维度上换算关系不同导致两者差了一个10log10(samplesPerBit/2)的偏移量。具体到这里samplesPerBit50所以偏移量是10log10(25)≈14dB。不搞清楚这个换算你会发现仿真曲线比理论曲线整体平移了十几dB误以为自己代码写错了。6.3 两种解调方式在低信噪比下的真实差距画完图你会看到非常典型的结果高信噪比下10dB以上两条曲线都几乎贴着零轴已经看不出明显差别但在0dB附近相干解调的误码率明显低于包络检波差距大约在1~2dB。这个结果和理论预期吻合因为相干解调利用相位信息获得了更多的有效信号能量而包络检波靠幅度信息天生处在劣势。蒙特卡洛仿真次数如果太少低误码率区域会出现曲线不平滑甚至某些点误码为0导致对数坐标画不出来。遇到这个问题不要紧张要么增加仿真次数要么减小低误码率区域的仿真点数范围。7. 代码搬运过程中的常见报错和解决方式写这篇博客之前我又把完整的仿真过程从头到尾跑了一遍过程中确实有几个小地方很容易踩坑这里集中列出来。7.1 滤波器输出长度与抽样索引不一致filter函数的输出长度默认和输入一致这点没问题。但如果有人用了conv或者conv2输出长度就变成了L_in L_filt - 1导致后续抽样索引超出数组边界。解决方案是使用filter而不是conv如果非要使用卷积记得截断输出到输入长度。7.2 随机种子问题导致误码率曲线抖动剧烈我之前调试时发现每个SNR点只用一次随机序列画出来的误码率曲线锯齿严重根本没法看。后来改成多次试验取平均曲线才平滑下来。这个方法在通信仿真里叫蒙特卡洛平均是处理随机数据必不可少的步骤。7.3 抽样时刻偏移导致误码率虚高滤波器是有群延迟的尤其IIR滤波器在中高频段会有明显的相位非线性。designfilt设计的滤波器虽然具有零相位响应使用filtfilt或者线性相位使用FIR滤波器但普通IIR滤波器的filter输出在起始阶段会有暂态过渡导致前几个码元的抽样结果异常。解决方式有两种仿真时丢弃前几个码元比如50个或者把抽样点适当往后移。实际工程里都会预留一段前导码或训练序列专供滤波器过渡仿真时也要有这个习惯。7.4 理论曲线画不出来或者错位最常见的原因是erfc函数里忘记除以2或者Eb/N0换算错误。我的建议是先把仿真曲线和理论曲线画在同一个坐标下找到偏移量后反向验证自己的换算是否正确。如果你发现曲线形状完全一致但水平方向差了一个固定值那基本就是换算关系没对上。8. 扩展思考代码还能往哪个方向迭代到这里2FSK的调制解调完整仿真已经跑通了。但收尾之前我想多说几句怎样在这个基础上继续往深处走据我经验很多刚接触通信仿真的人都栽在“只会跑通不会改”这个阶段。8.1 把相位不连续改成连续相位FSK上面的代码是直接把两个频率的独立余弦波拼在一起切换瞬间相位会有跳变。你可以试验在码元切换时调整下一段波的起点相位使其与上一段波的终点相位连续构造一个连续相位FSK波形对比两者的功率谱特性。实际操作起来非常有意思连续相位FSK的频谱旁瓣下降更快但实现上需要多一行相位累积的代码。8.2 加入频偏估计补偿仿真中的本地载波是完美的现实中接收端与发送端存在频率偏差可能是几赫兹到几十赫兹。你可以人为给接收信号加一个频率偏移再实现一个简单的频偏估计和校正算法看看误码率会恶化到什么程度。这个改进会让你对“载波同步为什么难”有切身体会。8.3 从2FSK到MFSK的推广2FSK只是最基础的二进制频率调制可以扩展为4FSK、8FSK甚至是16FSK每个码元携带更多比特信息。MFSK在低信噪比下表现优异但代价是占用更大的带宽仿真中对比不同M值的误码率曲线会很有趣但需要关注的参数和滤波器设计复杂程度也会上一个台阶。8.4 跟理论误码率公式印证的一个小技巧最后分享我在调这个仿真时用的小技巧先把仿真曲线画出来再画理论曲线如果两者对不上不要在第一时间怀疑公式推导而是检查你的Eb/N0换算是否准确。绝大多数情况下仿真代码只要逻辑没写错曲线应该和理论值非常接近。真正理解了这一点2FSK这块基本就吃透了。
