相位生成载波(PGC)调制与解调(包含MATLAB代码)

相位生成载波（PGC）调制与解调一、PGC调制干涉型光纤传感器的解调方法目前主要有：相位生成载波解调法、光路匹配差分干涉法、差分时延外差法。由于相位生成载波解调信号有动态范围大、灵敏度高、线性度好、测相精度高等优点，是目前光纤传感干涉领域工程上较为实用的解调方法。[1]相位生成载波的调制分为外调制和内调制。外调制一般采用压电陶瓷（PZT）作为相位调制器，假设调制信号频率为ω0，幅度为C，调制信号可以表示为（1）式：0(t)cos(t)C（1）则光纤干涉仪的输出的信号可表示为（2）式：00cos[(t)(t)]cos[cos(t)(t)]ssIABABC（2）式中，A为直流量，B为干涉信号幅度。s(t)Dcos(t)(t)s，其中，ϕs(t)不仅包含了待测信号Dcosωst，还包括了环境噪声引起的相位变化ψ(t)。将（2）式按Bessel函数展开，得到（3）式[2]：kk02k02k1010J(C)2(1)J(C)cos2ktcos(t)2(1)J(C)cos(2k1)tsin(t)sskkIAB(3)二、PGC解调微分交叉相乘（differentialandcross—multiply，DCM）算法和反正切算法是两种传统的PGC解调算法，此外,文献[1]中还介绍了三倍频DCM算法，基频混频PGC算法，基于反正切算法和基频混频算法的改进算法，反正切-微分自相乘算法（Arctan-DSM）算法。下面分别介绍DCM算法和反正切算法。2.1微分交叉相乘（DCM）算法DCM算法的原理图如图1所示：图1DCM算法原理图输入的干涉信号I分别与基频信号10cosSGt和二倍频信号20cos2SHt进行混频，再通过低通滤波器滤除高频成分，可以得到信号的正弦项（5）式和余弦项（6）式：1s(t)(C)sin(t)IBGJ（5）2sQ(t)(C)cos(t)BHJ（6）I(t)、Q(t)含有外界干扰，还不能直接提取待测信号，再通过微分交叉相乘（DCM）方法得到两个正交信号的平方项，利用sin2ϕs+cos2ϕs=1消除正交量，得到微分量（7）式：212s'(C)J(C)'(t)VBGHJ(7)经过积分运算再通过高通滤波器滤除缓慢变化的环境噪声，最终得到的解调信号为得到（8）式：2212s12s(C)J(C)(t)(C)J(C)Dcos(t)(t)VBGHJBGHJ（8）相位噪声项ψ(t)通常情况下为缓变信号，将V通过高通滤波器滤除相位噪声，就可以得到待测信号，实现传感信号的解调（9）式。212s(C)J(C)Dcos(t)outVBGHJ(9)由式（9）可以看出，最后的解调输出信号与待测信号成线性关系，因此与后面将要讨论的反正切算法相比，产生的非线性失真要小的多。但是由于输出信号中的干涉幅度B=κA，而κ又与光传输中偏振态的变化有关，A是与光源光功率的稳定度、光路中各环节光功率的衰减、光纤干涉仪输入的光强等因素有关的量。因此，解调信号幅度受调制深度、光强、光路损耗、耦合器分光比、偏振态等诸多因素的影响。[1]2.2反正切算法反正切算法的原理如图2所示：图2反正切算法原理图反正切算法前端部分均与DCM法相同，只是该法在两路信号分别通过低通滤波器后，I(t)、Q(t)进行相除得到（10）式[3]：12(C)(t)tan(t)(t)(C)sGJIQHJ（10）对式（10）进行反正切运算，得到（11）式：2s1(t)*(C)(t)Dcos(t)(t)arctan()(t)*(C)sIHJQGJ（11）噪声项ψ(t)通常情况下为缓变信号,通过高通滤波器滤除环境噪声ψ(t)即可得到待测信号Vout=Dcosωst。与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/2passff。D的取值范围应该满足0(1)/2sD。三、MATLAB仿真3.1PGC调制设置采样频率为4MHz，采样点数为10K。调制信号：C=2.37;%调制度f0=50000;%载波频率50KHzs1=C*cos(2*pi*f0/Fs*t);%调制信号待测信号：D=1.2;%待测信号幅度fs=1200;%待测信号频率1.2KHzfn=10;%假设噪声频率10Hzsn=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);DCMx_d=(x2_out(1:L-1)).*x11-(x1_out(1:L-1)).*x22;%积分forj=1:N-1020004000600080001000012000-505(a)调制信号s1020004000600080001000012000-202(b)待测信号s2020004000600080001000012000-202(c)信号s350010001500200025003000012(d)调制后的信号IX(j)=0;endfori=2:L-1X(i)=X(i-1)+(x_d(i)+x_d(i-1))/2;%积分end去除低频信号fori=2:L-1XX(i)=X(i)-mean(X);end解调前后的波形对比如图4所示：图4解调前后的波形信号3.3反正切算法解调信号载波信号：C=2.63;%调制度f0=50000;%载波频率50KHzs1=C*cos(2*pi*f0/Fs*t);%调制信号s1，s2，s3，I对应的波形如图5：020004000600080001000012000-1.5-1-0.500.511.5解调后的信号s3图5信号波形解调前后的波形对比如图6所示，由于使用卷积运算，信号的前后有失真，长度为滤波器长度的一半Number/2（Number=512）。图6解调前后的波形信号020004000600080001000012000-505(a)调制信号s1020004000600080001000012000-202(b)待测信号s2020004000600080001000012000-202(c)信号s350010001500200025003000012(d)调制后的信号I020004000600080001000012000-1.5-1-0.500.511.5解调后的信号s3参考文献:[1]张爱玲，王恺晗，等.干涉型光纤传感器PGC解调算法的研究[J].光电技术应用，2013，28（6）：49-52.[2]夏东明，娄淑琴，等.干涉型光纤传感器相位载波解调技术研究[J].光电技术应用，2011，26（5）：47-50.[3]GaoshengFang,TuanweiXu.Phase-sensitiveOpticalTimeDomainReflectometerBasedonPhaseGeneratedCarrierAlgorithm[J].JournalofLightwaveTechnology,2015.[4]王燕.干涉型光纤传感器及PGC解调技术研究[D].天津:天津理工大学,2013.附录：1.Bessel函数展开n02n1cos(zcos)J(z)2(1)J(z)cos(2n)nn2n11sin(zcos)2(1)J(z)cos[(2n1)]n02n1cos(zsin)J(z)2J(z)cos(2n)n2n11sin(zsin)2J(z)sin[(2n1)]n2.PGC调制+DCM算法解调信号(MATLAB程序)%%PGC调制%clc;clearall;closeall;Fs=4000000;%采样率4MHzN=10*1024;%采样点数10Kt=1:1:N;C=2.37;%调制度f0=50000;%载波频率50KHzs1=C*cos(2*pi*f0/Fs*t);%调制信号D=1.2;%待测信号幅度fs=1200;%待测信号频率1.2KHzfn=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);%DCMx_d=(x2_out(1:L-1)).*x11-(x1_out(1:L-1)).*x22;%积分forj=1:N-1X(j)=0;endfori=2:L-1X(i)=X(i-1)+(x_d(i)+x_d(i-1))/2;%积分end%去除低频信号fori=2:L-1XX(i)=X(i)-mean(X);endfigure(2);plot(XX./max(XX),'r');holdon;plot(s3,'b');legend('解调后