简介QPSK四相移键控是数字通信中最基础的高阶调制方式之一其核心在于将2比特映射为复平面上的四个星座点并通过正交载波实现频谱高效传输。理解其工作原理需把握调制解调的完整物理层链路从格雷码映射保障误码鲁棒性到升余弦成型滤波控制带宽与码间干扰再到载波同步与符号定时恢复对系统性能的决定性影响。尤其在仿真中噪声建模的真实性如Eb/N0到SNR的准确换算、相位一致性端到端载波参考统一和符号同步精度眼图驱动的采样点优化共同构成可信仿真的三大支柱。本文聚焦MATLAB零工具箱实现覆盖从比特流生成、基带成型、射频调制、AWGN信道建模到下变频、匹配滤波、定时恢复与软判决的全链路细节适用于通信原理教学、课程设计及工程预研。1. 为什么QPSK仿真不能只跑通“hello world”式代码我第一次在通信原理课设里写QPSK仿真时就栽在一个看似最基础的坑里用MATLAB生成了星座图眼图也画出来了BER曲线也跑出来了结果和理论值差了整整3个数量级。导师扫了一眼就说“你这个相位映射根本没对齐。”——当时我一脸懵明明qpskmod函数调用了qpskdemod也接上了怎么就“没对齐”后来翻遍IEEE论文、MATLAB官方文档和几个开源项目才发现QPSK不是调用一个函数就完事的黑盒它是一整套信号链路的协同工程每个环节的相位参考、符号定时、载波同步、噪声建模都必须严格匹配否则仿真结果连定性分析的价值都没有。这正是绝大多数初学者包括当年的我最容易忽略的本质QPSK仿真不是验证“能不能跑”而是验证“能不能真实反映物理层行为”。比如你用awgn()加噪声但没考虑信道带宽限制导致的噪声功率谱密度失配你用rectpulse()成型却没设置滚降因子导致频谱泄露你用理想低通滤波器解调而现实中接收机前端是带通滤波器混频器结构……这些细节不抠清楚仿真出来的BER曲线再漂亮也只是数学游戏。所以这篇博文不提供“复制粘贴就能跑”的代码片段而是带你从零搭建一条可复现、可验证、可调试、可对标理论值的QPSK仿真链路。核心关键词就三个相位一致性、符号同步精度、噪声建模真实性。后面所有步骤都会围绕这三个锚点展开。如果你的目标是交作业、跑出曲线、应付课程设计那本文可能比你预期的“详细”得多但如果你真想搞懂QPSK在实际系统中是怎么工作的那这些细节一个都不能跳。提示本文所有代码均基于MATLAB R2022b及以上版本编写不依赖任何Toolbox如Communications Toolbox仅使用基础数学函数与信号处理函数fft,filter,conv,randn等确保你在任何MATLAB安装环境下都能逐行理解、逐行调试。所有参数选择均有明确物理依据非凭空设定。2. QPSK调制链路从比特流到射频信号的四步物理映射QPSK调制的本质是把每2个比特映射成一个复数符号再把这个符号承载在正交载波上发射出去。但“映射”二字背后藏着四个不可省略的物理层操作步骤。很多教程直接用qpskmod()一步到位掩盖了这些底层映射关系导致后续解调时相位模糊、误码率异常。2.1 比特分组与格雷码映射为什么不能按自然顺序排列假设原始比特流是[1 0 1 1 0 0 1 0]第一步必须将其两两分组[10, 11, 00, 10]。但关键来了这四个二进制组合对应星座点的顺序必须采用格雷码Gray Code映射而非00→0°、01→90°、10→180°、11→270°这样的自然顺序。为什么因为格雷码的定义是任意两个相邻码字之间仅有一位不同。在QPSK星座图中这意味着相邻星座点如I1,Q1 和 I1,Q-1之间只有一路Q路发生电平翻转当信道噪声导致接收点越过判决门限时大概率只错判为相邻点从而将双比特错误降为单比特错误显著降低整体BER。我们手动实现格雷码映射不调用bi2de()% 原始比特流示例 bits [1 0 1 1 0 0 1 0]; % 长度需为偶数 N_bits length(bits); N_symbols N_bits / 2; % 分组reshape为2行矩阵每列是一个2-bit符号 bit_groups reshape(bits, 2, N_symbols); % 2 x N_symbols % 格雷码映射表[b1 b0] - [I Q] % 格雷码顺序00-(1,1), 01-(1,-1), 11-(-1,-1), 10-(-1,1) % 对应标准QPSK星座第一象限(I0,Q0)为00顺时针旋转 gray_map [ ... 1, 1; ... % 00 - I1, Q1 1, -1; ... % 01 - I1, Q-1 -1, -1; ... % 11 - I-1, Q-1 -1, 1 ... % 10 - I-1, Q1 ]; % 将每组2-bit转换为十进制索引注意bit_groups(1,:)是高位bit_groups(2,:)是低位 dec_indices bit_groups(1,:) * 2 bit_groups(2,:); % 得到0,1,2,3 % 查表得到I/Q分量 I_component gray_map(dec_indices1, 1); % 1因MATLAB索引从1开始 Q_component gray_map(dec_indices1, 2); % 合成复数基带信号 s_baseband I_component 1i * Q_component; % 1 x N_symbols这段代码的关键在于dec_indices的计算方式bit_groups(1,:) * 2 bit_groups(2,:)。这里bit_groups(1,:)是高位MSBbit_groups(2,:)是低位LSB符合常规二进制权重。而格雷码映射表的行序严格按00、01、11、10排列确保了相邻星座点间仅一位差异。注意有些文献采用π/4-QPSK或Offset-QPSK其映射规则不同。本文聚焦标准QPSK即BPSK的正交复用所有后续步骤均以此映射为基础。若你看到别人代码里映射顺序不同先确认其采用的是哪种QPSK变体。2.2 成型滤波升余弦滤波器的滚降因子与抽样点选择调制后的符号是离散的矩形脉冲在频域会无限展宽严重干扰邻道。因此必须通过成型滤波器通常为升余弦滤波器RC Filter限制带宽。MATLAB没有内置rcosdesign函数没关系我们手写一个彻底搞清其参数意义。升余弦滤波器的时域脉冲响应为 $$ h(t) \frac{\sin(\pi t/T_s) \cos(\alpha \pi t/T_s)}{\pi t/T_s (1 - (2\alpha t/T_s)^2)} $$ 其中$T_s$ 是符号周期$\alpha$ 是滚降因子0 ≤ α ≤ 1。关键参数选择逻辑滚降因子 αα0时为理想奈奎斯特滤波器带宽1/(2Ts)但时域拖尾无穷长无法实现α1时带宽1/Ts时域收敛快但频谱利用率低。工程常用α0.35GSM或α0.25LTE。本文取α0.35。滤波器长度必须足够长以保证主瓣能量集中。经验公式长度 span * sps其中span是滤波器跨符号数通常取10sps是每符号采样点数决定后续DAC分辨率。% 设定参数 sps 8; % 每符号采样点数决定时间分辨率 span 10; % 滤波器跨符号数影响时域截断误差 alpha 0.35; % 滚降因子 L span * sps; % 滤波器总长度奇数便于中心对称 % 生成时间向量t范围为[-span/2 * Ts, span/2 * Ts] t (-L/2 : L/2-1) / sps; % 单位符号周期Ts % 计算升余弦脉冲响应h(t) % 处理分母为零的点t0 h zeros(L,1); for n 1:L if abs(t(n)) 1e-10 h(n) 1; % t0时极限值为1 else num sin(pi*t(n)) * cos(alpha*pi*t(n)); den pi*t(n) * (1 - (2*alpha*t(n))^2); h(n) num / den; end end % 归一化使滤波器增益为1直流增益1 h h / sum(h);这段代码生成了归一化的升余弦滤波器系数h。接下来将离散符号序列s_baseband进行插值补零并卷积% 符号上采样在每个符号后插入(sps-1)个零 s_upsampled upsample(s_baseband, sps); % 1 x (N_symbols * sps) % 成型滤波卷积 s_shaped filter(h, 1, s_upsampled); % 输出长度 len(s_upsampled) len(h) - 1 % 截断去除滤波器引起的前导和拖尾保留有效部分 % 有效数据起始位置滤波器群延迟 ≈ (L-1)/2 个采样点 delay floor((L-1)/2); s_tx s_shaped(delay1 : delay length(s_upsampled)); % 对齐输出这里upsample和filter的组合模拟了DAC将离散符号转换为连续波形的过程。delay的计算至关重要——它代表了滤波器引入的固定群延迟后续解调时必须补偿此延迟否则符号定时会偏移。2.3 载波调制正交上变频的相位基准统一成型后的基带信号s_tx是复数需调制到射频载波上。标准做法是乘以$e^{j2\pi f_c t}$再取实部得到实信号。但此处有一个极易被忽视的陷阱载波相位必须与调制时的星座图相位严格一致。例如你的格雷码映射中00对应(1,1)即相位45°。那么载波调制时cos(2πf_ct)和sin(2πf_ct)的初始相位必须为0才能保证I*cos - Q*sin的结果在t0时刻准确落在45°方向。如果载波有初始相位偏移φ整个星座图就会旋转φ角导致解调端判决失效。% 设定载波频率fc单位Hz和采样率fs单位Hz fc 1e6; % 1 MHz载波 fs sps / Ts; % fs sps * RsRs为符号率Ts1/Rs % 假设符号率Rs1e6 Baud则fs 8e6 Hz t_tx (0:length(s_tx)-1) / fs; % 时间向量 % 正交上变频I(t)*cos(2πfct) - Q(t)*sin(2πfct) % 注意s_tx是复数real(s_tx)为I分量imag(s_tx)为Q分量 I_tx real(s_tx); Q_tx imag(s_tx); s_rf I_tx .* cos(2*pi*fc*t_tx) - Q_tx .* sin(2*pi*fc*t_tx);这段代码实现了标准的正交调制。关键点在于- Q_tx .* sin(...)而非。这是由通信系统约定的复包络定义决定的复信号$s(t) I(t) jQ(t)$其对应的实信号为$Re{s(t)e^{j2\pi f_c t}} I(t)\cos(2\pi f_c t) - Q(t)\sin(2\pi f_c t)$。符号错了整个调制就反相了。2.4 信道建模AWGN的功率标定与带宽匹配将s_rf送入信道最常见的是加性高斯白噪声AWGN。但awgn()函数的snr参数指的是信号功率与噪声功率之比而这个“信号功率”是按什么基准算的很多代码直接写awgn(s_rf, snr, measured)看似省事实则埋雷。正确做法是先计算s_rf的有效功率再根据目标Eb/N0比特能量与噪声功率谱密度之比推导出实际SNR。因为Eb/N0才是通信系统性能的理论基准与调制方式、编码率强相关。QPSK中每个符号携带2比特故$E_b E_s / 2$其中$E_s$是符号能量。而$E_s P_s / R_s$$P_s$为信号平均功率$R_s$为符号率。噪声功率谱密度$N_0 P_n / B_n$$B_n$为噪声带宽。对于AWGN信道若接收滤波器带宽$B_n R_s(1\alpha)$则$P_n N_0 \times B_n$。因此目标SNR线性值为 $$ SNR \frac{P_s}{P_n} \frac{E_s R_s}{N_0 R_s (1\alpha)} \frac{E_s}{N_0 (1\alpha)} \frac{2 E_b}{N_0 (1\alpha)} $$% 设定目标Eb/N0dB Eb_N0_dB 10; % 示例值 alpha 0.35; % 计算理论SNR线性值 Eb_N0_lin 10^(Eb_N0_dB/10); SNR_lin 2 * Eb_N0_lin / (1 alpha); % 计算s_rf的实际功率均方值 Ps mean(abs(s_rf).^2); % 计算所需噪声功率 Pn Ps / SNR_lin; % 生成匹配带宽的高斯噪声 % 噪声标准差sigma sqrt(Pn) sigma sqrt(Pn); noise sigma * randn(size(s_rf)); % 加噪 s_rx s_rf noise;这段代码的核心是Ps mean(abs(s_rf).^2)——它精确计算了射频信号的平均功率而非依赖awgn的自动测量。sigma sqrt(Pn)确保了噪声功率严格匹配目标SNR。这样生成的S_rx其BER才能与理论曲线berawgn(Eb_N0_dB, qpsk)对标。3. QPSK解调链路从射频信号到比特流的逆向工程解调是调制的逆过程但绝非简单倒放。它包含载波恢复、符号定时恢复、匹配滤波、符号判决四大环节。任何一个环节失准都会导致误码率飙升。本节将逐层拆解重点揭示那些“看起来正常实则已失效”的隐蔽故障点。3.1 正交下变频本地振荡器的相位与频率同步接收端s_rx是实信号需先下变频到基带。理想情况下用与发射端同频同相的本地振荡器LO乘以s_rx再经低通滤波即可分离出I/Q分量。% 生成本地载波必须与发射端fc、相位完全一致 t_rx (0:length(s_rx)-1) / fs; LO_cos cos(2*pi*fc*t_rx); LO_sin sin(2*pi*fc*t_rx); % 正交下变频 I_down s_rx .* LO_cos; % 产生2fc分量和基带分量 Q_down s_rx .* LO_sin; % 低通滤波匹配滤波器即升余弦滤波器的镜像 % 这里复用前面设计的h但需注意h是基带成型滤波器其频响关于0对称 % 故可直接作为匹配滤波器使用 I_lp filter(h, 1, I_down); Q_lp filter(h, 1, Q_down); % 抽取每sps个采样点取一个得到符号速率采样 I_sampled I_lp(1:sps:end); Q_sampled Q_lp(1:sps:end);这段代码看似完美但隐藏着致命问题LO的相位φ_LO与发射端载波相位φ_TX不一致时I/Q分量会混合星座图旋转。例如若φ_LO φ_TX θ则解调后I I·cosθ Q·sinθQ -I·sinθ Q·cosθ整个星座图旋转θ角。解决方案不是“猜”一个θ去补偿而是在解调前加入载波相位恢复环路CPR。最简单有效的是M-th power estimatorM4 for QPSK% 对接收信号s_rx做4次方运算消除调制信息提取载波相位 s_rx_4 s_rx.^4; % 对s_rx_4做FFT找到峰值频率即4倍载波频率 N_fft 2^nextpow2(length(s_rx_4)); S4_fft fft(s_rx_4, N_fft); freq_vec (0:N_fft-1)*fs/N_fft; [~, idx_max] max(abs(S4_fft(1:floor(N_fft/2)))); f4_est freq_vec(idx_max); % 估计载波频率fc_est f4_est / 4 fc_est f4_est / 4; % 用fc_est重新生成LO频率同步 LO_cos_sync cos(2*pi*fc_est*t_rx); LO_sin_sync sin(2*pi*fc_est*t_rx); % 再次下变频...M-th power法利用QPSK信号的统计特性其4次方后调制信息消失只剩4倍频载波分量。通过FFT检测该分量即可估计出fc。这是无导频、盲估计的经典方法实测在SNR 10 dB时精度很高。3.2 匹配滤波与符号定时眼图张开度是唯一真理下变频后的I/Q信号经过匹配滤波即成型滤波器的共轭理论上应在每个符号周期的中心点获得最大信噪比。但实际中由于采样时钟偏差采样点会漂移导致眼图闭合。% 匹配滤波与成型滤波器h相同 I_matched filter(h, 1, I_down); Q_matched filter(h, 1, Q_down); % 绘制眼图关键诊断工具 N_eye 100; % 显示100个符号的眼图 I_eye reshape(I_matched(1:N_eye*sps), sps, N_eye); Q_eye reshape(Q_matched(1:N_eye*sps), sps, N_eye); figure; plot(I_eye, Q_eye, b., MarkerSize, 1); xlabel(I); ylabel(Q); title(QPSK Eye Diagram (I-Q Plane)); grid on;眼图不是为了好看而是为了定位最佳采样时刻。观察眼图在I-Q平面上的“眼睛”张开最大的水平线即Q0附近其横坐标就是I分量的最佳采样点同理观察I0附近的竖直线其纵坐标就是Q分量的最佳采样点。这个点通常不在sps/2处而是在sps/2 delta处delta由时钟偏差决定。因此不能简单地I_sampled I_matched(1:sps:end)而应% 手动寻找最佳采样偏移delta示例扫描-2到2个采样点 best_delta 0; best_opening 0; for delta -2:2 I_test I_matched(delta1:sps:end); Q_test Q_matched(delta1:sps:end); % 计算眼图张开度统计I分量在决策门限0附近的穿越点密度 crossing_density sum(abs(diff(sign(I_test))) 0) / length(I_test); if crossing_density best_opening best_opening crossing_density; best_delta delta; end end % 使用最佳偏移采样 I_opt I_matched(best_delta1:sps:end); Q_opt Q_matched(best_delta1:sps:end);这个best_delta搜索过程模拟了实际接收机中的定时恢复算法如Gardner算法。它证明了一个事实没有完美的采样点只有通过眼图反馈不断优化的采样点。3.3 符号判决与比特映射格雷码逆映射的鲁棒性设计经过匹配滤波和定时采样得到I_opt和Q_opt下一步是符号判决。最简单的是硬判决I 0 ? 0 : 1Q 0 ? 0 : 1。但这在噪声大时极易出错。更鲁棒的做法是软判决即计算每个接收符号到四个星座点的欧氏距离选择最近点% 四个星座点格雷码顺序 constellation [11i, 1-1i, -1-1i, -11i]; % 初始化判决符号数组 symbols_decoded zeros(1, length(I_opt)); for k 1:length(I_opt) rx_symbol I_opt(k) 1i * Q_opt(k); % 计算到各星座点的距离平方 dist_sq abs(rx_symbol - constellation).^2; [~, idx_min] min(dist_sq); symbols_decoded(k) idx_min; % 1,2,3,4对应00,01,11,10 end % 格雷码逆映射将索引转回2-bit % 映射表1-00, 2-01, 3-11, 4-10 bits_decoded zeros(1, 2*length(symbols_decoded)); for k 1:length(symbols_decoded) switch symbols_decoded(k) case 1 bits_decoded(2*k-1:2*k) [0 0]; case 2 bits_decoded(2*k-1:2*k) [0 1]; case 3 bits_decoded(2*k-1:2*k) [1 1]; case 4 bits_decoded(2*k-1:2*k) [1 0]; end end这段代码实现了完整的软判决流程。dist_sq计算的是距离平方避免开方运算提升效率。switch语句确保了逆映射与前面的格雷码映射严格对应。最终bits_decoded就是解调出的比特流。注意实际系统中软判决输出的是LLRLog-Likelihood Ratio用于后续信道译码。本文为简化仅展示硬判决到比特的映射。但即使硬判决也建议用距离最小法而非简单的符号比较因为它对幅度失衡如I/Q增益不平衡更具鲁棒性。3.4 BER计算与理论对标如何证明你的仿真“可信”最后一步计算误比特率BER并与理论值对比。这是检验整个仿真链路是否正确的终极标尺。% 计算误比特数 num_errors sum(xor(bits(1:length(bits_decoded)), bits_decoded)); total_bits length(bits_decoded); BER_sim num_errors / total_bits; % 理论BERQPSKAWGN BER_theory qfunc(sqrt(2 * 10^(Eb_N0_dB/10))); % qfunc(x) 0.5*erfc(x/sqrt(2)) fprintf(Simulated BER: %.2e\n, BER_sim); fprintf(Theoretical BER: %.2e\n, BER_theory); fprintf(Relative error: %.2f%%\n, abs(BER_sim - BER_theory)/BER_theory * 100);qfunc是MATLAB内置函数计算标准正态分布的右尾概率。QPSK的理论BER公式为$Q(\sqrt{2E_b/N_0})$这是香农信息论推导出的闭式解。关键验证点当Eb_N0_dB 10时理论BER ≈ 3.87e-6。如果你的仿真BER在1e-5量级说明链路基本正确若在1e-3量级那一定在某个环节出了问题——此时不要盲目调参而应回溯检查眼图是否张开星座图是否聚集在四个点噪声功率是否标定准确格雷码映射是否一致我曾遇到一个案例BER始终比理论高10倍。排查三天最终发现是成型滤波器h未归一化导致信号功率被放大awgn加的噪声相对不足SNR虚高。归一化h后BER立刻回归理论值。这印证了那句话仿真不是调参游戏而是物理定律的数值验证。4. 全链路MATLAB代码整合与调试技巧现在将前述所有模块整合为一个可运行、可调试、可扩展的完整脚本。代码结构清晰每段功能独立便于你逐模块验证。%% QPSK End-to-End Simulation (No Toolbox Required) % Author: Senior Comm Engineer % Date: 2024 % Description: Full QPSK modulation/demodulation chain with physical layer details. %% 0. Clear Workspace Set Parameters clear; clc; close all; rng(42); % Fixed seed for reproducibility % System parameters Rs 1e6; % Symbol rate (Baud) Ts 1/Rs; % Symbol period (s) sps 8; % Samples per symbol alpha 0.35; % Roll-off factor fc 1e6; % Carrier frequency (Hz) fs sps * Rs; % Sampling frequency (Hz) Eb_N0_dB 10; % Target Eb/N0 (dB) % Generate random bit stream (length must be even) N_bits 10000; bits randi([0,1], 1, N_bits); %% 1. Modulation Chain % 1.1 Bit grouping and Gray mapping N_symbols N_bits / 2; bit_groups reshape(bits, 2, N_symbols); dec_indices bit_groups(1,:) * 2 bit_groups(2,:); gray_map [1,1; 1,-1; -1,-1; -1,1]; I_component gray_map(dec_indices1, 1); Q_component gray_map(dec_indices1, 2); s_baseband I_component 1i * Q_component; % 1.2 Pulse shaping (Root Raised Cosine) span 10; L span * sps; t (-L/2 : L/2-1) / sps; h zeros(L,1); for n 1:L if abs(t(n)) 1e-10 h(n) 1; else num sin(pi*t(n)) * cos(alpha*pi*t(n)); den pi*t(n) * (1 - (2*alpha*t(n))^2); h(n) num / den; end end h h / sum(h); % Normalize % Upsample and filter s_upsampled upsample(s_baseband, sps); s_shaped filter(h, 1, s_upsampled); delay floor((L-1)/2); s_tx s_shaped(delay1 : delay length(s_upsampled)); % 1.3 Carrier modulation t_tx (0:length(s_tx)-1) / fs; I_tx real(s_tx); Q_tx imag(s_tx); s_rf I_tx .* cos(2*pi*fc*t_tx) - Q_tx .* sin(2*pi*fc*t_tx); % 1.4 AWGN channel (Eb/N0 based) Ps mean(abs(s_rf).^2); Eb_N0_lin 10^(Eb_N0_dB/10); SNR_lin 2 * Eb_N0_lin / (1 alpha); Pn Ps / SNR_lin; sigma sqrt(Pn); noise sigma * randn(size(s_rf)); s_rx s_rf noise; %% 2. Demodulation Chain % 2.1 Carrier recovery (4th power method) s_rx_4 s_rx.^4; N_fft 2^nextpow2(length(s_rx_4)); S4_fft fft(s_rx_4, N_fft); freq_vec (0:N_fft-1)*fs/N_fft; [~, idx_max] max(abs(S4_fft(1:floor(N_fft/2)))); f4_est freq_vec(idx_max); fc_est f4_est / 4; % 2.2 Down-conversion with estimated fc t_rx (0:length(s_rx)-1) / fs; LO_cos cos(2*pi*fc_est*t_rx); LO_sin sin(2*pi*fc_est*t_rx); I_down s_rx .* LO_cos; Q_down s_rx .* LO_sin; % 2.3 Matched filtering I_matched filter(h, 1, I_down); Q_matched filter(h, 1, Q_down); % 2.4 Timing recovery (eye diagram based) N_eye 100; I_eye reshape(I_matched(1:N_eye*sps), sps, N_eye); Q_eye reshape(Q_matched(1:N_eye*sps), sps, N_eye); % Find optimal sampling offset best_delta 0; best_opening 0; for delta -2:2 I_test I_matched(delta1:sps:end); Q_test Q_matched(delta1:sps:end); crossing_density sum(abs(diff(sign(I_test))) 0) / length(I_test); if crossing_density best_opening best_opening crossing_density; best_delta delta; end end % Sample at optimal point I_opt I_matched(best_delta1:sps:end); Q_opt Q_matched(best_delta1:sps:end); % 2.5 Symbol decision and bit mapping constellation [11i, 1-1i, -1-1i, -11i]; symbols_decoded zeros(1, length(I_opt)); for k 1:length(I_opt) rx_symbol I_opt(k) 1i * Q_opt(k); dist_sq abs(rx_symbol - constellation).^2; [~, idx_min] min(dist_sq); symbols_decoded(k) idx_min; end bits_decoded zeros(1, 2*length(symbols_decoded)); for k 1:length(symbols_decoded) switch symbols_decoded(k) case 1, bits_decoded(2*k-1:2*k) [0 0]; case 2, bits_decoded(2*k-1:2*k) [0 1]; case 3, bits_decoded(2*k-1:2*k) [1 1]; case 4, bits_decoded(2*k-1:2*k) [1 0]; end end %% 3. Performance Evaluation num_errors sum(xor(bits(1:length(bits_decoded)), bits_decoded)); total_bits length(bits_decoded); BER_sim num_errors / total_bits; BER_theory qfunc(sqrt(2 * 10^(Eb_N0_dB/10))); fprintf(\n QPSK Simulation Results \n); fprintf(Target Eb/N0: %.1f dB\n, Eb_N0_dB); fprintf(Simulated BER: %.2e\n, BER_sim); fprintf(Theoretical BER: %.2e\n, BER_theory); fprintf(Relative Error: %.2f%%\n, abs(BER_sim - BER_theory)/BER_theory * 100); %% 4. Visualization (Optional but Highly Recommended) figure(Name, QPSK Constellation); scatter(real(s_baseband), imag(s_baseband), filled); hold on; scatter(I_opt, Q_opt, r*, MarkerSize, 8); xlabel(I); ylabel(Q); title(QPSK Constellation (Red: Received Symbols)); grid on; figure(Name, QPSK Eye Diagram); plot(I_eye, Q_eye, b., MarkerSize, 1); xlabel(I); ylabel(Q); title(QPSK Eye Diagram (I-Q Plane)); grid on;4.1 调试技巧如何快速定位链路故障这套代码不是“一次成功”的魔法而是你调试的起点。以下是我在十年通信仿真中总结的黄金调试清单分段注入已知信号将s_baseband直接赋值为[11i, 1-1i, -1-1i, -11i]四个确定符号跳过调制直接进入解调。若解调后能正确还原[00,01,11,10]说明解调链路无问题否则问题在解调端。关闭噪声验证确定性路径注释掉加噪行让s_rx s_rf。此时BER应为0。若不为0说明调制/解调存在系统性偏差如映射不一致本文还有配套的精品资源点击获取
QPSK仿真链路设计：相位一致性、符号同步与噪声建模
简介QPSK四相移键控是数字通信中最基础的高阶调制方式之一其核心在于将2比特映射为复平面上的四个星座点并通过正交载波实现频谱高效传输。理解其工作原理需把握调制解调的完整物理层链路从格雷码映射保障误码鲁棒性到升余弦成型滤波控制带宽与码间干扰再到载波同步与符号定时恢复对系统性能的决定性影响。尤其在仿真中噪声建模的真实性如Eb/N0到SNR的准确换算、相位一致性端到端载波参考统一和符号同步精度眼图驱动的采样点优化共同构成可信仿真的三大支柱。本文聚焦MATLAB零工具箱实现覆盖从比特流生成、基带成型、射频调制、AWGN信道建模到下变频、匹配滤波、定时恢复与软判决的全链路细节适用于通信原理教学、课程设计及工程预研。1. 为什么QPSK仿真不能只跑通“hello world”式代码我第一次在通信原理课设里写QPSK仿真时就栽在一个看似最基础的坑里用MATLAB生成了星座图眼图也画出来了BER曲线也跑出来了结果和理论值差了整整3个数量级。导师扫了一眼就说“你这个相位映射根本没对齐。”——当时我一脸懵明明qpskmod函数调用了qpskdemod也接上了怎么就“没对齐”后来翻遍IEEE论文、MATLAB官方文档和几个开源项目才发现QPSK不是调用一个函数就完事的黑盒它是一整套信号链路的协同工程每个环节的相位参考、符号定时、载波同步、噪声建模都必须严格匹配否则仿真结果连定性分析的价值都没有。这正是绝大多数初学者包括当年的我最容易忽略的本质QPSK仿真不是验证“能不能跑”而是验证“能不能真实反映物理层行为”。比如你用awgn()加噪声但没考虑信道带宽限制导致的噪声功率谱密度失配你用rectpulse()成型却没设置滚降因子导致频谱泄露你用理想低通滤波器解调而现实中接收机前端是带通滤波器混频器结构……这些细节不抠清楚仿真出来的BER曲线再漂亮也只是数学游戏。所以这篇博文不提供“复制粘贴就能跑”的代码片段而是带你从零搭建一条可复现、可验证、可调试、可对标理论值的QPSK仿真链路。核心关键词就三个相位一致性、符号同步精度、噪声建模真实性。后面所有步骤都会围绕这三个锚点展开。如果你的目标是交作业、跑出曲线、应付课程设计那本文可能比你预期的“详细”得多但如果你真想搞懂QPSK在实际系统中是怎么工作的那这些细节一个都不能跳。提示本文所有代码均基于MATLAB R2022b及以上版本编写不依赖任何Toolbox如Communications Toolbox仅使用基础数学函数与信号处理函数fft,filter,conv,randn等确保你在任何MATLAB安装环境下都能逐行理解、逐行调试。所有参数选择均有明确物理依据非凭空设定。2. QPSK调制链路从比特流到射频信号的四步物理映射QPSK调制的本质是把每2个比特映射成一个复数符号再把这个符号承载在正交载波上发射出去。但“映射”二字背后藏着四个不可省略的物理层操作步骤。很多教程直接用qpskmod()一步到位掩盖了这些底层映射关系导致后续解调时相位模糊、误码率异常。2.1 比特分组与格雷码映射为什么不能按自然顺序排列假设原始比特流是[1 0 1 1 0 0 1 0]第一步必须将其两两分组[10, 11, 00, 10]。但关键来了这四个二进制组合对应星座点的顺序必须采用格雷码Gray Code映射而非00→0°、01→90°、10→180°、11→270°这样的自然顺序。为什么因为格雷码的定义是任意两个相邻码字之间仅有一位不同。在QPSK星座图中这意味着相邻星座点如I1,Q1 和 I1,Q-1之间只有一路Q路发生电平翻转当信道噪声导致接收点越过判决门限时大概率只错判为相邻点从而将双比特错误降为单比特错误显著降低整体BER。我们手动实现格雷码映射不调用bi2de()% 原始比特流示例 bits [1 0 1 1 0 0 1 0]; % 长度需为偶数 N_bits length(bits); N_symbols N_bits / 2; % 分组reshape为2行矩阵每列是一个2-bit符号 bit_groups reshape(bits, 2, N_symbols); % 2 x N_symbols % 格雷码映射表[b1 b0] - [I Q] % 格雷码顺序00-(1,1), 01-(1,-1), 11-(-1,-1), 10-(-1,1) % 对应标准QPSK星座第一象限(I0,Q0)为00顺时针旋转 gray_map [ ... 1, 1; ... % 00 - I1, Q1 1, -1; ... % 01 - I1, Q-1 -1, -1; ... % 11 - I-1, Q-1 -1, 1 ... % 10 - I-1, Q1 ]; % 将每组2-bit转换为十进制索引注意bit_groups(1,:)是高位bit_groups(2,:)是低位 dec_indices bit_groups(1,:) * 2 bit_groups(2,:); % 得到0,1,2,3 % 查表得到I/Q分量 I_component gray_map(dec_indices1, 1); % 1因MATLAB索引从1开始 Q_component gray_map(dec_indices1, 2); % 合成复数基带信号 s_baseband I_component 1i * Q_component; % 1 x N_symbols这段代码的关键在于dec_indices的计算方式bit_groups(1,:) * 2 bit_groups(2,:)。这里bit_groups(1,:)是高位MSBbit_groups(2,:)是低位LSB符合常规二进制权重。而格雷码映射表的行序严格按00、01、11、10排列确保了相邻星座点间仅一位差异。注意有些文献采用π/4-QPSK或Offset-QPSK其映射规则不同。本文聚焦标准QPSK即BPSK的正交复用所有后续步骤均以此映射为基础。若你看到别人代码里映射顺序不同先确认其采用的是哪种QPSK变体。2.2 成型滤波升余弦滤波器的滚降因子与抽样点选择调制后的符号是离散的矩形脉冲在频域会无限展宽严重干扰邻道。因此必须通过成型滤波器通常为升余弦滤波器RC Filter限制带宽。MATLAB没有内置rcosdesign函数没关系我们手写一个彻底搞清其参数意义。升余弦滤波器的时域脉冲响应为 $$ h(t) \frac{\sin(\pi t/T_s) \cos(\alpha \pi t/T_s)}{\pi t/T_s (1 - (2\alpha t/T_s)^2)} $$ 其中$T_s$ 是符号周期$\alpha$ 是滚降因子0 ≤ α ≤ 1。关键参数选择逻辑滚降因子 αα0时为理想奈奎斯特滤波器带宽1/(2Ts)但时域拖尾无穷长无法实现α1时带宽1/Ts时域收敛快但频谱利用率低。工程常用α0.35GSM或α0.25LTE。本文取α0.35。滤波器长度必须足够长以保证主瓣能量集中。经验公式长度 span * sps其中span是滤波器跨符号数通常取10sps是每符号采样点数决定后续DAC分辨率。% 设定参数 sps 8; % 每符号采样点数决定时间分辨率 span 10; % 滤波器跨符号数影响时域截断误差 alpha 0.35; % 滚降因子 L span * sps; % 滤波器总长度奇数便于中心对称 % 生成时间向量t范围为[-span/2 * Ts, span/2 * Ts] t (-L/2 : L/2-1) / sps; % 单位符号周期Ts % 计算升余弦脉冲响应h(t) % 处理分母为零的点t0 h zeros(L,1); for n 1:L if abs(t(n)) 1e-10 h(n) 1; % t0时极限值为1 else num sin(pi*t(n)) * cos(alpha*pi*t(n)); den pi*t(n) * (1 - (2*alpha*t(n))^2); h(n) num / den; end end % 归一化使滤波器增益为1直流增益1 h h / sum(h);这段代码生成了归一化的升余弦滤波器系数h。接下来将离散符号序列s_baseband进行插值补零并卷积% 符号上采样在每个符号后插入(sps-1)个零 s_upsampled upsample(s_baseband, sps); % 1 x (N_symbols * sps) % 成型滤波卷积 s_shaped filter(h, 1, s_upsampled); % 输出长度 len(s_upsampled) len(h) - 1 % 截断去除滤波器引起的前导和拖尾保留有效部分 % 有效数据起始位置滤波器群延迟 ≈ (L-1)/2 个采样点 delay floor((L-1)/2); s_tx s_shaped(delay1 : delay length(s_upsampled)); % 对齐输出这里upsample和filter的组合模拟了DAC将离散符号转换为连续波形的过程。delay的计算至关重要——它代表了滤波器引入的固定群延迟后续解调时必须补偿此延迟否则符号定时会偏移。2.3 载波调制正交上变频的相位基准统一成型后的基带信号s_tx是复数需调制到射频载波上。标准做法是乘以$e^{j2\pi f_c t}$再取实部得到实信号。但此处有一个极易被忽视的陷阱载波相位必须与调制时的星座图相位严格一致。例如你的格雷码映射中00对应(1,1)即相位45°。那么载波调制时cos(2πf_ct)和sin(2πf_ct)的初始相位必须为0才能保证I*cos - Q*sin的结果在t0时刻准确落在45°方向。如果载波有初始相位偏移φ整个星座图就会旋转φ角导致解调端判决失效。% 设定载波频率fc单位Hz和采样率fs单位Hz fc 1e6; % 1 MHz载波 fs sps / Ts; % fs sps * RsRs为符号率Ts1/Rs % 假设符号率Rs1e6 Baud则fs 8e6 Hz t_tx (0:length(s_tx)-1) / fs; % 时间向量 % 正交上变频I(t)*cos(2πfct) - Q(t)*sin(2πfct) % 注意s_tx是复数real(s_tx)为I分量imag(s_tx)为Q分量 I_tx real(s_tx); Q_tx imag(s_tx); s_rf I_tx .* cos(2*pi*fc*t_tx) - Q_tx .* sin(2*pi*fc*t_tx);这段代码实现了标准的正交调制。关键点在于- Q_tx .* sin(...)而非。这是由通信系统约定的复包络定义决定的复信号$s(t) I(t) jQ(t)$其对应的实信号为$Re{s(t)e^{j2\pi f_c t}} I(t)\cos(2\pi f_c t) - Q(t)\sin(2\pi f_c t)$。符号错了整个调制就反相了。2.4 信道建模AWGN的功率标定与带宽匹配将s_rf送入信道最常见的是加性高斯白噪声AWGN。但awgn()函数的snr参数指的是信号功率与噪声功率之比而这个“信号功率”是按什么基准算的很多代码直接写awgn(s_rf, snr, measured)看似省事实则埋雷。正确做法是先计算s_rf的有效功率再根据目标Eb/N0比特能量与噪声功率谱密度之比推导出实际SNR。因为Eb/N0才是通信系统性能的理论基准与调制方式、编码率强相关。QPSK中每个符号携带2比特故$E_b E_s / 2$其中$E_s$是符号能量。而$E_s P_s / R_s$$P_s$为信号平均功率$R_s$为符号率。噪声功率谱密度$N_0 P_n / B_n$$B_n$为噪声带宽。对于AWGN信道若接收滤波器带宽$B_n R_s(1\alpha)$则$P_n N_0 \times B_n$。因此目标SNR线性值为 $$ SNR \frac{P_s}{P_n} \frac{E_s R_s}{N_0 R_s (1\alpha)} \frac{E_s}{N_0 (1\alpha)} \frac{2 E_b}{N_0 (1\alpha)} $$% 设定目标Eb/N0dB Eb_N0_dB 10; % 示例值 alpha 0.35; % 计算理论SNR线性值 Eb_N0_lin 10^(Eb_N0_dB/10); SNR_lin 2 * Eb_N0_lin / (1 alpha); % 计算s_rf的实际功率均方值 Ps mean(abs(s_rf).^2); % 计算所需噪声功率 Pn Ps / SNR_lin; % 生成匹配带宽的高斯噪声 % 噪声标准差sigma sqrt(Pn) sigma sqrt(Pn); noise sigma * randn(size(s_rf)); % 加噪 s_rx s_rf noise;这段代码的核心是Ps mean(abs(s_rf).^2)——它精确计算了射频信号的平均功率而非依赖awgn的自动测量。sigma sqrt(Pn)确保了噪声功率严格匹配目标SNR。这样生成的S_rx其BER才能与理论曲线berawgn(Eb_N0_dB, qpsk)对标。3. QPSK解调链路从射频信号到比特流的逆向工程解调是调制的逆过程但绝非简单倒放。它包含载波恢复、符号定时恢复、匹配滤波、符号判决四大环节。任何一个环节失准都会导致误码率飙升。本节将逐层拆解重点揭示那些“看起来正常实则已失效”的隐蔽故障点。3.1 正交下变频本地振荡器的相位与频率同步接收端s_rx是实信号需先下变频到基带。理想情况下用与发射端同频同相的本地振荡器LO乘以s_rx再经低通滤波即可分离出I/Q分量。% 生成本地载波必须与发射端fc、相位完全一致 t_rx (0:length(s_rx)-1) / fs; LO_cos cos(2*pi*fc*t_rx); LO_sin sin(2*pi*fc*t_rx); % 正交下变频 I_down s_rx .* LO_cos; % 产生2fc分量和基带分量 Q_down s_rx .* LO_sin; % 低通滤波匹配滤波器即升余弦滤波器的镜像 % 这里复用前面设计的h但需注意h是基带成型滤波器其频响关于0对称 % 故可直接作为匹配滤波器使用 I_lp filter(h, 1, I_down); Q_lp filter(h, 1, Q_down); % 抽取每sps个采样点取一个得到符号速率采样 I_sampled I_lp(1:sps:end); Q_sampled Q_lp(1:sps:end);这段代码看似完美但隐藏着致命问题LO的相位φ_LO与发射端载波相位φ_TX不一致时I/Q分量会混合星座图旋转。例如若φ_LO φ_TX θ则解调后I I·cosθ Q·sinθQ -I·sinθ Q·cosθ整个星座图旋转θ角。解决方案不是“猜”一个θ去补偿而是在解调前加入载波相位恢复环路CPR。最简单有效的是M-th power estimatorM4 for QPSK% 对接收信号s_rx做4次方运算消除调制信息提取载波相位 s_rx_4 s_rx.^4; % 对s_rx_4做FFT找到峰值频率即4倍载波频率 N_fft 2^nextpow2(length(s_rx_4)); S4_fft fft(s_rx_4, N_fft); freq_vec (0:N_fft-1)*fs/N_fft; [~, idx_max] max(abs(S4_fft(1:floor(N_fft/2)))); f4_est freq_vec(idx_max); % 估计载波频率fc_est f4_est / 4 fc_est f4_est / 4; % 用fc_est重新生成LO频率同步 LO_cos_sync cos(2*pi*fc_est*t_rx); LO_sin_sync sin(2*pi*fc_est*t_rx); % 再次下变频...M-th power法利用QPSK信号的统计特性其4次方后调制信息消失只剩4倍频载波分量。通过FFT检测该分量即可估计出fc。这是无导频、盲估计的经典方法实测在SNR 10 dB时精度很高。3.2 匹配滤波与符号定时眼图张开度是唯一真理下变频后的I/Q信号经过匹配滤波即成型滤波器的共轭理论上应在每个符号周期的中心点获得最大信噪比。但实际中由于采样时钟偏差采样点会漂移导致眼图闭合。% 匹配滤波与成型滤波器h相同 I_matched filter(h, 1, I_down); Q_matched filter(h, 1, Q_down); % 绘制眼图关键诊断工具 N_eye 100; % 显示100个符号的眼图 I_eye reshape(I_matched(1:N_eye*sps), sps, N_eye); Q_eye reshape(Q_matched(1:N_eye*sps), sps, N_eye); figure; plot(I_eye, Q_eye, b., MarkerSize, 1); xlabel(I); ylabel(Q); title(QPSK Eye Diagram (I-Q Plane)); grid on;眼图不是为了好看而是为了定位最佳采样时刻。观察眼图在I-Q平面上的“眼睛”张开最大的水平线即Q0附近其横坐标就是I分量的最佳采样点同理观察I0附近的竖直线其纵坐标就是Q分量的最佳采样点。这个点通常不在sps/2处而是在sps/2 delta处delta由时钟偏差决定。因此不能简单地I_sampled I_matched(1:sps:end)而应% 手动寻找最佳采样偏移delta示例扫描-2到2个采样点 best_delta 0; best_opening 0; for delta -2:2 I_test I_matched(delta1:sps:end); Q_test Q_matched(delta1:sps:end); % 计算眼图张开度统计I分量在决策门限0附近的穿越点密度 crossing_density sum(abs(diff(sign(I_test))) 0) / length(I_test); if crossing_density best_opening best_opening crossing_density; best_delta delta; end end % 使用最佳偏移采样 I_opt I_matched(best_delta1:sps:end); Q_opt Q_matched(best_delta1:sps:end);这个best_delta搜索过程模拟了实际接收机中的定时恢复算法如Gardner算法。它证明了一个事实没有完美的采样点只有通过眼图反馈不断优化的采样点。3.3 符号判决与比特映射格雷码逆映射的鲁棒性设计经过匹配滤波和定时采样得到I_opt和Q_opt下一步是符号判决。最简单的是硬判决I 0 ? 0 : 1Q 0 ? 0 : 1。但这在噪声大时极易出错。更鲁棒的做法是软判决即计算每个接收符号到四个星座点的欧氏距离选择最近点% 四个星座点格雷码顺序 constellation [11i, 1-1i, -1-1i, -11i]; % 初始化判决符号数组 symbols_decoded zeros(1, length(I_opt)); for k 1:length(I_opt) rx_symbol I_opt(k) 1i * Q_opt(k); % 计算到各星座点的距离平方 dist_sq abs(rx_symbol - constellation).^2; [~, idx_min] min(dist_sq); symbols_decoded(k) idx_min; % 1,2,3,4对应00,01,11,10 end % 格雷码逆映射将索引转回2-bit % 映射表1-00, 2-01, 3-11, 4-10 bits_decoded zeros(1, 2*length(symbols_decoded)); for k 1:length(symbols_decoded) switch symbols_decoded(k) case 1 bits_decoded(2*k-1:2*k) [0 0]; case 2 bits_decoded(2*k-1:2*k) [0 1]; case 3 bits_decoded(2*k-1:2*k) [1 1]; case 4 bits_decoded(2*k-1:2*k) [1 0]; end end这段代码实现了完整的软判决流程。dist_sq计算的是距离平方避免开方运算提升效率。switch语句确保了逆映射与前面的格雷码映射严格对应。最终bits_decoded就是解调出的比特流。注意实际系统中软判决输出的是LLRLog-Likelihood Ratio用于后续信道译码。本文为简化仅展示硬判决到比特的映射。但即使硬判决也建议用距离最小法而非简单的符号比较因为它对幅度失衡如I/Q增益不平衡更具鲁棒性。3.4 BER计算与理论对标如何证明你的仿真“可信”最后一步计算误比特率BER并与理论值对比。这是检验整个仿真链路是否正确的终极标尺。% 计算误比特数 num_errors sum(xor(bits(1:length(bits_decoded)), bits_decoded)); total_bits length(bits_decoded); BER_sim num_errors / total_bits; % 理论BERQPSKAWGN BER_theory qfunc(sqrt(2 * 10^(Eb_N0_dB/10))); % qfunc(x) 0.5*erfc(x/sqrt(2)) fprintf(Simulated BER: %.2e\n, BER_sim); fprintf(Theoretical BER: %.2e\n, BER_theory); fprintf(Relative error: %.2f%%\n, abs(BER_sim - BER_theory)/BER_theory * 100);qfunc是MATLAB内置函数计算标准正态分布的右尾概率。QPSK的理论BER公式为$Q(\sqrt{2E_b/N_0})$这是香农信息论推导出的闭式解。关键验证点当Eb_N0_dB 10时理论BER ≈ 3.87e-6。如果你的仿真BER在1e-5量级说明链路基本正确若在1e-3量级那一定在某个环节出了问题——此时不要盲目调参而应回溯检查眼图是否张开星座图是否聚集在四个点噪声功率是否标定准确格雷码映射是否一致我曾遇到一个案例BER始终比理论高10倍。排查三天最终发现是成型滤波器h未归一化导致信号功率被放大awgn加的噪声相对不足SNR虚高。归一化h后BER立刻回归理论值。这印证了那句话仿真不是调参游戏而是物理定律的数值验证。4. 全链路MATLAB代码整合与调试技巧现在将前述所有模块整合为一个可运行、可调试、可扩展的完整脚本。代码结构清晰每段功能独立便于你逐模块验证。%% QPSK End-to-End Simulation (No Toolbox Required) % Author: Senior Comm Engineer % Date: 2024 % Description: Full QPSK modulation/demodulation chain with physical layer details. %% 0. Clear Workspace Set Parameters clear; clc; close all; rng(42); % Fixed seed for reproducibility % System parameters Rs 1e6; % Symbol rate (Baud) Ts 1/Rs; % Symbol period (s) sps 8; % Samples per symbol alpha 0.35; % Roll-off factor fc 1e6; % Carrier frequency (Hz) fs sps * Rs; % Sampling frequency (Hz) Eb_N0_dB 10; % Target Eb/N0 (dB) % Generate random bit stream (length must be even) N_bits 10000; bits randi([0,1], 1, N_bits); %% 1. Modulation Chain % 1.1 Bit grouping and Gray mapping N_symbols N_bits / 2; bit_groups reshape(bits, 2, N_symbols); dec_indices bit_groups(1,:) * 2 bit_groups(2,:); gray_map [1,1; 1,-1; -1,-1; -1,1]; I_component gray_map(dec_indices1, 1); Q_component gray_map(dec_indices1, 2); s_baseband I_component 1i * Q_component; % 1.2 Pulse shaping (Root Raised Cosine) span 10; L span * sps; t (-L/2 : L/2-1) / sps; h zeros(L,1); for n 1:L if abs(t(n)) 1e-10 h(n) 1; else num sin(pi*t(n)) * cos(alpha*pi*t(n)); den pi*t(n) * (1 - (2*alpha*t(n))^2); h(n) num / den; end end h h / sum(h); % Normalize % Upsample and filter s_upsampled upsample(s_baseband, sps); s_shaped filter(h, 1, s_upsampled); delay floor((L-1)/2); s_tx s_shaped(delay1 : delay length(s_upsampled)); % 1.3 Carrier modulation t_tx (0:length(s_tx)-1) / fs; I_tx real(s_tx); Q_tx imag(s_tx); s_rf I_tx .* cos(2*pi*fc*t_tx) - Q_tx .* sin(2*pi*fc*t_tx); % 1.4 AWGN channel (Eb/N0 based) Ps mean(abs(s_rf).^2); Eb_N0_lin 10^(Eb_N0_dB/10); SNR_lin 2 * Eb_N0_lin / (1 alpha); Pn Ps / SNR_lin; sigma sqrt(Pn); noise sigma * randn(size(s_rf)); s_rx s_rf noise; %% 2. Demodulation Chain % 2.1 Carrier recovery (4th power method) s_rx_4 s_rx.^4; N_fft 2^nextpow2(length(s_rx_4)); S4_fft fft(s_rx_4, N_fft); freq_vec (0:N_fft-1)*fs/N_fft; [~, idx_max] max(abs(S4_fft(1:floor(N_fft/2)))); f4_est freq_vec(idx_max); fc_est f4_est / 4; % 2.2 Down-conversion with estimated fc t_rx (0:length(s_rx)-1) / fs; LO_cos cos(2*pi*fc_est*t_rx); LO_sin sin(2*pi*fc_est*t_rx); I_down s_rx .* LO_cos; Q_down s_rx .* LO_sin; % 2.3 Matched filtering I_matched filter(h, 1, I_down); Q_matched filter(h, 1, Q_down); % 2.4 Timing recovery (eye diagram based) N_eye 100; I_eye reshape(I_matched(1:N_eye*sps), sps, N_eye); Q_eye reshape(Q_matched(1:N_eye*sps), sps, N_eye); % Find optimal sampling offset best_delta 0; best_opening 0; for delta -2:2 I_test I_matched(delta1:sps:end); Q_test Q_matched(delta1:sps:end); crossing_density sum(abs(diff(sign(I_test))) 0) / length(I_test); if crossing_density best_opening best_opening crossing_density; best_delta delta; end end % Sample at optimal point I_opt I_matched(best_delta1:sps:end); Q_opt Q_matched(best_delta1:sps:end); % 2.5 Symbol decision and bit mapping constellation [11i, 1-1i, -1-1i, -11i]; symbols_decoded zeros(1, length(I_opt)); for k 1:length(I_opt) rx_symbol I_opt(k) 1i * Q_opt(k); dist_sq abs(rx_symbol - constellation).^2; [~, idx_min] min(dist_sq); symbols_decoded(k) idx_min; end bits_decoded zeros(1, 2*length(symbols_decoded)); for k 1:length(symbols_decoded) switch symbols_decoded(k) case 1, bits_decoded(2*k-1:2*k) [0 0]; case 2, bits_decoded(2*k-1:2*k) [0 1]; case 3, bits_decoded(2*k-1:2*k) [1 1]; case 4, bits_decoded(2*k-1:2*k) [1 0]; end end %% 3. Performance Evaluation num_errors sum(xor(bits(1:length(bits_decoded)), bits_decoded)); total_bits length(bits_decoded); BER_sim num_errors / total_bits; BER_theory qfunc(sqrt(2 * 10^(Eb_N0_dB/10))); fprintf(\n QPSK Simulation Results \n); fprintf(Target Eb/N0: %.1f dB\n, Eb_N0_dB); fprintf(Simulated BER: %.2e\n, BER_sim); fprintf(Theoretical BER: %.2e\n, BER_theory); fprintf(Relative Error: %.2f%%\n, abs(BER_sim - BER_theory)/BER_theory * 100); %% 4. Visualization (Optional but Highly Recommended) figure(Name, QPSK Constellation); scatter(real(s_baseband), imag(s_baseband), filled); hold on; scatter(I_opt, Q_opt, r*, MarkerSize, 8); xlabel(I); ylabel(Q); title(QPSK Constellation (Red: Received Symbols)); grid on; figure(Name, QPSK Eye Diagram); plot(I_eye, Q_eye, b., MarkerSize, 1); xlabel(I); ylabel(Q); title(QPSK Eye Diagram (I-Q Plane)); grid on;4.1 调试技巧如何快速定位链路故障这套代码不是“一次成功”的魔法而是你调试的起点。以下是我在十年通信仿真中总结的黄金调试清单分段注入已知信号将s_baseband直接赋值为[11i, 1-1i, -1-1i, -11i]四个确定符号跳过调制直接进入解调。若解调后能正确还原[00,01,11,10]说明解调链路无问题否则问题在解调端。关闭噪声验证确定性路径注释掉加噪行让s_rx s_rf。此时BER应为0。若不为0说明调制/解调存在系统性偏差如映射不一致本文还有配套的精品资源点击获取
![[光学原理与应用-521]：对光的错误理解与纠偏](http://pic.xiahunao.cn/yaotu/[光学原理与应用-521]：对光的错误理解与纠偏)
