简介：基于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_{k=1}^{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_{k=1}^{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 越大，去调制后信噪比损失越惨。这也是为什么 QPSK（M=4）用这个算法效果不错，而 16PSK、32PSK 用这个算法几乎没法工作——高阶调制下噪声放大太狠，块内平均根本压不回来。

还有一个被放大的因素是残余频偏。从上面的推导能看到，残余频偏在 M 次方后变成了原来的 M 倍。本来残余频偏小到几十赫兹，M 次方后变成几百赫兹，块长一长，块内的相位转一圈，平均结果直接趋近于零，估计值变成随机数。这个特性直接影响下一节要讲的块长选择。

所以这个算法的适用边界很清楚：M-PSK 调制（通常是 M=2 或 M=4），残余频偏已经做过粗校准，信噪比在高斯信道下大约高于 0~5dB 的场景。超出这个范围，要么结果方差爆炸，要么相位模糊误判，都需要其他方法来兜底。

3. 用 Python 从零实现 Viterbi-Viterbi 相位估计

3.1 生成带相位偏移的 QPSK 测试信号

先把测试平台搭起来。用 NumPy 生成一段 QPSK 基带符号，加载波相位偏移和复高斯噪声，这样后续每个环节都能对照真实相位验证。

import numpy as np def generate_qpsk_symbols(n_symbols, seed=42): """生成 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, seed=42): """叠加固定相位偏移和复高斯白噪声""" 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 为例，M=4 时等效频偏是 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 用 M=4，BPSK 用 M=2，8PSK 用 M=8。

如果把 M 设成了实际调制阶数的倍数，比如 QPSK 信号用 M=8，会发生什么？QPSK 符号的四次方是 1，但八次方就不是了。设符号相位是 (45^\circ)，八次方相位是 (360^\circ)，OK；但符号相位是 (135^\circ) 时，八次方相位就是 (1080^\circ)，取相位后是 (0^\circ) 和 (0^\circ)？仔细算一下：(135^\circ \times 8 = 1080^\circ)，减去 (2 \times 360^\circ=720^\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)。看起来 M=8 对 QPSK 也能去调制。

那如果把 M 设小了？QPSK 用 M=2，平方之后星座变成两个点在 0 和 180 度交替，调制信息并没有被完全去除，平均结果在两个方向上互相拉扯，相位估计会在两个值之间跳。

结论：M 是调制阶数的整数倍时，去调制依然有效，但相位估计范围缩小了。M=8 的估计范围变成 ([-22.5^\circ, +22.5^\circ])，需要消除的模糊数从 4 重变成 8 重，独特字判决的出错概率也变大。M=4 是 QPSK 的最佳选择，范围合适、模糊数少、噪声放大最小。BPSK 则要压到 M=2。

4.3 滑动平均 vs 块处理：流式场景怎么选

前面代码实现的是块处理：把一段符号切成不重叠的块，每块出一个相位估计值。这在突发通信里很自然，来一个 burst 估一次，burst 内部共享同一个相位。但连续传输场景下，逐符号的相位变化也需要跟踪，这个时候要换成滑动平均结构。

滑动平均的代码逻辑是把前面的核心函数应用在一个滑窗上：

def vv_sliding_phase_estimate(rx, M, window_len): """滑窗式 Viterbi-Viterbi，输出每个符号位置的相位估计""" n = len(rx) est = np.zeros(n, dtype=float) 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（多半是角度包裹没处理好）；如果低频区突然翘头，说明进入了相位模糊区，独特字长度不够。这个方法帮我排查过好几次模糊消除出错的问题，每次都能精准定位到是独特字长度还是候选相位生成的问题。

这个算法用了快十年，我的感受是：它简单、快、不用训练序列，但"盲"的代价就是相位模糊和门限效应这两堵墙躲不掉。遇到这两种问题别硬调参数，先确认是哪种墙，再决定是加独特字还是加级联。希望帮到你。

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