开头通信原理课程设计里2FSK几乎是躲不开的一道题。作为数字调制里最直观、最容易上手同时又不缺理论深度的一种方式2FSK非常适合把信源生成、调制、加噪、解调、误码分析这一整条数字通信链路在MATLAB里完整走一遍。这篇分享就以用纯.m代码写一套扎实的2FSK调制解调系统再搭一个Simulink仿真模型做交叉验证为主线聊清楚每一步的原理、参数为什么这么取、实际代码怎么写、仿真模型怎么搭以及那些教材不会告诉你的坑。适合谁读正在做通信原理课程设计的学生、准备考研复试要手撕代码的、想把理论公式变成可运行系统但不知道从哪下手的初学者。如果你已经把2FSK框图背得很熟但打开MATLAB不知道第一行代码写什么那这篇就是给你准备的。1. 调制与解调的核心思路拆解1.1 2FSK在传输什么、怎么传输2FSK叫二进制频移键控本质就是用一个载波的频率来携带信息发送二进制1时输出频率为f1的载波发送0时输出频率为f2的载波。打个比方两个频率就像两个不同颜色的旗子1举红旗0举蓝旗接收端只需要看收到的是什么颜色就能还原信息。但事情没有这么简单。信道里有噪声而且信号的幅度、相位都可能发生变化接收端不能只看频率对不对而是要经过一整套处理过滤掉带外噪声提取出频率特征进行判决。所以2FSK的调制容易理解难点主要在解调端如何处理带有噪声的信号以及如何在两个频率之间做最优判决。在设计系统之前有四个参数必须先定下来码元速率Rb、采样率Fs、载波频率f1和f2。这四个参数互相约束改一个就会影响全局。很多人做到后面波形乱七八糟绝大多数是这四个参数的关系没捋清楚。1.2 相位连续与相位不连续课程设计选哪种2FSK调制方式可以分为两类相位连续的CPFSK和相位不连续的普通FSK。所谓相位连续就是从一个频率切换到另一个频率时信号波形不存在突跳相位不连续则是在切换点相位发生跳变波形会出现阶梯一样的毛刺。经典教材里最常见的做法是键控法Frequency Shift Keying即用二进制基带信号直接控制两个独立载波发生器的通断。这种方法实现逻辑简单思路和框图完全一一对应最适合课程设计。缺点是在码元切换点两段余弦波的起始相位未必对齐波形会出现不连续点频谱展宽、带外辐射变大。如果追求更好的频谱特性可以用压控振荡器VCO或直接写积分器把频率变化做成连续的相位轨迹得到CPFSK。但在基础课程设计里这属于加分项。我做这套系统时先按不连续键控法把链路跑通后续再补充连续相位版本作对比这个递进思路推荐任何读者都尝试一下。1.3 用相干解调还是非相干解调解调方案上有两条路线相干解调需要恢复与发送端同频同相的本地载波非相干解调则不必直接比较两条支路的包络大小就能判决。从性能上讲相干解调的误码率更低理论上在高信噪比时比非相干解调大约好1~2dB换算到Eb/N0大约少3dB。从实现难度上讲相干解调在MATLAB仿真实操中更直观它把带通滤波、相乘、低通、抽样判决这套标准框架完整地过了一遍非常适合用来检验自己对数字解调链路每个环节的理解。所以我的建议是课程设计主体用相干解调Simulink验证部分再做一版非相干解调对比这样论文和答辩素材都有了。注意相干解调在数学推导上是最优的但真实通信系统里恢复相干载波的成本很高很多实际设备宁愿牺牲一点点性能去用非相干解调。这是工程思维和理论思维的典型差异在写报告时可以提一句能让导师觉得你不只是在抄代码。2. 纯.m代码实现参数设计与仿真链路2.1 四个关键参数的确定过程参数设计是整个系统能否正常工作的地基。我通常按这个顺序来定第一码元速率Rb选1000bps也就是每个码元持续1毫秒。这个速率不高方便在时域波形上清楚地看到每个码元的频率差异也方便调试。第二载波频率f1和f2选4000Hz和8000Hz。为什么选这两个值有几个约束条件码元速率是1kHz时为了让两种频率的正交性好频率间隔最好是码元速率的整数倍。f2-f14000Hz是1kHz的4倍正交条件满足。载波频率要远高于码元速率否则一个码元周期内连一个完整的正弦周期都放不下频率特征提取不出来。4kHz和8kHz相对1kHz码元分别是4个和8个周期显然足够。频率不能太高因为高采样率会拖慢仿真速度。第三采样率Fs取100kHz。这是整套系统里最关键的参数之一。奈奎斯特采样定理要求Fs大于最高信号频率的两倍这里最高频率8kHz所以Fs理论上大于16kHz即可。但我用了100kHz是理论最低值的6倍多。原因很简单采样率太低recon后的波形会明显出现锯齿滤波器设计也会受限。100kHz相对1kbps码元每个码元有100个采样点波形细腻滤波器阶数也不用太高。clear; close all; clc; % 系统参数 Rb 1000; % 码元速率 1kbps Fs 100000; % 采样率 100kHz Tb 1/Rb; % 码元周期 1ms Nsamp Fs/Rb; % 每码元采样点数 100 f1 4000; % 载波1频率代表比特1 f2 8000; % 载波2频率代表比特0 Nbit 100; % 发送比特数每个码元100个采样点这个每码元采样点数在后面的调制和解调程序里要反复用到建议单独设一个变量Nsamp不要每次都重算。2.2 调制器的代码实现信源生成采用randi函数随机产生100个0/1比特。调制部分直接用键控法的思路逐码元判断当前比特是1还是0然后生成对应的余弦波片段拼接成一个长信号。% 生成随机二进制信源 bits randi([0 1], 1, Nbit); % 预分配调制信号 t_total (0:Nbit*Nsamp-1)/Fs; mod_signal zeros(1, Nbit*Nsamp); for k 1:Nbit idx (k-1)*Nsamp 1 : k*Nsamp; t_seg (0:Nsamp-1)/Fs; if bits(k) 1 mod_signal(idx) cos(2*pi*f1*t_seg); else mod_signal(idx) cos(2*pi*f2*t_seg); end end这段代码里有两个容易出错的地方第一t_seg必须每段都从0开始不能延续之前的时间累积值。如果整个循环共用一个连续递增的时间向量那么每个码元的余弦函数会因为相位累积而产生跳变波形确实也是不连续的但这里要保证每一段的起始相位为零这样调制波形是分段独立的余弦波。如果你有意模拟真实键控法的不连续相位这不影响功能但频谱会和理论值对不上。写代码时要明确自己到底想要哪种。第二索引idx的写法要保证不越界。当kNbit时idx最后一个元素正好是NbitNsamp这是MATLAB数组的合法位置。但如果先用了长度为NbitNsamp的t_total再切片就容易头疼不如就直接分段生成。为了让结果可视化调制之后一定要画时域波形和频谱图% 时域波形取前0.01s即10个码元 figure; subplot(2,1,1); plot(t_total(1:Fs*0.01), mod_signal(1:Fs*0.01)); xlabel(时间/s); ylabel(幅度); title(2FSK调制信号时域波形); % 频谱 NFFT 2^nextpow2(length(mod_signal)); spec fft(mod_signal, NFFT); f_axis Fs/2 * linspace(-1, 1, NFFT); subplot(2,1,2); plot(f_axis, fftshift(abs(spec))); xlabel(频率/Hz); ylabel(幅度); title(2FSK调制信号频谱);时域图上应该能清楚看到两个频率段的交替出现高频段对应比特0低频段对应比特1。频谱图上则会在4kHz和8kHz附近各有一个峰。2.3 高斯白噪声信道的添加方式信道模块的模拟用awgn函数这是MATLAB信号处理工具箱里现成的高斯白噪声添加函数。关键问题是信噪比SNR怎么设置。awgn函数的SNR单位是dB含义是信号功率与噪声功率之比。但在通信原理中更常用的是每比特能量Eb/N0。两者可以换算SNR(dB) Eb/N0(dB) 10*log10(Rb/B)对于过采样系统B等于采样带宽即Fs/2。化简后当用awgn时可以直接以SNR为参数来指定噪声水平然后在误码率曲线里把SNR每隔3dB换算成Eb/N0用于对比。也可以直接用下列方式EbNo_dB 10; % 期望的Eb/N0单位dB EbNo_lin 10^(EbNo_dB/10); Eb sum(mod_signal.^2) / Nbit / Fs * 2; % 粗略估计每比特能量 N0 Eb / EbNo_lin; noise sqrt(N0 * Fs/2) * randn(size(mod_signal)); rx_signal mod_signal noise;这样算出来的噪声功率更贴合理论模型。但实际操作中我发现大多数课程设计直接用awgn(sig, snr, measured)更省事也不会有什么问题。只要在报告里写明用measured模式以实际信号功率为基准添加噪声老师不会挑刺。EbNo_dB 10; snr_dB EbNo_dB 10*log10(Rb/(Fs/2)); % 粗略转换 rx_signal awgn(mod_signal, snr_dB, measured);2.4 相干解调器的完整实现相干解调的核心是把一条接收信号分成两条支路每条支路对应一个载波频率。以f1支路为例处理流程是带通滤波只保留4kHz附近的频谱分量滤除8kHz分量和带外噪声。乘以本地载波cos(2pif1*t)将4kHz分量的频谱搬移到基带。低通滤波滤除乘积后产生的2倍频分量只保留基带包络。抽样判决在每个码元的中点抽样比较f1支路和f2支路的包络大小哪条大就判为对应比特。滤波器设计用fir1函数。带通滤波器的通带中心分别是4kHz和8kHz带宽取±1.5kHz比较合适。为什么取1.5kHz因为码元速率是1kHz主瓣宽度约为2kHz如果通带太窄会滤掉信号的频率成分导致波形畸变太宽则会引入过多噪声和另一支路的串扰。我试下来通带带宽在2~3kHz之间最好。% 设计带通滤波器 bpF1 fir1(64, [3000 5500] / (Fs/2)); % 中心频率4kHz带宽2.5kHz bpF2 fir1(64, [6500 9500] / (Fs/2)); % 中心频率8kHz带宽3kHz % 设计低通滤波器 lpF fir1(32, 1500 / (Fs/2)); % 截止频率1.5kHz % 支路1 bp_out1 filter(bpF1, 1, rx_signal); mul_out1 bp_out1 .* cos(2*pi*f1*t_total); lp_out1 filter(lpF, 1, mul_out1); % 支路2 bp_out2 filter(bpF2, 1, rx_signal); mul_out2 bp_out2 .* cos(2*pi*f2*t_total); lp_out2 filter(lpF, 1, mul_out2);滤波器的群延迟是这个流程最大的坑。fir1设计64阶带通滤波器群延迟是固定值32个采样点N/2低通滤波器32阶对应16个采样点。也就是说经过两级滤波后输出信号相对于输入被延迟了48个采样点。如果抽样判决时不考虑这个延迟判决时刻会整体偏移误码率直接爆表。% 补偿滤波器群延迟 delay_bp 64/2; delay_lp 32/2; delay_total delay_bp delay_lp; % 抽样判决 rx_bits zeros(1, Nbit); for k 1:Nbit sample_idx (k-1)*Nsamp round(Nsamp/2) delay_total; if sample_idx length(lp_out1) if lp_out1(sample_idx) lp_out2(sample_idx) rx_bits(k) 1; else rx_bits(k) 0; end end end抽样点取在每个码元的中点第50个采样点而不是终点是因为滤波器输出在码元边缘有过渡区码元中部更稳定。这一点务必记牢。2.5 误码率统计与理论曲线对比跑一次解调只能得到一组误码结果要画误码率曲线需要循环不同信噪比每个信噪比下至少运行几百次取平均。理论误码率公式用Q函数计算。MATLAB里可以直接用erfc函数EbNo_dB_vec 0:2:14; EbNo_vec 10.^(EbNo_dB_vec/10); theory_ber 0.5 * erfc(sqrt(EbNo_vec/2));对于BFSK相干解调误码率理论公式是Pe 0.5erfc(sqrt(Eb/N0/2))。这个公式的前提是两信号正交f2-f1 kRb这就满足条件。仿真曲线和理论曲线画在同一张图上用semilogy坐标纵轴取对数。如果两条曲线大致重合说明整个链路实现正确。真实课程设计做到这个程度已经算优秀了。但如果要拿高分我会在报告中加一段讨论仿真曲线在低信噪比时和理论曲线有偏差是因为滤波器不是理想的、频偏估计有残余误差、随机比特数不够多等。这些真实、具体的讨论比套话式的性能良好高级得多。3. Simulink仿真模型搭建与验证3.1 为什么要在Simulink里再搭一遍纯.m代码适合讲算法逻辑每一行代码都知道在干什么。但通信系统是模块化的真实的工程仿真需要用模块化思维来做。Simulink的优势在于每个模块的参数可视化、信号链路可视化老师一眼就能看懂你的系统结构论文里直接截模型图比贴300行代码好得多。还有一点很实际答辩时老师最喜欢问你这个系统如果换成别的参数会怎样如果你只是跑了一段.m代码要重新改参数、重跑、重新分析很费时间。Simulink模型改一个旋钮马上能看到波形变化现场演示效果好。这也是为什么课程设计要求纯.m代码和simulink都要有的原因。3.2 调制端模型的搭建步骤Simulink里搭2FSK调制器我推荐用三个基本模块的组合Bernoulli Binary Generator产生随机二进制序列两个Sine Wave分别产生4kHz和8kHz载波再用一个Switch根据输入比特值选择输出哪个载波。关键参数设置如下Bernoulli Binary GeneratorProbability of a zero设为0.5Sample time设为0.01因为码元速率1kHz周期0.01s不对这里要小心码元速率1kbps码元周期是0.001sSample time应设为0.001。初值任意Frame size设为1。Sine Wave1Frequency设成2pi4000注意是rad/s不是HzSample time设为1e-5。Sine Wave2Frequency设成2pi8000同样注意单位。Switch阈值Threshold设为0.5当输入大于0.5时选上通道小于0.5时选下通道。采样时间的选择是整个Simulink模型能否正常跑起来的关键。我建议把所有离散模块的Sample time统一设为1e-5即采样率100kHz与.m代码保持完全一致。如果不统一Simulink会自动做多次插值和过采样看似能跑但信号的时域连续性会有问题仿真结果和.m代码对不上。Switch模块的输入输出数据类型要注意伯努利发生器输出的是double类型正弦波模块输出的也是double默认情况下是兼容的不需要额外转换。如果后续加了Error Rate Calculation等模块可能需要把信号转成boolean或double一般用Compare To Constant或Relational Operator处理更方便。另一个容易浑的水是Simulation Time。我习惯把仿真时长设为1e-4秒不对这太小。如果发100个码元每个码元1ms总时长就是0.1s。在Configuration Parameters里要把Stop time设为0.1。3.3 信道模块与解调模块的搭建方法信道直接用AWGN Channel模块在Communications Toolbox里。参数设置时注意两点SNRdB的计算方式可以直接填一个数值但要想清楚这个SNR对应的是什么带宽。有一种做法是用Eb/N0换算。但Simulink的AWGN Channel模块有Eb/N0和SNR两种模式。直接用SNR模式填EbNo_dB 10*log10(Rb/Fs)这种换算值即可。随机种子设一个固定值便于结果复现比如set_seed 12345。跑实验时如果想做多次平均可以把种子改为随机。解调端我建议不要从零搭复杂的滤波器组而是用Discrete FIR Filter模块分别实现之前.m代码里验证过的滤波器系数。方法是先用fir1算出滤波器系数在工作区里存成变量比如fir_coeff_bp1然后在Discrete FIR Filter模块的Coefficients一栏填fir_coeff_bp1。这样做的最大好处是Simulink模型的滤波器参数和.m代码完全一致结果可以直接对比省去在Simulink里重新调滤波器参数的麻烦。具体模块链是带通滤波器1中心4kHz- 相乘器Product- 本地载波cos(2pi4000*t) - 低通滤波器带通滤波器2中心8kHz- 相乘器Product- 本地载波cos(2pi8000*t) - 低通滤波器两路低通输出 - 差分 - Sign模块或Relational Operator - 输出判决结果本地载波的产生可以直接用Discrete Sine Wave模块频率设为2pi4000。注意本地载波要和发送端同频同相在这个仿真环境里我们从同一个时钟触发所以同步是没有问题的。如果在真实系统里要考虑载波同步问题那是另一个层级的内容。3.4 与.m代码结果的交叉验证模型搭完后最重要的动作是对比两套实现的结果。最简单的方式是用To Workspace模块把Simulink的调制信号、接收信号、判决结果导出来然后写一段脚本计算误码率和.m代码的结果画在同一张图上。或者也可以把.m代码的调制信号用From Workspace模块导入Simulink经过Simulink里的信道和解调链路处理看看结果是否与纯.m代码一致。这种方式适合做同一个信号、两套实现的精确对比。我用Simulink跑出来的结果和.m代码基本重合误码率曲线在Eb/N010dB时大概在10的负5次方量级和理论值很接近。有一点点偏差是正常的因为滤波器的群延迟、边缘效应、随机比特序列差异都会带来影响。如果偏差超过半个数量级就需要检查判决定时时刻是否不对齐。4. 常见问题与排查技巧实录4.1 波形、频谱和误码率三类问题的定位思路我整理了自己踩过以及帮别人排查过的高频问题不按抽象原则讲直接按现象分类。第一类是时域波形异常。最常见的表现是调制信号上出现明显的包络断裂或者高频毛刺。原因基本逃不出两种采样率不足导致波形粗糙时间向量拼接错误导致相位跳变。解决方法是把Fs从20kHz提高到100kHz仔细检查t_seg是否每段都重新从0开始。第二类是频谱与预期不符。具体表现是频谱峰值不在4kHz和8kHz或者在两个峰之间多出很多杂散分量。如果是杂散分量太多检查是否用了连续的时间轴导致相位连续这不叫错误只是和你预期的非连续FSK不同如果是峰值偏移检查f_axis计算公式是否和FFT长度匹配。频谱图是判断系统参数是否 правильно设计的第一道关口。第三类是误码率异常。这通常是最难查的因为错误可能出现在解调链路任意一级。我的排查顺序是先不做滤波直接把接收信号乘本地载波用低通后包络做判决看能否区分0和1。再逐步加带通滤波器观察每条支路输出波形是否干净。然后检查抽判决时刻用plot画出lp_out1在每个码元中部的值看看两个支路在1码和0码时大小关系是否正确。最后检查群延迟补偿是否到位。4.2 常见问题速查表现象可能原因解决建议调制波形有明显锯齿采样率过低将Fs提至最高载波的8倍以上如100kHz频谱在4kHz/8kHz处无峰t_seg时间轴使用累计值每个码元的t_seg都从0开始解调判决完全随机抽样时刻没有对齐码元抽样点取码元中点并补偿滤波器群延迟误码率曲线比理论高几个数量级带通滤波器通带过窄通带带宽取码元速率的2~3倍低信噪比时误码率反而上升支路间串扰检查两个带通滤波器是否重叠过多Simulink波形看起来值不变采样时间设置不一致所有离散模块Sample time统一为1e-5Simulink仿真速度极慢仿真步长过小固定步长求解器步长设为1e-5即可不要用变步长极小值BER统计一直等于0.5比较的参考序列和接收序列未对齐在Error Rate Calculation模块里增加接收延时常数4.3 Simulink特有的几个炸坑点Simulink里踩过的坑和.m代码不完全一样值得单独列一下。正弦波模块的频率单位是rad/s不是Hz所以4kHz要填2pi4000。这大概是Simulink新手最经典的一个错误填4000意味着6.37Hz波形看起来完全不对。Switch模块的阈值。伯努利二进制发生器的输出是0和1Switch的阈值设在0.5当输入大于0.5输出上支路对应比特1的f1小于0.5输出下支路对应比特0的f2。这个逻辑是对的。但如果阈值设在2那就永远选择下支路调制输出全是8kHz的f2载波。可以先用Scope看调制输出的频率是否随比特变化这一步就能发现问题。变步长求解器的隐患。如果Simulink设置成变步长正弦波模块在非均匀采样下的输出看起来会有抖动Scope上显示的是折线而不是平滑正弦。虽然不影响最终统计结果但答辩时可以被人挑刺。我通常把求解器设置成固定步长、步长1e-5这样Scope显示更美观和.m代码的采样也完全对应。4.4 调试一个案例的过程记录有一次帮A同学排查他的误码率曲线在Eb/N08dB时依然有0.2明显不对。我先让他画出两条支路的低通滤波输出结果发现在一个比特1的码元内f2支路的输出在码元后半段逐渐上升导致抽样点误判。初步判断是f2带通滤波器的通带太宽把f1的一部分能量也滤进来了。检查发现他设置的f2带通滤波器为[7000 9000]Hz而f1是4kHz看似没有重叠。但进一步看他的f1带通滤波器通带是[2000 6000]Hz带宽4kHz2kHz范围内与f2支路的过渡带重叠。滤波器不是砖墙4kHz的能量在6kHz处仍有残余如果f2支路的通带边缘在7kHz起步有一定抑制但残余已经比噪声大了。把f1带通改成[3000 5500]Hzf2改成[6500 9000]Hz后误码率立即恢复到和理论曲线重合的水平。这就是滤波器带宽不是越大越好的实际案例。4.5 提高系统性能的扩展改造如果做完了基础版本想在报告里加亮点可以从以下几个方向扩展把相位不连续FSK改成CPFSK。实现方式是在调制端使用频率控制振荡器频率切换时保持相位连续比较前后两种方案在频谱上的带外衰减差异。实现非相干解调。替换掉本地载波相乘部分用包络检测取绝对值后低通或直接用Hilbert变换求包络对比相干和非相干的误码率曲线能直观看到性能差距。加入成形滤波。发送端加根升余弦滤波器观察带宽压缩和ISI变化这个方向能把2FSK和数字基带传输的知识打通。我实际做的时候把CPFSK和非相干解调都分别跑通了。三份结果放在报告里整套系统的理论-仿真-工程折衷逻辑非常完整。结尾这套系统做完之后我对调制解调这件事的理解从公式记忆变成了链路直觉。以前背下相干解调框图但不知道哪一级最重要自己写过一遍之后才明白抽样判决时刻和滤波器群延迟才是决定误码率高低的魔鬼细节。如果你也在做2FSK或者类似的数字调制课程设计我建议一定不要只去网上找一份代码跑通就交差试试改改载波频率、试试去掉带通滤波器、试试用非相干解调观察每个改动对最终误码率曲线的影响。折腾几次之后你收获的就不只是一个仿真项目而是对整个数字通信链路每一个环节在抗噪声中到底扮演什么角色的直观认知。另外一个很实用的经验是写完第一个版本之后把所有参数都整理成一个独立的参数脚本文件比如param_init.mSimulink模型也从这个脚本加载参数。这样无论是跑参数扫描、做对比实验还是写报告配图都会省下大量重复操作时间。这个习惯我在后来的很多仿真项目里都在用推荐你从现在就开始养成。