简介基于Viterbi-Viterbi的载波相位估计算法MATLAB实现面向光通信、数字信号处理领域的研究者与工程师解决高精度载波相位恢复问题适用于CAP调制等相干通信场景。压缩包共1个文件为CPE.m脚本大小899B代码虽短但结构完整包含观测模型定义、距离计算、基于动态规划的路径跟踪、回溯及相位补偿等核心环节。该资源目前已有2143人学习可见此算法在相位估计方向的应用热度。通过阅读脚本可直观理解Viterbi-Viterbi如何从离散码字估计扩展到连续相位空间掌握最大似然序列估计在相位恢复中的实际实现思路为优化解调性能提供可运行的参考代码。1. 这是载波同步里的老牌算法先搞懂它在解决什么做过突发通信解调的人都有过这种经历定时同步做完了频偏估得只剩几十赫兹残差但星座图就是整体转了一个角度怎么判决都错。这个角度就是载波初始相位它既不知道大小又没法用环路慢慢追——突发就那么短几百个符号必须一次估计出来。Viterbi-Viterbi 算法就是干这个的它把接收到的 M-PSK 信号连续做 M 次方让所有符号的调制相位全部折叠到同一个点之后剩下的就只有载波相位和噪声再做块内平均把噪声压下去最后取相位除以 M就得到了载波相位估计。整个过程是盲估计不需要已知符号计算量也小QPSK 场景下几个几十行的函数就能跑通适合做突发解调、卫星通信接收机和相干光通信 DSP 的工程师也适合正在搭载波同步链路的通信方向学生。这篇文章就把它从原理到代码到参数踩坑完整拆开。2. 算法原理M 次方非线性把调制信息剥掉2.1 M-PSK 信号为什么让相位估计这件事变难先看接收信号模型。发送端发出的是 M-PSK 符号 (a_k e^{j(2\pi m/M)})(m) 取 (0,1,\dots,M-1)信息就藏在离散的相位里。经过信道后接收信号可以写成[ r_k a_k \cdot e^{j(2\pi \Delta f k T \theta)} n_k ]其中 (\theta) 就是我们要估计的载波初始相位(\Delta f) 是残余频偏(n_k) 是复高斯白噪声。如果天线端下行转换时本地振荡器和发端不同步(\theta) 可以是任意值从 (-180^\circ) 到 (180^\circ) 都有可能。这时候问题就来了如果直接对 (r_k) 取平均再取相位(\theta) 是会被淹没在调制信息里的。举个例子QPSK 的四个符号相位分别是 (45^\circ, 135^\circ, 225^\circ, 315^\circ)它们均匀分布在单位圆上平均之后趋近于零载波相位的痕迹完全消失。哪怕用最大似然准则去估计在未知符号序列的情况下也需要对 (M^N) 种符号组合做联合搜索N 一长就完全不可行。所以核心矛盾是想从接收信号里提取载波相位但调制相位本身是个离散随机量它和载波相位叠加在一起必须先把调制分量去掉剩下的才是干净的单频信号这样再做平均才有意义。Viterbi-Viterbi 的巧妙之处就在于用了一个固定的非线性运算就把调制分量彻底消除复杂度只有 O(N)。2.2 去调制M 次方运算的数学本质抓住一个关键性质任何一个 M-PSK 符号它的 M 次方恒等于 1。[ a_k^M \left(e^{j(2\pi m/M)}\right)^M e^{j2\pi m} 1 ]不管 (m) 取什么值乘以 M 之后都落在复平面的实轴 1 上。这个性质就是整个算法的地基。把接收信号 (r_k) 直接做 M 次方展开忽略噪声只看信号项[ r_k^M \left(a_k e^{j(2\pi \Delta f kT \theta)}\right)^M a_k^M e^{j M (2\pi \Delta f kT \theta)} e^{j M \theta} \cdot e^{j M \cdot 2\pi \Delta f kT} ]调制项 (a_k^M) 消失了剩下的是两项的乘积第一项是常数 (e^{jM\theta})与符号索引 (k) 无关第二项是残余频偏引起的相位积累与 (k) 成正比。如果残余频偏已经被压得很小第二项在短块内近似为常数那么 (r_k^M) 本质上就是一个被复噪声污染了的恒定复数。对这个复数做块内平均[ \frac{1}{N}\sum_{k1}^{N} r_k^M \approx e^{jM\theta} \frac{1}{N}\sum n_k ]噪声项被平均掉约 (1/\sqrt{N}) 的幅度所以平均结果的相位就很接近 (M\theta)。取相位再除以 M[ \hat{\theta} \frac{1}{M} \cdot \arg\left(\sum_{k1}^{N} r_k^M\right) ]于是相位估计完成。整个流程就是三步M 次方去调制、块内平均压噪声、取相位还原。2.3 M 次方会放大噪声这个算法不是万能的展开 (r_k^M) 时噪声项可不是老老实实的。二项式展开写出来[ (a_k e^{j\varphi} n_k)^M a_k^M e^{jM\varphi} M a_k^{M-1} e^{j(M-1)\varphi} n_k \binom{M}{2} a_k^{M-2} e^{j(M-2)\varphi} n_k^2 \cdots n_k^M ]第一项是信号但后面每一项都有噪声而且噪声的功率被放大了。最关键的是线性项 (M a_k^{M-1} n_k)它的功率是原噪声功率的 (M^2) 倍。这意味着 M 次方之后信号的等效信噪比显著下降信号项功率不变噪声功率却膨胀了 (M^2) 量级。M 越大去调制后信噪比损失越惨。这也是为什么 QPSKM4用这个算法效果不错而 16PSK、32PSK 用这个算法几乎没法工作——高阶调制下噪声放大太狠块内平均根本压不回来。还有一个被放大的因素是残余频偏。从上面的推导能看到残余频偏在 M 次方后变成了原来的 M 倍。本来残余频偏小到几十赫兹M 次方后变成几百赫兹块长一长块内的相位转一圈平均结果直接趋近于零估计值变成随机数。这个特性直接影响下一节要讲的块长选择。所以这个算法的适用边界很清楚M-PSK 调制通常是 M2 或 M4残余频偏已经做过粗校准信噪比在高斯信道下大约高于 0~5dB 的场景。超出这个范围要么结果方差爆炸要么相位模糊误判都需要其他方法来兜底。3. 用 Python 从零实现 Viterbi-Viterbi 相位估计3.1 生成带相位偏移的 QPSK 测试信号先把测试平台搭起来。用 NumPy 生成一段 QPSK 基带符号加载波相位偏移和复高斯噪声这样后续每个环节都能对照真实相位验证。import numpy as np def generate_qpsk_symbols(n_symbols, seed42): 生成 QPSK 符号序列星座点旋转 45 度能量归一化 rng np.random.default_rng(seed) idx rng.integers(0, 4, n_symbols) # 星座点在 45/135/225/315 度避免落到实轴上 symbols np.exp(1j * (np.pi / 4 idx * np.pi / 2)) return symbols def add_phase_noise(symbols, theta_deg, snr_db, seed42): 叠加固定相位偏移和复高斯白噪声 rng np.random.default_rng(seed) theta np.deg2rad(theta_deg) snr_linear 10 ** (snr_db / 10) noise_power 1.0 / snr_linear noise np.sqrt(noise_power / 2) * ( rng.standard_normal(len(symbols)) 1j * rng.standard_normal(len(symbols)) ) rx symbols * np.exp(1j * theta) noise return rx逻辑说明generate_qpsk_symbols里把星座点旋转了 45 度这样四个符号都不在实轴或虚轴上后期画星座图更容易观察旋转量。add_phase_noise里的noise_power 1.0 / snr_linear是因为符号能量归一到 1噪声功率就是信噪比的倒数。参数说明theta_deg是真实相位偏移建议测试时分别用 10 度、80 度、150 度覆盖算法能估计和无模糊的不同区间snr_db建议在 0 到 20 之间调整观察低信噪比时的退化。3.2 核心估计函数非线性变换与块内平均核心代码只有十几行但拆开讲每一步的意图。def vv_phase_estimate(rx, M, block_len): Viterbi-Viterbi 盲载波相位估计 参数: rx: 接收复数基带序列 M: PSK 调制阶数QPSK 取 4 block_len: 每个估计块的符号数必须大于 8推荐 32~256 返回: est_phase: 每个块的模糊相位估计弧度范围 [-pi/M, pi/M) n_blocks len(rx) // block_len est_phase np.zeros(n_blocks) for i in range(n_blocks): block rx[i * block_len : (i 1) * block_len] y block ** M # 1. M 次方去调制 avg np.mean(y) # 2. 块内平均压制噪声 est_phase[i] np.angle(avg) / M # 3. 取相位并还原 return est_phase逻辑说明block ** M在 Python 里是逐元素的复数幂运算等价于对每个符号的幅度取 M 次方、相位乘 M。对 M-PSK 信号这一步把所有星座点压到同一个相位上。np.mean(y)对块内所有去调制后的复数取平均噪声的实部和虚部相互抵消信号项因为是常数所以保留下来。np.angle(avg)返回平均值的相位范围是 ([-\pi, \pi))除以 M 之后得到的就是模糊的相位估计。参数说明block_len是最关键的自由参数。它太小噪声平均不充分估计方差大它太大残余频偏在块内积累的相位旋转会让平均结果衰减。典型范围是 32 到 256突发场景下直接用整个突发做一次估计也行。3.3 相位模糊消除用独特字把 90 度歧义掰回来上一节代码输出的相位范围是 ([-\pi/M, \pi/M))对 QPSK 也就是 (-45^\circ) 到 (45^\circ) 之间。但真实相位可以是任意值。假设真实相位是 (80^\circ)(M\theta 320^\circ)取相位后得到的是 (-40^\circ)除以 4 得到 (-10^\circ)和真实值差 90 度整数倍。这就是相位模糊。消除相位模糊的常见做法是在突发帧里放一段已知符号作为独特字UW然后用候选相位旋转接收信号和这段已知符号做相关相关峰最大的候选就是真实相位。def resolve_ambiguity(rx_uw, M, est_phase, uw_symbols): 用独特字相关性消除 M 重相位模糊 参数: rx_uw: 接收端提取的独特字符号段含相位偏移和噪声 M: 调制阶数 est_phase: VV 估计的模糊相位标量 uw_symbols: 已知独特字复数符号序列 返回: phase_corrected: 消除模糊后的相位估计弧度 candidates est_phase 2 * np.pi * np.arange(M) / M best_phase candidates[0] best_score -np.inf for c in candidates: rotated rx_uw * np.exp(-1j * c) # 内积的绝对值作为相关度量 score np.abs(np.dot(rotated, uw_symbols.conj())) / len(uw_symbols) if score best_score: best_score score best_phase c return best_phase逻辑说明候选相位是模糊估计值加上 (0, 2\pi/M, 4\pi/M, \dots) 共 M 个。对每个候选相位旋转接收的 UW 符号如果旋转量正确信号在原始星座点上对齐内积的模会显著大于其他候选反之信号被转到错误位置内积模很小。最后选出相关峰最大的候选相位作为最终结果。参数说明uw_symbols的长度直接影响模糊判决成功率经验上至少 8 到 16 个符号少于 8 个在低信噪比下很容易选错候选。这段 UW 是通信帧结构的一部分会占用有效载荷所以实际系统要在同步可靠性和开销之间取舍。4. 三个必调参数块长度、M 值、实现结构怎么定4.1 块长度 N估计精度与容忍频偏的折衷块长度 N 是 Viterbi-Viterbi 里影响最大、也最需要根据场景反复试的参数。先看噪声对相位估计的影响。块内平均的作用是把噪声通过统计平均压低平均后噪声功率近似降为原来的 (1/N)因此估计相位误差的标准差大致与 (1/\sqrt{N}) 成正比。把块长从 16 提到 64估计误差理论上减少一半。所以从纯噪声角度看N 越大越好。但残余频偏完全改变了这个单调关系。上一章推导过M 次方之后残余频偏被放大 M 倍每个符号的残留相位旋转是 (M \cdot 2\pi \Delta f T)。如果块长 N 大到让整个块内累计旋转明显超过一个弧度那么块首和块尾的复数向量方向差别很大平均结果互相抵消平均幅值急剧下降等效信噪比反而变差。以 10 Mbaud 符号率、残余频偏 100 Hz、QPSK 为例M4 时等效频偏是 400 Hz一个符号的旋转是 (2\pi \times 400 / 10^7 \approx 0.00025) 弧度。块长 256 时累计旋转约 0.065 弧度影响很小块长 4096 时累计旋转约 1 弧度平均幅值已经明显衰减。所以这个场景下 256 很安全4096 就不行了。实测调参建议先在高信噪比比如 15dB下固定频偏为零逐渐加大 N 直到估计误差不再明显减小记下 N1然后固定 N 逐渐加大频偏观察估计值开始跳变的频偏值记下这个边界最后取一个同时满足精度和频偏容忍的折衷。一般来说 32 到 256 是绝大多数 PSK 场景的安全区间。4.2 调制阶数 M设错了全盘皆输M 值看起来是算法本身的固有参数但实际工程中经常有人在这里翻车。M 必须精确等于 PSK 调制阶数。QPSK 用 M4BPSK 用 M28PSK 用 M8。如果把 M 设成了实际调制阶数的倍数比如 QPSK 信号用 M8会发生什么QPSK 符号的四次方是 1但八次方就不是了。设符号相位是 (45^\circ)八次方相位是 (360^\circ)OK但符号相位是 (135^\circ) 时八次方相位就是 (1080^\circ)取相位后是 (0^\circ) 和 (0^\circ)仔细算一下(135^\circ \times 8 1080^\circ)减去 (2 \times 360^\circ720^\circ)剩 (360^\circ)正好还是 1。等等QPSK 四个符号的八次方也都等于 1因为 ((e^{j45^\circ})^8 e^{j360^\circ}1)((e^{j135^\circ})^8 e^{j1080^\circ}e^{j360^\circ}1)。看起来 M8 对 QPSK 也能去调制。那如果把 M 设小了QPSK 用 M2平方之后星座变成两个点在 0 和 180 度交替调制信息并没有被完全去除平均结果在两个方向上互相拉扯相位估计会在两个值之间跳。结论M 是调制阶数的整数倍时去调制依然有效但相位估计范围缩小了。M8 的估计范围变成 ([-22.5^\circ, 22.5^\circ])需要消除的模糊数从 4 重变成 8 重独特字判决的出错概率也变大。M4 是 QPSK 的最佳选择范围合适、模糊数少、噪声放大最小。BPSK 则要压到 M2。4.3 滑动平均 vs 块处理流式场景怎么选前面代码实现的是块处理把一段符号切成不重叠的块每块出一个相位估计值。这在突发通信里很自然来一个 burst 估一次burst 内部共享同一个相位。但连续传输场景下逐符号的相位变化也需要跟踪这个时候要换成滑动平均结构。滑动平均的代码逻辑是把前面的核心函数应用在一个滑窗上def vv_sliding_phase_estimate(rx, M, window_len): 滑窗式 Viterbi-Viterbi输出每个符号位置的相位估计 n len(rx) est np.zeros(n, dtypefloat) half window_len // 2 for k in range(n): left max(0, k - half) right min(n, k half) block rx[left:right] y block ** M avg np.mean(y) est[k] np.angle(avg) / M return est逻辑说明对每个符号取其前后各half个符号构成的窗口做平均。这样每个符号都用自己的邻域信息估计相位适用于频偏随时间缓慢漂移的连续流场景。窗口越大估计越平滑但跟踪频偏变化的能力越差。注意这个实现的性能问题每个符号都做一次 M 次方和平均复杂度是 O(N·window_len)。实际工程中一般不会逐符号滑动而是用半重叠块处理——块长 128、步长 64每两个相邻块之间共享 64 个符号然后对两个估计值做线性插值。这既保留了跟踪能力又把运算量压缩到块处理的两倍左右。5. 避坑指南Viterbi-Viterbi 的五个翻车现场5.1 相位模糊估计结果永远差 90 度的整数倍现象星座图收敛得非常好四个点清晰分明但解调误码率始终在 30% 到 50% 之间徘徊怎么调都压不下去。把星座图旋转 90 度再和发射星座对齐发现完全吻合。原因这是相位模糊最典型的表现。(M) 次方取相位后估计值被压缩进 ([-\pi/M, \pi/M))真实相位超过这个范围的部分被折叠。QPSK 的真实相位 (80^\circ) 在估计时变成了 (-10^\circ)两者相差 90 度恰好是整数倍模糊。解决用独特字做相关判决或者改用差分 QPSK 编码。差分编码下 90 度旋转不会影响相邻符号的相位差解调结果不受模糊影响代价是误码率稍有上升。独特字方案需要额外开销但对非差分系统是必须的。5.2 残余频偏引起的块内相位旋转现象把残余频偏从小到大慢慢加一开始估计值稳定加到某个值后估计结果突然变成随机跳变而且每块之间毫无关联。回头看块内平均值的幅度发现小得可怜。原因残余频偏在 M 次方后被放大了 M 倍块内累计相位旋转过大块首和块尾的向量方向相反平均后互相抵消。幅度趋近于零时相位由噪声主导估计值就是随机数。解决在做 Viterbi-Viterbi 之前先用频偏粗估计算法把残余频偏压下去。常见做法是用接收信号的 M 次方做 FFT 测频偏或者用延迟相关法把频偏压到符号率的百分之一以内。另一个应急手段是缩短块长但要注意这会同时抬高噪声方差属于两难取舍。5.3 低信噪比下的估计方差爆炸现象SNR 在 10dB 以上时估计误差很小把 SNR 降到 3dB误码率没有渐进恶化而是突然跳到一个很高的水平。观察估计值和真实值的差不是均匀分布的而是集中在几个固定偏差上。原因M 次方运算对噪声项有放大效应低信噪比下 M 次方后的等效信噪比进一步恶化。当这个等效 SNR 跌到某个门限以下平均值的相位开始被少数大噪声样本支配出现类似判决错误的跳变。这就是非线性估计的门限效应。解决增加块长是最直接的手段但要受频偏约束。另一个思路是放宽工作门限允许第一级 Viterbi-Viterbi 输出一个不那么精确的相位旋转信号后用判决引导做第二级细化。两级级联可以把工作 SNR 往下压几个 dB。如果预算允许也可以在帧结构里加入短训练序列做数据辅助估计。5.4 M 值设错导致星座不收敛现象明明发送的是 QPSK星座图却出现了 8 个点或者 4 个点但每个点周围都有两个小簇。估计值和真实相位的关系完全对不上。原因M 和实际调制阶数不匹配。M 设成 8 对 QPSK 依然能去调制但相位估计范围缩小了一半模糊数翻倍独特字判决更容易出错M 设成 2 则去调制不彻底星座点还保留着调制相位的痕迹平均结果在两个方向之间摆动。解决在工程实现里调制阶数应该在解调器的参数配置里写死并和发射端约定一致。如果系统支持多种调制方式自适应切换要让 VV 的 M 值跟随调制方式配置同步变化而不是固定一个值。调试时先打印星座图确认收敛点数等于调制阶数再往下走。5.5 独特字长度不足导致模糊判决错误现象高信噪比时一切正常SNR 降到 8dB 以下误码率出现周期性突增。统计发现每几百个突发就有几个突发的星座整体旋转了 90 度。原因独特字长度太短相关峰不够尖锐。候选相位之间的相关性只有微弱差异噪声一上来就容易选到错误的候选。8 个符号的独特字在 10dB 下选错概率已经不可忽略。解决独特字长度至少取 16 个符号条件允许时取 32。加长独特字的收益在低信噪比下非常显著因为相关判决的错误概率随长度指数下降。如果帧长度受限可以改为用导频符号在多个段之间做联合判决但这会显著增加实现复杂度。6. 进阶玩法级联结构与 MSE 验证曲线6.1 两级级联粗估计加判决引导细化Viterbi-Viterbi 用块内平均压噪声但块内平均终究是个线性操作简单平均的残差在高阶调制下仍然偏大。常见做法是在它后面再接一级判决引导估计第一级用 VV 估出粗相位对信号做旋转旋转后星座点已经大致落在正确位置对每个符号做硬判决把判决结果当作已知符号然后用最大似然估计公式做第二级细化。需要注意第二级细化只能在第一级误差小于 ±π/M 时才有效。如果第一级已经发生相位模糊错误第二级会沿着错误方向继续优化。所以第二级输入要加一个保护计算判决后的残差模值如果平均残差超过阈值就判定第一级失败丢弃这次估计而不是输出错上加错的结果。6.2 验证方法用 MSE 随 SNR 的曲线判断实现是否正确拿到一个 Viterbi-Viterbi 实现最快验证它是否正确的方法是跑蒙特卡洛仿真画 MSE 对 SNR 的曲线。在固定 SNR 下生成几千次独立的相位偏移测试每次记录估计值和真实值的相位差把所有差值取平方再平均得到该 SNR 下的 MSE。把 SNR 从 0dB 扫到 20dB画出的 MSE 曲线应该呈现一条平滑下降的直线对数坐标下没有突变。我自己的习惯是同时画出理论参考线不考虑频偏、假设模糊已正确消除时相位估计误差近似与 (1/(N \cdot SNR_{eff})) 成正比其中 (SNR_{eff}) 是 M 次方后的等效信噪比。把实际曲线和这条参考线放一起如果高频区偏离超过 2~3 倍说明实现里有 bug多半是角度包裹没处理好如果低频区突然翘头说明进入了相位模糊区独特字长度不够。这个方法帮我排查过好几次模糊消除出错的问题每次都能精准定位到是独特字长度还是候选相位生成的问题。这个算法用了快十年我的感受是它简单、快、不用训练序列但盲的代价就是相位模糊和门限效应这两堵墙躲不掉。遇到这两种问题别硬调参数先确认是哪种墙再决定是加独特字还是加级联。希望帮到你。本文还有配套的精品资源点击获取
