光纤干涉PGC-DCM相位解调：MATLAB实现、参数标定与避坑指南

简介：一套基于PGC相位生成载波调制与微分交叉相乘(DCM)解调算法的MATLAB代码，面向光纤传感、干涉测量领域的算法学习者和信号处理工程师。代码围绕相位载波生成、信号混频、微分交叉相乘、低通滤波等关键环节展开，完整呈现PGC调制与DCM解调的处理逻辑，方便读者理解从原始干涉信号到相位信息的恢复过程。资源包共3个文件，包含两个.m脚本（分别实现DCM解调主流程和低通滤波器设计）与一个对低通滤波器设计过程及要点进行总结的PDF文档，整体压缩包仅382KB，轻量便于离线查阅。目前已有2257人学习或下载，具备较好的参考价值。通过这份资源，读者可获得可直接运行的MATLAB解调示例、滤波器设计思路以及配套说明文档，在仿真中验证PGC载波调制与DCM解调算法的可行性和参数影响，并据此调整载波频率、滤波截止频率等关键参数，向实际应用系统扩展。 做光纤干涉信号处理的朋友，应该都清楚：干涉仪输出的原始信号是一个随相位周期性变化的余弦波，直接对这个波做换算，得到的相位不是折叠在周期里，就是被条纹对比度变化搞得忽大忽小。早期我被这个问题卡了挺久，后来系统地把PGC相位生成载波调制和DCM微分交叉相乘解调整条链路在MATLAB里完整跑通，很多之前模糊的地方才真正清晰起来。这套方法在光纤水听器、光纤传感解调、激光干涉测振这些场景里非常常见，属于绕不开的基础功。

这篇就把我整理好的MATLAB实现思路、参数标定方法和实际踩坑记录一次性写清楚。代码结构我会拆开讲，每个环节为什么这么设计、参数为什么取这个值，都会给到判断依据，方便你拿到之后直接改成自己的仿真或者工程原型。

1. 为什么干涉仪非要绕一个"载波"来解相位

1.1 干涉信号的数学模型与直接解调的困境

先看干涉仪最典型的输出形式。忽略偏振等次要因素，光电探测器上得到的信号可以写成：

I(t) = A + B·cos(φ(t))

其中A是直流偏置，B与入射光功率、干涉条纹可见度有关，φ(t)是待测的外界信号引起的相位差。理论上只要对这个式子做反余弦运算，就能把φ(t)提取出来。但问题在于：

• cos函数在φ接近0或π的整数倍时斜率接近0，反余弦结果对噪声极其敏感，直接解调会看到严重的非线性失真。
• B不是常数，光源波动、光纤扰动都会让B变化，反余弦解调的误差随B剧烈变化。
• 相位本身可能超过2π范围，直接反余弦会遇到周期性折叠。

这就是人们常说的干涉信号衰落与相位模糊。解决思路之一是让系统工作在一个更有利于线性解调的工作点上，但环境扰动会缓慢改变工作点，于是就有了更工程化的做法：主动在干涉仪里注入一个高频相位载波，把待测信号整体搬到高频段，再通过混频和滤波把信号分量提取出来。这就是PGC（Phase Generated Carrier）的基本思想。

1.2 PGC调制把一个"非线性测量问题"变成"线性传输问题"

PGC的做法其实和广播电台调频发射很像。我们不在基带直接解调，而是给干涉仪加上一个幅度为C、频率为fc的高频载波调制项，让干涉输出变成：

I(t) = A + B·cos( C·cos(2π·fc·t) + φ(t) )

其中C称为调制深度，单位是rad。待测信号φ(t)原本让输出在cos的非线性区间里来回摆动，加了载波之后，φ(t)的信息会被调制到fc及其谐波分量上。后续用本地参考信号做混频、低通，就能分别在基频和二倍频处得到携带sinφ(t)和cosφ(t)的缓变分量。

这一招的关键意义在于：它把难以处理的cos反演问题，转化成了两个缓变正交分量的提取问题。后面的DCM解调和PGC-Arctan解调，本质都是在处理这两个正交分量。所以整条链路的第一步，是先理解载波在频域上干了什么事。

2. DCM解调链条：从Bessel函数到微分交叉相乘

2.1 Bessel展开与正交混频

干涉输出中带有cos(C·cosωc·t)和sin(C·cosωc·t)这两类项，它们是周期信号，理论上可以用Bessel函数展开成无穷多个谐波分量。对cos项和sin项分别展开后，可以看到基频分量正比于J1(C)·sinφ(t)，二倍频分量正比于J2(C)·cosφ(t)。数学推导很多教材都有，我这里直接给结论，方便后面写代码时对照概念：

• 用本地参考信号cos(ωc·t)与I(t)混频，再经过低通滤波，得到S1 ≈ -2B·J1(C)·sinφ(t)。
• 用本地参考信号cos(2ωc·t)与I(t)混频，再经过低通滤波，得到S2 ≈ -2B·J2(C)·cosφ(t)。

这里J1(C)、J2(C)分别是一阶和二阶Bessel函数在C处的取值。注意本地参考信号的频率和相位必须与注入的相位载波严格一致，否则混频后会产生额外的低频误差。

2.2 微分交叉相乘：消除信号衰落的关键一步

得到S1和S2之后，最直观的思路是直接反正切：φ(t) = arctan(S1/S2)。这就是PGC-Arctan方案。DCM不这样做，它对S1和S2做微分和交叉相乘：

S1'·S2 - S1·S2'

把S1 = k1·sinφ、S2 = k2·cosφ代入，展开后能看到sinφ、cosφ项全部被三角恒等式消掉，只剩下K·φ'(t)，其中K = 4·B²·J1(C)·J2(C)。再对这个结果做一次积分，就恢复出φ(t)。

这步操作的价值在于：B作为一个公共因子进入增益K，当干涉可见度波动时，K会变化，但解调得到的相位波形形状不会发生周期性衰落或畸变。相比于直接针对cos信号求反演，DCM在B随机波动的场景下有天然优势。代价是输出相位前面乘了一个标定系数K，需要做一次增益标定。

2.3 积分与增益标定

DCM最后一步是积分。MATLAB里用cumtrapz就能完成，但要注意积分结果里面包含K这个增益倍数，实际工程上必须标定：输入一个已知幅度的相位信号φ0，测量解调输出幅度φ_est，然后得到标定系数gain = φ0 / φ_est，后续所有解调结果除以gain即可。

还有一个很容易被忽略的点：J1(C)和J2(C)的大小与调制深度C强相关。C如果漂移，K就会漂移，导致解调结果幅度不准。所以无论是仿真还是工程，都应该先做一遍C值标定，尽量把C稳定在J1(C)和J2(C)都比较大的区间，常见推荐值是C约2.6到3.1 rad之间，经典经验值是2.63 rad。

3. MATLAB工程化实现：参数怎么定、代码怎么组织

3.1 仿真信号生成

先把干涉仪输出仿真出来。为了覆盖常见使用场景，我建议把采样率、载波频率、信号频率、调制深度、干涉可见度参数全部做成变量，方便后面做参数扫描。

%% 参数定义 fs = 5e6; % 采样率 5 MHz fc = 200e3; % 相位生成载波频率 200 kHz f0 = 10e3; % 待测信号频率 10 kHz amp = 1.0; % 待测信号幅度，单位 rad C = 2.63; % 调制深度，单位 rad B = 1.0; % 与干涉条纹可见度相关 DC = 2.0; % 直流偏置 N = fs / 10; % 仿真0.1秒 t = (0:N-1) / fs; %% 生成干涉输出 phi = amp * sin(2*pi*f0*t); carrier = C * cos(2*pi*fc*t); I = DC + B * cos(carrier + phi);

3.2 DCM解调主函数

解调部分我建议封装成独立函数，方便对多组参数重复调用。核心步骤就是混频、低通、微分、交叉相乘、相减、积分。

function phi_est = dcm_demod(I, t, fc, C, B) % DCM解调 % 输入： % I 干涉仪输出信号 % t 时间序列 % fc 载波频率 % C 调制深度 % B 干涉可见度相关项 % 输出： % phi_est 解调得到的相位信号 fs = 1 / (t(2) - t(1)); % 1. 与本地参考混频 x1 = I .* 2 .* cos(2*pi*fc*t); x2 = I .* 2 .* cos(4*pi*fc*t); % 2. 低通滤波提取低频分量 % 截止频率需要低于载波频率，同时高于信号最高频率 fcut = min(60e3, fs/10); % 示例值，需根据信号带宽调整 [b, a] = butter(4, fcut/(fs/2), 'low'); sinphi = filtfilt(b, a, x1); cosphi = filtfilt(b, a, x2); % 3. 微分（用中心差分，避免相位延迟） d_sin = gradient(sinphi) * fs; d_cos = gradient(cosphi) * fs; % 4. 交叉相乘相减 term = d_sin .* cosphi - sinphi .* d_cos; % 5. 积分恢复 phase_scaled = cumtrapz(t, term); % 6. 增益标定 J1 = besselj(1, C); J2 = besselj(2, C); K = 4 * B^2 * J1 * J2; phi_est = phase_scaled / K; end

调用之后，可以直接对比解调输出和真实相位：

phi_est = dcm_demod(I, t, fc, C, B); plot(t*1e3, phi, t*1e3, phi_est);

如果参数合理，两条曲线几乎重合。需要注意，滤波器和微分器都会在信号起始段产生边界效应，观察解调结果时建议丢掉开头一小段再看。

3.3 滤波器设计与微分器的工程坑

低通滤波器参数是整个解调质量的关键。混频之后，我们需要保留的是低频的sinφ、cosφ分量，但同时要抑制fc±信号带宽、2fc±信号带宽这些高频分量。如果LPF的截止频率取得太高，载波泄漏会让解调输出带上周期性纹波；取得太低，则信号的高频成分被衰减，动态范围变差。

我这里用butter(4, fcut/(fs/2))是IIR滤波，配合filtfilt做零相位处理。仿真场景下没问题，但filtfilt是离线处理，只能用于后处理分析。如果要做实时解调，需要换成filter函数，或者用DSP/FPGA上的实时滤波结构，那时候还要额外考虑滤波器群延迟，不能直接套用离线仿真的参数。

微分环节更是一个经典坑。直接对S1、S2做数值微分，会把高频噪声放大——尤其是低通滤波没有完全抑制掉的微弱载波残留，经过微分后会变成明显的毛刺。我的建议是：

• 如果能提高采样率，就尽量提高，让差分近似的精度更好。
• 使用gradient而非diff进行微分，避免信号长度缩短和相位偏移。
• 在实际传感器数据中，可以在微分之后再加一个低通平滑，或者在微分之前使用更高阶的低通滤波器。
• 如果噪声仍然明显，可以考虑用Savitzky-Golay滤波器的导数滤波功能替代纯差分，效果会好很多。

4. 参数实验与解调效果验证

4.1 调制深度C的扫描实验

DCM解调的增益K与J1(C)·J2(C)成正比，因此不同C值下，同一幅度相位信号的解调输出幅度会明显不同。我建议在仿真中做一次C值扫描，把C从0.5到5按0.1步进扫一遍，分别解调相同输入，记录输出幅度。你会看到输出幅度随C变化出现明显起伏，在C约2.63附近有一个峰值区间，低于1或者高于4之后输出幅度显著下降，甚至被噪声淹没。

这组实验的意义是：在搭建实际系统时，你要通过调整注入载波幅度来观察解调输出，找到输出幅度最大且对幅度扰动不敏感的工作点。如果发现输出幅度相对C的斜率过大，说明系统工作点选择不好，对载波幅度漂移过于敏感。

4.2 用外部采集数据验证：把CSV导入MATLAB跑仿真

仿真合拍之后，下一步自然是拿真实采集数据验证。很多朋友习惯直接把示波器或采集卡导出的CSV文件拖进MATLAB，这个流程本身没问题，但有几个细节要注意：

• 确认CSV中的时间列和信号列的单位，示波器导出的电压值如果没除以探测器跨阻增益，B值和量纲就和仿真对不上。
• 确认采样率是否准确，混频参考的频率是绝对频率，采样率偏差会导致解调相位出现斜漂。
• 实测数据普遍带工频干扰和高频噪声，建议先做一次带通滤波预处理，再进入DCM解调。

data = readmatrix('practical_data.csv'); t_real = data(:, 1); I_real = data(:, 2); % 剔除均值、去除趋势 I_real = detrend(I_real, 'constant'); % 进入解调函数前确认fs fs_real = 1 / (t_real(2) - t_real(1));

实测数据跑出来的解调结果不会像仿真那么干净，但波形趋势和幅度量级可以作为系统验证依据。如果解调结果出现明显漂移，优先检查参考信号的频率准确度和低通滤波器的截止频率选择。

5. DCM与PGC-Arctan方案：怎么选、怎么切换

5.1 两种方案的数学对比

DCM之外，PGC-Arctan是另一种常用解调方法。它的做法更加直接：混频低通得到S1和S2之后，直接用atan2(S1, S2)提取相位。两者对比如下：

DCM最大的优势是动态范围天然更大，因为它整个处理链里没有周期性折叠，相位超过π也不会跳变，微分交叉相乘后积分出来还是连续的。Arctan方案虽然实现直观，但反正切结果会被限制在±π/2之内，相位超出这个范围就必须配合相位解缠绕算法，复杂度和误判风险都上来了。

5.2 选型建议

我的实践经验是：如果待测信号本身是小相位、高精度的慢变化信号，PGC-Arctan配合解缠绕往往更稳，因为反正切天然做了归一化，对光强波动不那么敏感。如果是大动态范围的冲击、振动、水声信号，尤其相位可能大幅超过π的场景，DCM会更省心，省掉了相位解缠绕这个麻烦。

另外一个常见误区是觉得DCM不用反正切所以更简单。实际上DCM的增益标定比Arctan麻烦，因为K同时包含B²和Bessel函数项，如果B不稳定，就得经常重新标定。Arctan里面S1和S2做比值后，B的影响会直接消掉，这是它的隐藏优势。

6. 实际调试中的避坑清单

下面这些是我不止一次踩过、也在代际迭代中反复验证过的典型问题，按出现频率从高到低列出来：

1. 微分造成的高频噪声毛刺。可以先用filtfilt把混频低通结果再平滑一次，再微分；或者换用SG导数滤波器。不要试图在微分之后再加大低通，那样相位滞后会很难补偿。

2. 载波频率与信号带宽间隔不够。载波至少取信号最高频率的5到10倍，否则LPF无法在保留信号的同时有效滤除混频杂散分量。比如信号带宽100kHz时，载波最少给到500kHz以上，采样率再留足3到5倍余量。

3. 本地参考相位失配。混频参考信号的相位如果和载波不一致，S1和S2中会混入正交泄漏项，表现为解调结果出现非线性串扰。工程实现时通常需要加锁相环或对参考相位做在线搜索校准。

4. 直流分量滤除不彻底。如果LPF不能完全抑制直流偏置A混频后产生的残余项，积分输出会叠加一个斜坡漂移。这种情况可以在混频之前先对I做一次去直流，或者把LPF阶数提高。

5. 调制深度C被忽略。很多初学者把C当成无关紧要的常数，结果系统换一个光功率后解调幅度明显变化。正确做法是，在系统联调时用已知幅度的标准信号源做C值标定，并把C的工作点记录成系统参数。

6. 滤波器阶数过高引入群延迟。离线仿真无所谓，实时系统里滤波器的群延迟会造成解调相位滞后，评估动态响应时要记得补偿。用Bessel滤波器代替Butterworth可以改善群延迟特性，但阻带衰减会差一些。

整个PGC-DCM链路在MATLAB里跑通之后，后续迁移到嵌入式或者FPGA平台时，可以把解调核心模块直接翻译成定点运算结构。微分、乘法、积分这些操作在FPGA上都有现成IP核，真正需要重新做的大头其实是滤波器系数定标和实时性能优化。就我个人的经验来说，先把MATLAB仿真做扎实，把每个参数的作用和坑都摸清楚，再往下走会省掉非常多返工时间。

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