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星座图中，这意味着：

• 相邻星座点（如I=1,Q=1 和 I=1,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星座：第一象限(I>0,Q>0)为00，顺时针旋转 gray_map = [ ... 1, 1; ... % 00 -> I=1, Q=1 1, -1; ... % 01 -> I=1, Q=-1 -1, -1; ... % 11 -> I=-1, Q=-1 -1, 1 ... % 10 -> I=-1, Q=1 ]; % 将每组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_indices+1, 1); % +1因MATLAB索引从1开始 Q_component = gray_map(dec_indices+1, 2); % 合成复数基带信号 s_baseband = I_component + 1i * Q_component; % 1 x N_symbols

这段代码的关键在于dec_indices的计算方式：bit_groups(1,:) * 2 + bit_groups(2,:)。这里bit_groups(1,:)是高位（MSB），bit_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.35（GSM）或α=0.25（LTE）。本文取α=0.35。
• 滤波器长度：必须足够长以保证主瓣能量集中。经验公式：长度 =span * sps，其中span是滤波器跨符号数（通常取10），sps是每符号采样点数（决定后续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) % 处理分母为零的点（t=0） h = zeros(L,1); for n = 1:L if abs(t(n)) < 1e-10 h(n) = 1; % t=0时，极限值为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(delay+1 : 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的结果在t=0时刻准确落在45°方向。如果载波有初始相位偏移φ，整个星座图就会旋转φ角，导致解调端判决失效。

% 设定载波频率fc（单位：Hz）和采样率fs（单位：Hz） fc = 1e6; % 1 MHz载波 fs = sps / Ts; % fs = sps * Rs，Rs为符号率，Ts=1/Rs % 假设符号率Rs=1e6 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/N0（dB） 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 estimator（M=4 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平面上的“眼睛”张开最大的水平线（即Q=0附近），其横坐标就是I分量的最佳采样点；同理，观察I=0附近的竖直线，其纵坐标就是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(delta+1:sps:end); Q_test = Q_matched(delta+1: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_delta+1:sps:end); Q_opt = Q_matched(best_delta+1:sps:end);

这个best_delta搜索过程，模拟了实际接收机中的定时恢复算法（如Gardner算法）。它证明了一个事实：没有完美的采样点，只有通过眼图反馈不断优化的采样点。

3.3 符号判决与比特映射：格雷码逆映射的鲁棒性设计

经过匹配滤波和定时采样，得到I_opt和Q_opt，下一步是符号判决。最简单的是硬判决：I > 0 ? 0 : 1，Q > 0 ? 0 : 1。但这在噪声大时极易出错。

更鲁棒的做法是软判决，即计算每个接收符号到四个星座点的欧氏距离，选择最近点：

% 四个星座点（格雷码顺序） constellation = [1+1i, 1-1i, -1-1i, -1+1i]; % 初始化判决符号数组 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就是解调出的比特流。

注意：实际系统中，软判决输出的是LLR（Log-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; % 理论BER（QPSK，AWGN） 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_indices+1, 1); Q_component = gray_map(dec_indices+1, 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(delay+1 : 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(delta+1:sps:end); Q_test = Q_matched(delta+1: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_delta+1:sps:end); Q_opt = Q_matched(best_delta+1:sps:end); % 2.5 Symbol decision and bit mapping constellation = [1+1i, 1-1i, -1-1i, -1+1i]; 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 调试技巧：如何快速定位链路故障

这套代码不是“一次成功”的魔法，而是你调试的起点。以下是我在十年通信仿真中总结的黄金调试清单：

1. 分段注入已知信号：

   • 将s_baseband直接赋值为[1+1i, 1-1i, -1-1i, -1+1i]（四个确定符号），跳过调制，直接进入解调。若解调后能正确还原[00,01,11,10]，说明解调链路无问题；否则，问题在解调端。

2. 关闭噪声，验证确定性路径：

   • 注释掉加噪行，让s_rx = s_rf。此时BER应为0。若不为0，说明调制/解调存在系统性偏差（如映射不一致

本文还有配套的精品资源，点击获取