相位生成载波（PGC）调制与解调讲解

相位⽣成载波（PGC）调制与解调讲解
相位⽣成载波（PGC ）调制与解调
⼀、 PGC 调制
⼲涉型光纤传感器的解调⽅法⽬前主要有：相位⽣成载波解调法、光路匹配差分⼲涉法、差分时延外差法。

由于相位⽣成载波解调信号有动态范围⼤、灵敏度⾼、线性度好、测相精度⾼等优点，是⽬前光纤传感⼲涉领域⼯程上较为实⽤的解调⽅法。

[1]
相位⽣成载波的调制分为外调制和内调制。

外调制⼀般采⽤压电陶瓷（PZT ）作为相位调制器，假设调制信号频率为ω0 ，幅度为C ，调制信号可以表⽰为（1）式：
0(t)cos(t)C φω= （1）
则光纤⼲涉仪的输出的信号可表⽰为（2）式：
00cos[(t)(t)]cos[cos(t)(t)]s s I A B A B C φφωφ=++=++ （2）
式中，A 为直流量， B 为⼲涉信号幅度。

s (t)Dcos(t)(t)s φωψ=+，其中，?s (t) 不仅包含了待测信号D cos ωs t ，还包括了环境噪声引起的相位变化ψ(t)。

将（2）式按 Bessel 函数展开，得到（3）式[2]：
k k 02k 02k 1010J (C)2(1)J (C)cos 2k t cos (t)2(1)J (C)cos(2k 1)t sin (t)s s k k I A B ωφωφ∞
∞+==
=++---+??
∑∑ (3)
⼆、 PGC 解调
微分交叉相乘（differential and cross —multiply ，DCM ）算法和反正切算法是
两种传统的 PGC 解调算法，此外,⽂献[1]中还介绍了三倍频DCM 算法，基频混频PGC 算法，基于反正切算法和基频混频算法的改进算法，反正切-微分⾃相乘算法（Arctan-DSM ）算法。

下⾯分别介绍DCM 算法和反正切算法。

2.1 微分交叉相乘（DCM ）算法 DCM 算法的原理图如图1所⽰：
图1 DCM 算法原理图
输⼊的⼲涉信号I 分别与基频信号10cos S G t ω=和⼆倍频信号20cos 2S H t ω=进⾏混频，再通过低通滤波器滤除⾼频成分，可以得到信号的正弦项（5）式和余弦项（6）
式：
1s (t)(C)sin (t)I BGJ φ=- （5） 2s Q(t)(C)cos (t)BHJ φ=- （6）
I(t) 、Q(t) 含有外界⼲扰，还不能直接提取待测信号，再通过微分交叉相乘（DCM ）⽅法得到两个正交信号的平⽅项，利⽤sin 2?s + cos 2?s = 1消除正交量，得到微分量（7）式：
212s '(C)J (C)'(t)V B GHJ φ= (7)
经过积分运算再通过⾼通滤波器滤除缓慢变化的环境噪声，最终得到的解调信号为得到（8）式：
[]2212s 12s (C)J (C)(t)(C)J (C)Dcos(t)(t)V B GHJ B GHJ φωψ==+ （8）
相位噪声项ψ(t) 通常情况下为缓变信号，将V 通过⾼通滤波器滤除相位噪声，就可以得到待测信号，实现传感信号的解调（9）式。

212s (C)J (C)Dcos(t)out V B GHJ ω= (9)
由式（8）可以看出，最后的解调输出信号与待测信号成线性关系，因此与后⾯将要讨论的反正切算法相⽐，产⽣的⾮线性失真要⼩的多。

但是由于输出信号中的⼲涉幅度B = κA ，⽽κ⼜与光传输中偏振态的变化有关，A 是与光源光功率的稳定度、光路中各环节光功率的衰减、光纤⼲涉仪输⼊的光强等因素有关的量。

因此，解调信号幅度受调制深度、光强、光路损耗、耦合器分光⽐、偏振态等诸多因素的影响。

[1] 2.2 反正切算法
反正切算法的原理如图2所⽰：
图2反正切算法原理图
反正切算法前端部分均与 DCM 法相同，只是该法在两路信号分别通过低通滤波器后，I(t) 、Q(t)进⾏相除得到（10）式[3]：
12(C)
(t)tan (t)(t)(C)
s GJ I Q HJ φ= （10）对式（10）进⾏反正切运算，得到（11）式：
2s 1(t)*(C)
(t)Dcos(t)(t)arctan(
)(t)*(C)
s I HJ Q GJ φωψ=+= （11）
噪声项ψ(t) 通常情况下为缓变信号,通过⾼通滤波器滤除环境噪声ψ(t) 即可得到待测信号Vout=D cos ωs t 。

与 DCM 算法相⽐，反正切算法通过除法运算，消除了B 对解调结果的影响，⽽且，如果令 G=H ,那么G 、H 对解调结果的影响也会被消除。

同时，反正切算法⽐ DCM 算法原理相对简单，使得其解调算法⽐较简单，从⽽缩短了系统信号处理的时间，使系统的实时性得到了显著地提⾼。

但是，由于调制深度（C 值）的偏差，使得J1(C) /J2(C)不等于1，从⽽使解调结果产⽣了⾮线性，同时带来了严重的谐波失真以及总谐波失真。

[1]
C 选取的原则是C 的值尽可能⼩且是的J1(C ）J2(C)的变化趋势趋于0。

当C=2.37时，J1(C ）J2(C)的导数为0，J1(C ）J2(C)取得极⼤值。

因此，对于PGC-DCM 算法的最佳调制度C=2.37.
⽽对于PGC-Arctan 算法来讲，解调的结果是与J1(C ）/J2(C)相关。

当C=2.63时，J1(C ）/J2(C)=1，故对于传统的PGC-Arctan 算法来讲，C=2.63是最佳调制度。

⽆论是PGC-DCM 算法还是PGC-Atctan 算法，C 选取的原则是选取⼀个恰当的值使得解调结果随着与C 相关的贝塞尔函数想的变化趋于稳定。

[4]
根据⽂献[4],低通滤波器的截⽌频率应该满⾜：0/2pass f f <。

D 的取值范围应该满⾜0(1)/2s D ωω+<。

三、 MATLAB 仿真
3.1 PGC 调制
设置采样频率为4MHz ，采样点数为10K 。

调制信号： C=2.37;%调制度
f0=50000;%载波频率50KHz
s1=C*cos(2*pi*f0/Fs*t);%调制信号待测信号：
D=1.2;%待测信号幅度
fs=1200;%待测信号频率1.2KHz fn=10;%假设噪声频率10Hz sn=0.1;%噪声信号
s2=D*cos(2*pi*fs/Fs*t);%待测信号调制后的信号： A=1; B=1;
s3=s2+sn;
I=A+B*(cos(C*cos(2*pi*f0/Fs*t)+s3));%调制后的信号
s1，s2，s3，I 对应的波形如图3：
图3 信号波形
3.2 微分交叉相乘（DCM ）算法解调信号分别乘以基频信号和2倍频信号 G=5; H=5;
mod1=G*cos(2*pi*f0/Fs*t); mod2=H*cos(2*pi*f0*2/Fs*t); x1=I.*mod1;%与基频信号相乘 x2=I.*mod2;%与2倍频信号相乘低通滤波 % x1低通 % x2低通
LP=lowpassfilter;
x1_out=conv(x1,LP.Numerator);%I(t) x2_out=conv(x2,LP.Numerator);%Q(t) 微分
x11=diff(x1_out); x22=diff(x2_out); L=length(x1_out); DCM
x_d=(x2_out(1:L-1)).*x11-(x1_out(1:L-1)).*x22; %积分
for j=1:N-1
2000
4000
6000800010000
12000
-50
5(a)调制信号s1
2000
4000
60008000
10000
12000
-20
2(b)待测信号s2
2000
4000
60008000
10000
12000
-20
2(c)信号s3
500
1000
1500
2000
2500
3000
1
2(d)调制后的信号I
X(j)=0; end
for i=2:L-1
X(i)=X(i-1)+(x_d(i)+x_d(i-1))/2; %积分 end 去除低频信号 for i=2:L-1
XX(i)=X(i)-mean(X); end
解调前后的波形对⽐如图4所⽰：
图4 解调前后的波形信号
3.3 反正切算法解调信号载波信号： C=2.63;%调制度
f0=50000;%载波频率50KHz
s1=C*cos(2*pi*f0/Fs*t);%调制信号
s1，s2，s3，I 对应的波形如图5：
20004000600080001000012000
-1.5-1
-0.5
0.5
1
1.5
图5 信号波形
解调前后的波形对⽐如图6所⽰，由于使⽤卷积运算，信号的前后有失真，长度为滤波器长度的⼀半Number/2（Number=256）。

图6 解调前后的波形信号
2000
4000
6000
800010000
12000
-50
5(a)调制信号s1
2000
4000
60008000
10000
12000
-20
2(b)待测信号s2
2000
4000
6000
800010000
12000
-20
2(c)信号s3
500
1000
1500
2000
2500
3000
1
2(d)调制后的信号I
20004000600080001000012000
-1.5-1
-0.5
0.5
1
1.5
参考⽂献:
[1] 张爱玲，王恺晗，等. ⼲涉型光纤传感器PGC解调算法的研究[J]. 光电技术应⽤，2013，28 （6）：49-52.
[2] 夏东明，娄淑琴，等. ⼲涉型光纤传感器相位载波解调技术研究[J]. 光电技术应⽤，2011，26 （5）：47-50.
[3] Gaosheng Fang, Tuanwei Xu. Phase-sensitive Optical Time Domain Reflectometer Based on Phase Generated Carrier
Algorithm [J]. Journal of Lightwave Technology, 2015.
[4] 王燕. ⼲涉型光纤传感器及PGC解调技术研究[D]. 天津: 天津理⼯⼤学, 2013.
附录：
1. Bessel 函数展开
n 02n 1
cos(zcos )J (z)2(1)J (z)cos(2n )n θθ∞
==+-∑
n 2n 11
sin(zcos )2(1)J (z)cos[(2n 1)]n θθ∞
-==---∑
02n 1
cos(z sin )J (z)2J (z)cos(2n )n θθ∞
==+∑
2n 11
sin(zsin )2J (z)sin[(2n 1)]n θθ∞
-==-∑
2. PGC 调制+ DCM 算法解调信号(MATLAB 程序) %%PGC 调制% clc;
clear all; close all;
Fs=4000000;%采样率4MHz N=10*1024;%采样点数10K t=1:1:N;
C=2.37;%调制度
f0=50000;%载波频率50KHz
s1=C*cos(2*pi*f0/Fs*t);%调制信号
D=1.2;%待测信号幅度
fs=1200;%待测信号频率1.2KHz fn=10;%假设噪声频率10Hz
%sn=0.05*cos(2*pi*fn/Fs);%噪声信号 sn=0.1;%噪声信号
s2=D*cos(2*pi*fs/Fs*t);%待测信号
%调制 A=1; B=1;
s3=s2+sn;
I=A+B*(cos(C*cos(2*pi*f0/Fs*t)+s3));%调制后的信号 figure(1);
subplot(4,1,1),plot(s1);title('(a)调制信号s1'); subplot(4,1,2),plot(s2);title('(b)待测信号s2'); subplot(4,1,3),plot(s3);title('(c)信号s3'); subplot(4,1,4),plot(I),axis([1,3000,-0.2,2.2]);title('(d)调制后的信号I');
%分别乘以基频信号和2倍频信号
G=5;
H=5;
mod1=G*cos(2*pi*f0/Fs*t);
mod2=H*cos(2*pi*f0*2/Fs*t);
x1=I.*mod1;%与基频信号相乘
x2=I.*mod2;%与2倍频信号相乘
% x1低通
% x2低通
LP=lowpassfilter;
x1_out=conv(x1,LP.Numerator);%I(t)
x2_out=conv(x2,LP.Numerator);%Q(t)
%微分
x11=diff(x1_out);
x22=diff(x2_out);
L=length(x1_out);
%DCM
x_d=(x2_out(1:L-1)).*x11-(x1_out(1:L-1)).*x22;
%积分
for j=1:N-1
X(j)=0;
end
for i=2:L-1
X(i)=X(i-1)+(x_d(i)+x_d(i-1))/2; %积分
end
%去除低频信号
for i=2:L-1
XX(i)=X(i)-mean(X);
end
figure(2);
plot(XX./max(XX),'r');hold on;
plot(s3,'b');
legend('解调后的信号','s3');
3. PGC调制+ 反正切算法解调信号(MATLAB程序) %%PGC调制% clc;
clear all;
close all;
Fs=4000000;%采样率4MHz
N=10*1024;%采样点数10K
t=1:1:N;
C=2.63;%调制度
f0=50000;%载波频率50KHz
s1=C*cos(2*pi*f0/Fs*t);%调制信号
D=1.2;%待测信号幅度
fs=1200;%待测信号频率1.2KHz
fn=10;%假设噪声频率10Hz
%sn=0.05*cos(2*pi*fn/Fs);%噪声信号
sn=0.1;%噪声信号
s2=D*cos(2*pi*fs/Fs*t);%待测信号
%调制
A=1;
B=1;
s3=s2+sn;
I=A+B*(cos(C*cos(2*pi*f0/Fs*t)+s3));%调制后的信号
figure(1);
subplot(4,1,1),plot(s1);title('(a)调制信号s1');
subplot(4,1,2),plot(s2);title('(b)待测信号s2');
subplot(4,1,3),plot(s3);title('(c)信号s3');
subplot(4,1,4),plot(I),axis([1,3000,-0.2,2.2]);title('(d)调制后的信号I'); %分别乘以基频信号和2倍频信号
G=5;
H=5;
mod1=G*cos(2*pi*f0/Fs*t);
mod2=H*cos(2*pi*f0*2/Fs*t);
x1=I.*mod1;%与基频信号相乘
x2=I.*mod2;%与2倍频信号相乘
% x1低通
% x2低通
LP=lowpassfilter;
x1_out=conv(x1,LP.Numerator);%I(t)
x2_out=conv(x2,LP.Numerator);%Q(t)
L=length(x1_out);
for i=1:L
XX(i)=atan(x1_out(i)/x2_out(i));
end
figure(2);
plot(XX./max(XX),'r');hold on;
plot(s3,'b');
legend('解调后的信号','s3');
4. 低通滤波器设计（使⽤fdatool设计滤波器）
完成之后点击 File->Generate MATLAB code ；⽣成滤波器lowpassfilter.m⽂件。

.M⽂件内容：function Hd = lowpassfilter
%LOWPASSFILTER Returns a discrete-time filter object.
%
% MATLAB Code
% Generated by MATLAB(R) 7.13 and the Signal Processing Toolbox 6.16.
%
% Generated on: 27-Oct-2016 16:10:36
%
% FIR Window Lowpass filter designed using the FIR1 function.
% All frequency values are in Hz.
Fs = 10000000; % Sampling Frequency
N = 512; % Order
Fc = 20000; % Cutoff Frequency
flag = 'scale'; % Sampling Flag
% Create the window vector for the design algorithm. win = hamming(N+1); % Calculate the coefficients using the FIR1 function.
b = fir1(N, Fc/(Fs/2), 'low', win, flag);
Hd = dfilt.dffir(b);
% [EOF]。