ARTICLE DETAIL

资讯详情

深耕网站视觉设计与运营推广的一线实战洞察。

基于MATLAB的数字调制解调系统仿真:从2ASK到QPSK全解析

基于MATLAB的数字调制解调系统仿真:从2ASK到QPSK全解析 搞通信或者做信号处理方向的人大概率都绕不过“数字调制解调”这个坎。不管是刚入门《通信原理》的学生还是准备课程设计、毕业设计的同学终归要面对一个问题书上讲的2ASK、2FSK、BPSK、QPSK到底怎么在工作站上真实跑起来之前带过几个师弟做这个题目也帮不少人看过代码发现很多人的坑都出在同一个地方——理论公式能默写但一让写MATLAB仿真就卡壳。这篇文章我就把这个“基于MATLAB的基本数字调制解调系统”彻底打开揉碎从设计思路、核心代码、误码率验证到高频报错一步不落给你讲清楚你照着敲就能调通。为什么要专门写这个因为很多课程设计、综合实验都会选这个题目——既能锻炼信号处理的基本功又不至于难到无从下手。但网上搜到的资料要么只给一个“能出图”的脚本要么是纯理论推导没有可执行代码两头不挨着。我在这里会给你一整套能直接运行的MATLAB实现并且解释每一步为什么要这么写参数为什么这么取遇到报错怎么排查。适合正在做通信原理课程设计、MATLAB仿真大作业、或者想快速构建数字调制仿真平台的读者。1. 整体设计思路先搭框架再填细节1.1 为什么选MATLAB而不是Python或Verilog先说结论做算法的验证和演示MATLAB依然是通信仿真场景下效率最高的选择没有之一。虽然Python的SciPy和commpy也能做但MATLAB的通信工具箱Communications Toolbox把很多底层函数都封装好了你不需要手写滤波器系数不需要自己实现卷积一条rcosdesign就能生成升余弦滤波器一个berawgn函数就能直接拿理论误码率做对照这在课程设计级别的项目中能省掉大量验证理论正确性的时间。但这里有个关键点很多教程会直接教你用comm.ASKModulator、comm.PSKModulator这些系统对象几行代码就搞定调制解调。不是不行而是对理解原理没有帮助也容易被老师怀疑是“调包的”。更好的做法是核心调制和解调部分用基础语法手写一遍——生成载波、相乘、判决——然后再用工具箱函数做交叉验证。这样既有学习的深度又有工程上的可靠度汇报的时候也讲得出细节。1.2 系统整体架构与模块划分一个完整的数字调制解调仿真系统逻辑上其实就五块信源、调制、加噪、解调、性能分析。信源负责产生二进制数据流这里用randi([0 1], 1, N)生成0/1序列调制部分根据要仿真的方式ASK/FSK/PSK把比特映射成对应的波形信道部分这里是仿真模型用AWGN加性高斯白噪声来模拟真实信道中的噪声叠加信噪比用SNR参数控制解调部分做的是逆过程——相关解调或者相干解调最后输出判决结果性能分析则是统计误比特率并且和理论的误码率曲线做对比。实际写代码的时候我建议你按功能拆成多个函数文件而不是全挤在一个脚本里。比如ask_mod.m、fsk_mod.m、psk_mod.m对应调制端ask_demod.m、fsk_demod.m、psk_demod.m对应解调端主脚本run_simulation.m负责调用和绘图。这样后续要增加调制方式比如加一个16QAM非常方便排查问题的时候也能直接锁定哪个函数出了问题。2. 三大基础调制方式的原理与MATLAB实现2.1 2ASK最直观的“开关键控”2ASK二进制振幅键控的本质就是用二进制数据去控制载波的幅度发送“1”的时候输出载波发送“0”的时候不输出载波。数学表达式写出来就是s_ASK(t) b(t) * Ac * cos(2πfct)其中b(t)就是映射后的电平序列通常把0映射为01映射为1。理解了这个原理代码就顺理成章了。给定比特序列bits、载波频率fc、采样率fs和每个比特的采样点数sps我们先要把数据扩展成和载波时间轴对应的波形序列。这里有一个新手经常踩的坑把“过采样”和“采样率”搞混。每个比特持续时间为Tb采样率是fs那么一个比特内就有fs*Tb个采样点这个值就是我们说的spssamples per symbol。仿真的时间轴t要按总的采样点数来生成而不是按比特索引来生成。代码如下function [s_ask, t] ask_mod(bits, fc, fs, sps) % 2ASK调制 % bits: 0/1比特序列 % fc: 载波频率 (Hz) % fs: 采样率 (Hz) % sps: 每个比特的采样点数 N length(bits); total_samples N * sps; t (0:total_samples-1) / fs; % 将比特序列扩展为波形电平上行到采样域 level repelem(bits, sps); % 生成载波并完成相乘 carrier cos(2*pi*fc*t); s_ask level .* carrier; end注意这里用了repelem函数它是把每个元素重复sps次一步就完成了零阶保持扩展。如果你使用的是老版本MATLAB2015a之前没有repelem可以用reshape(repmat(bits, sps, 1), 1, [])替代效果完全一样。载波频率fc的选取有个约束为了让仿真波形看得清楚fc至少要大于fs/(2*sps)且尽量保证每个比特内有整数个完整的载波周期。比如fs100e3、sps100那么比特率就是fs/sps1000bps载波频率取5e3到10e3都是合适的。取fc和sps满足整数倍关系还有一个额外好处解调端抽样判决时能落在载波波形的固定相位上不容易出现抽样点恰好在过零处的情况。解调端比较常用的是相干解调法。既然是课程设计级别的仿真我们就用理想同步的假设——即接收端知道载波的频率和相位直接用本地载波相乘再低通滤波。低通滤波器可以用FIR也可以用最简单的移动平均滤波。移动平均在码元速率远低于采样率时效果足够好而且不用调滤波器参数对于新手更友好。function bits_hat ask_demod(rx_signal, fc, fs, sps, threshold) % 2ASK相干解调 N_symbols length(rx_signal) / sps; t (0:length(rx_signal)-1) / fs; % 本地载波相乘 mix rx_signal .* cos(2*pi*fc*t); % 移动平均滤波窗口长度等于sps % 这样每个抽样点输出的是该比特内的平均能量 window ones(1, sps) / sps; filtered conv(mix, window, same); % 按比特中心位置抽样 sample_idx round(sps/2) : sps : length(filtered); sampled filtered(sample_idx); % 阈值判决对于OOK最佳判决门限约为幅度的一半 bits_hat sampled threshold; bits_hat double(bits_hat(:)); end这里conv用的是same参数保证输出长度和输入相同。移动平均的窗口设为sps作用是对每个比特内的采样点求平均等效于一个截止频率约为0.5/sps*fs的低通滤波器能把倍频分量滤掉。阈值判决这里理论上OOK的最佳判决门限是A/2A为接收信号幅度但在仿真中信噪比已知直接取接收到波形幅值包络的一半就可以。实际写的时候可以先用max(filtered)/2估计或者直接在调制的时候把1码的幅度规定为10码幅度为0那样阈值就可以固定为0.5。2.2 2FSK用频率承载信息2FSK二进制频移键控看名字就知道用两个不同的载波频率表达0和1。fc1对应比特1fc2对应比特0。实现上比ASK多一步频率切换但思路依然非常直白。function s_fsk fsk_mod(bits, fc1, fc2, fs, sps) % 2FSK调制非连续相位FSK N length(bits); total_samples N * sps; t (0:total_samples-1) / fs; % 生成两种频率的载波整段 carrier1 cos(2*pi*fc1*t); carrier2 cos(2*pi*fc2*t); % 初始化输出 s_fsk zeros(1, total_samples); % 对每个比特选择对应频段的载波 for k 1:N idx (k-1)*sps 1 : k*sps; if bits(k) 1 s_fsk(idx) carrier1(idx); else s_fsk(idx) carrier2(idx); end end end这里用的是“分段拼接”的方式而不是把两个载波全乘上一个0/1门控序列再相加那种方式会出现两个频率叠加的问题导致信号功率翻倍实现上更直观且是教科书上最标准的2FSK相位不连续模型。FSK的解调方法有两种思路一种是相干解调类似ASK那样分别和f1、f2相关比较两路能量大小另一种是非相干解调直接用包络检波后比较。在实际通信系统中非相干解调因为没有载波同步的要求往往更常用。在MATLAB仿真里我们可以用带通滤波器分成两路再求包络比较。但更简单的做法是用相关运算把接收信号分别和cos(2πf1t)、cos(2πf2t)相乘再积分哪个积分值大就判哪个。function bits_hat fsk_demod(rx_signal, fc1, fc2, fs, sps) % 2FSK非相关解调并行能量比较 N_symbols length(rx_signal) / sps; t (0:length(rx_signal)-1) / fs; % 与两个频率分别做相关 corr1 reshape(rx_signal, sps, N_symbols) .* reshape(cos(2*pi*fc1*t), sps, N_symbols); corr2 reshape(rx_signal, sps, N_symbols) .* reshape(cos(2*pi*fc2*t), sps, N_symbols); energy1 sum(corr1, 1); energy2 sum(corr2, 1); bits_hat energy1 energy2; bits_hat double(bits_hat(:)); end这里reshape把接收信号按比特切分然后一次性做内积运算避免循环。注意cos(2*pi*fc2*t)是整段的但reshape之后每一列正好是一个比特的载波片段所以内积结果就是该比特内信号和对应载波的相关系数。对于FSK两个频率的间隔有讲究如果fc2 - fc1等于比特率的整数倍那么两个频率在判决时刻是正交的误码性能最好。这也是为什么仿真中常取fc120e3、fc230e3、比特率10e3这类倍数关系——保证正交性。2.3 BPSK与QPSK相位调制的基础形态BPSK二进制相移键控是最基础的相位调制0和1分别用载波的0相位和π相位表示表达式为s_BPSK(t) A * cos(2πfct φk)其中φk为0或π。在MATLAB实现中有一个更简洁的做法——电平映射把0映射为-11映射为1然后直接乘载波function s_bpsk bpsk_mod(bits, fc, fs, sps) N length(bits); total_samples N * sps; t (0:total_samples-1) / fs; % BPSK: 0 - -1, 1 - 1也可以反着映射关键是差分 level repelem(bits, sps); level 2 * level - 1; % 把0/1变成-1/1 s_bpsk level .* cos(2*pi*fc*t); end解调端同样做相干解调乘本地载波后低通滤波然后抽样判决。这里的判决门限是0大于0判为1即原比特1小于0判为-1即原比特0。代码可以和ASK共用一个框架区别只在最后阈值判断等于0而不是0.5。QPSK则更进一步每2个比特映射成一个符号有四个相位状态π/4、3π/4、5π/4、7π/4或者偏移版本。映射表可以自己定比如格雷映射会让相邻符号只有一位不同误码性能更好。QPSK实现的关键是用reshape把比特流两两分组然后用查表或者公式映射成复信号再和载波相乘取实部。做QPSK最推荐的方式是使用复数基带等效而不是直接产生带通波形这样实现和理解都简单。function symbols qpsk_mod(bits) % QPSK调制基带等效格雷映射 % 输入bits长度为偶数 even bits(1:2:end); odd bits(2:2:end); % 00 - -1-j, 01 - -1j, 10 - 1-j, 11 - 1j % 注意这里要归一化功率 symbols ((even*2-1) 1j*(odd*2-1)) / sqrt(2); end这个基带模型后续如果要加频偏、相偏可以对symbols乘一个exp(1j*(2*pi*f_offset*t phase_offset))如果要生成实信号上变频再乘以载波取实部即可。因为QPSK的符号速率是比特率的一半所以基带仿真里时间和采样率的设置也要相应调整。3. 系统级联调与误码率性能分析3.1 加性高斯白噪声信道建模与信噪比换算有了调制端和解调端中间必须经过“信道”这步。实际仿真里AWGN信道是绝对的主流因为它数学上可解析很多通信系统的理论误码率都是在AWGN下推导的。MATLAB里加噪声可以用awgn函数但很多新手不注意信噪比的单位换算。awgn(x, snr, measured)里的snr单位是dB且是信号功率与噪声功率的比值和Eb/N0不同后者是每比特能量与噪声功率谱密度的比值。你可能需要的是直接以Eb/N0为横轴绘制误码率曲线这也是通信原理课程里标准画法。那就要手动换算对于BPSK/2ASK比特率和符号率相等对于QPSK符号率是比特率的一半。信号平均功率P_signal可以用mean(abs(signal).^2)计算但这里面还牵扯到过采样的问题。最简单的做法是用awgn把不同Eb/N0转成对应的SNR公式为SNR_dB EbN0_dB 10*log10(Rb/Bn)其中Rb为比特率Bn为噪声带宽仿真中等于采样率/2。转换后调用awgn(rx, SNR_dB, measured)就能精确控制信噪比。在课程设计报告中这个换算过程通常是必考/必写内容最好自己写个脚本验证几组数字。一个简便做法是在仿真时设置采样率fs1归一化此时一个符号对应一个采样点即sps1这种情况下SNR和Eb/N0的换算关系就变成了SNR_dB EbN0_dB 10*log10(Rb/fs)只要Rbfs两者直接相等。不搞过采样运行速度快到飞起很适合蒙特卡洛批量仿真。3.2 蒙特卡洛仿真与误码率曲线绘制蒙特卡洛仿真本质就是“跑很多次实验统计错误比例”。对每个Eb/N0点我们生成几万到几十万个比特经过调制、加噪、解调比较收发比特统计误码个数除以总比特数得到误码率。理论上仿真次数越多误码率越接近真实值。一般误码率低到1e-5时至少需要传输1e6比特才有把握否则统计误差太大。主脚本可以这样组织clear; close all; clc; EbN0_dB 0:2:12; N_bits 1e6; % 仿真比特数 sps 1; % 基带等效模型不用过采样 fc 0; % 等效载频0基带 ber_ask zeros(size(EbN0_dB)); ber_bpsk zeros(size(EbN0_dB)); ber_fsk zeros(size(EbN0_dB)); for idx 1:length(EbN0_dB) EbN0 10^(EbN0_dB(idx)/10); bits randi([0 1], 1, N_bits); % ASK s ask_mod_baseband(bits); % 基带等效返回幅度序列 noise_var 1 / (2*EbN0); noise sqrt(noise_var) * randn(1, N_bits); r s noise; bits_hat r 0.5; ber_ask(idx) sum(bits ~ bits_hat) / N_bits; % BPSK s 2*bits - 1; noise_var 1 / (2*EbN0); r s sqrt(noise_var)*randn(1, N_bits); bits_hat r 0; ber_bpsk(idx) sum(bits ~ bits_hat) / N_bits; % FSK正交非相干 s zeros(1, N_bits); s(bits1) 1; s(bits0) 1; % 能量归一化 % 非相干正交FSK理论BER: 0.5*exp(-0.5*EbN0) ber_fsk(idx) 0.5 * exp(-0.5*EbN0); end figure; semilogy(EbN0_dB, ber_ask, o-, LineWidth, 1.5); hold on; semilogy(EbN0_dB, ber_bpsk, s-, LineWidth, 1.5); semilogy(EbN0_dB, ber_fsk, ^-, LineWidth, 1.5); grid on; xlabel(Eb/N0 (dB)); ylabel(误码率 (BER)); legend(2ASK, BPSK, 2FSK(非相干理论));这段代码可能不是最规范的但帮你理清了蒙特卡洛的本质跑大量随机数据统计错误比例。实际做的时候我建议你写一个循环把调制方式也做成参数把不同调制方式的仿真曲线和理论曲线画在同一张图里能直观看到“仿真散点压在理论曲线上”的效果。理论误码率公式也需要提前准备好用在报告里2ASK相干解调P_b Q(sqrt(Eb/N0))BPSK相干解调P_b Q(sqrt(2*Eb/N0))2FSK相干解调P_b Q(sqrt(Eb/N0))2FSK非相干解调P_b 0.5*exp(-0.5*Eb/N0)QPSK与BPSK在相同Eb/N0下误码率相同给定格雷映射其中Q(x)0.5*erfc(x/sqrt(2))MATLAB里可以直接用qfunc。4. 实操过程中的高频报错与排错经验4.1 维度不匹配、索引越界与数据类型问题维度不匹配是MATLAB新手第一大问题几乎90%的报错都是Matrix dimensions must agree。出现原因多半是时间向量长度和信号向量长度对不上或者repelem之后的长度不是载波长度的整数倍。排查方法是在报错行之前用disp(length(t))、disp(length(signal))打印长度肉眼对一下。很多时候是因为sps设置导致bits长度与N*sps不一致比如赋值的索引范围(k-1)*sps1 : k*sps中的sps写错成了变量fs。索引越界的报错Index exceeds array bounds则往往出在循环里。比如对FSK解调时reshape(rx_signal, sps, N_symbols)要求length(rx_signal)整除sps。如果接收信号的长度因为卷积操作发生了变化用conv时的same和full差异就很容易越界。我的习惯是每次卷积、滤波之后都检查一下length确保它是sps的整数倍不行就手动截断到整数倍。还有一个数据类型的大坑MATLAB的ber统计结果如果做除法两个整数相除不会自动变成浮点数其实会但如果你使用bits_hat r 0.5得到的逻辑数组和原始bitsdouble数组做不等比较是合法的但直接求和后如果其中一个矩阵是逻辑型、另一个是double型有时会产生隐式转换导致结果不对。稳妥的做法是统一转成double再运算。4.2 采样率、码元速率与载波频率的取值技巧参数设置是仿真能不能“好看”的关键。如果载波频率太低比如比码元速率还低波形上根本无法区分出一个比特内的多个载波周期解调出来自然不对如果载波频率太高在每个比特内的采样点数又不够仿真精度不够。经验法则是载波频率取码元速率的5到10倍且一个比特内的采样点数至少要有20个点中频采样或者100个点显示波形用。具体来说如果比特率Rb1000bps、载波fc10kHz、每个比特采样100点那么采样率fs100kHz。载波每周期有fs/fc10个采样点足够平滑每个比特里有10个完整载波周期视觉上能明显看到“每个比特内有10个正弦波”。这是展示课设波形最合适的参数组合。如果是基带等效模型跑性能仿真就不需要这些约束直接sps1效率最高。还有一个小细节用plot画调制波形时如果数据点太多图形会糊成一团可以只画出前几十个比特对应的片段同时把Marker设置成.之类的让波形更清晰。用stem画比特序列时也一样取前20个比特足够了。4.3 误码率曲线高频波动或“下不去”的排查方法很多同学跑出来的误码率曲线在高信噪比时出现平台期不再下降或者抖动非常厉害。出现平台期通常不是调制解调代码的随机噪声问题而是“残余的固定干扰”——常见原因有三个一是滤波不够干净相干解调后信号里还残留二倍频分量但抽样恰好抽到了残留较大的位置二是本地载波和发送端相位没有完全同步在相干解调里哪怕差一点相位都会让有效信号幅度衰减三是阈值不随信噪比调整在高信噪比时信号幅度可能因滤波器暂态而畸变固定阈值不再最优。抖动厉害则多半是仿真比特数不够。比如在Eb/N010dB时误码率约1e-5跑1e5比特理论上只能期望1个误码这时的点毫无统计意义。我建议写成自适应循环误码率低于某个阈值就自动加大仿真比特数直到误码数达到至少50~100个这样曲线才平滑。还有一个经验技巧不要在每个信噪比点从零重新生成随机种子比如用rng(idx)让每个点固定不同的种子这样跑出来的曲线更平滑复现性也好。5. 工具箱函数对照与进阶扩展方向5.1 用Communications Toolbox做交叉验证前面我强调用手写代码理解原理但工程验证阶段完全可以用工具箱函数来提高效率。比如用comm.PSKModulator、comm.PSKDemodulator、comm.AWGNChannel这些系统对象做一遍完整的蒙特卡洛仿真当参考曲线。如果手写代码的曲线和工具箱的曲线对不上就说明手写实现里有bug。用工具箱做BPSK加噪声的代码片段很短M 2; pskMod comm.BPSKModulator; pskDemod comm.BPSKDemodulator; channel comm.AWGNChannel(NoiseMethod, Signal to noise ratio (SNR), SNR, 10); errorRate comm.ErrorRate; txSig pskMod(bits); rxSig channel(txSig); demodBits pskDemod(rxSig); errStats errorRate(bits, demodBits);注意pskm的输入要按列传入比特也必须是列向量这是工具箱和手写代码习惯的差异。交叉验证的思路是先把工具箱跑通记录结果再把手写代码的结果叠在同一个图上如果偏差在0.5dB以内就可以认为实现正确。课程设计的评审老师看到你既能手写、又能用工具箱交叉验证印象分会很加分。5.2 从单载波扩展到多载波与成型滤波如果你学有余力做完基本调制解调系统后可以往两个方向扩展一是加成型滤波。把RZ归零的方波脉冲改用升余弦脉冲rcosdesign能明显降低频谱旁瓣和码间串扰ISI。加入成型滤波之后收发两端要用同样的滤波器做匹配滤波整个系统的抗噪声性能会更好眼图也会变得清晰。这是从“能工作”到“像回事”的关键一步。二是扩展到调制阶数更高的方式比如8PSK、16QAM。16QAM的实现思路和QPSK差别不大只是星座图上的映射点变多了判决区域更复杂。做这些扩展时你会发现星座图scatterplot是特别好的调试工具它能一次性告诉你相位偏了多少、幅度有没有归一化、噪声方差多大——比光看BER数字直观得多。课程设计里加一张漂亮的星座图比对报告质感瞬间上来。我个人在做这些仿真时最大的体会是通信仿真里的坑大多是一些“微观”的维度、索引问题而不是理论问题。理论和代码之间隔着的一层正是这些细节。比如本地的载波要不要带初相、滤波器的延迟怎么补偿、蒙特卡洛的帧怎么分组每一样都会影响最后的曲线。把这些细节一个个踩平你对整个系统的理解才算真正到位。最后再分享一个小习惯写仿真脚本时每个模块用单独的cell块两个百分号%%开头分隔跑完一段就CtrlEnter执行一段。这样一旦结果不对可以直接定位到是从哪个环节开始歪掉的。比起一次性写完整个脚本再痛苦地debug这种边写边验的方式会舒服得多。希望这篇文章能帮你把MATLAB数字调制解调仿真一把跑通有卡壳的地方欢迎来聊。
返回列表