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

《相位生成载波(PGC)调制与解调(包含MATLAB代码).pdf》由会员分享，可在线阅读，更多相关《相位生成载波(PGC)调制与解调(包含MATLAB代码).pdf（12页珍藏版）》请在三一文库上搜索。

1、相位生成载波相位生成载波（PGC）调制与）调制与解调解调 一、一、 PGCPGC 调制调制 干涉型光纤传感器的解调方法目前主要有：相位生成载波解调法、光路匹配差分干 涉法、差分时延外差法。由于相位生成载波解调信号有动态范围大、灵敏度高、线性度 好、测相精度高等优点，是目前光纤传感干涉领域工程上较为实用的解调方法。 1 相位生成载波的调制分为外调制和内调制。外调制一般采用压电陶瓷（PZT）作为 相位调制器，假设调制信号频率为0 ，幅度为C，调制信号可以表示为（1）式： 0 (t)cos(t)C （1） 则光纤干涉仪的输出的信号可表示为（2）式： 00 cos(t)(t)coscos(t)(t)

2、ss IABABC （2） 式中，A 为直流量， B 为干涉信号幅度。 s (t)Dcos(t)(t) s ，其中，s(t) 不 仅包含了待测信号D cos st ，还包括了环境噪声引起的相位变化(t)。 将（2）式按 Bessel 函数展开，得到（3）式 2： kk 02k02k 10 10 J (C)2( 1) J(C)cos2kt cos(t)2( 1) J(C)cos(2k 1)t sin(t) ss kk IAB (3) 二、二、 PGCPGC 解调解调 微分交叉相乘（differential and crossmultiply，DCM）算法和反正切算法是 两种传统的 PGC 解调算

3、法，此外,文献1中还介绍了三倍频 DCM 算法，基频混频 PGC 算法，基于反正切算法和基频混频算法的改进算法，反正切-微分自相乘算法 （Arctan-DSM）算法。下面分别介绍 DCM 算法和反正切算法。 2.12.1 微分交叉相乘（微分交叉相乘（DCMDCM）算法）算法 DCM 算法的原理图如图 1 所示： 图 1 DCM 算法原理图 输入的干涉信号 I 分别与基频信号 10 cosSGt和二倍频信号 20 cos2SHt进 行混频， 再通过低通滤波器滤除高频成分， 可以得到信号的正弦项 （5） 式和余弦项 （6） 式： 1s (t)(C)sin(t)IBGJ （5） 2s Q(t)(C)

4、cos(t)BHJ （6） I(t) 、Q(t) 含有外界干扰，还不能直接提取待测信号，再通过微分交叉相乘（DCM） 方法得到两个正交信号的平方项， 利用sin2s + cos2s = 1消除正交量， 得到微分量 （7） 式： 2 12s (C)J (C) (t)VB GHJ (7) 经过积分运算再通过高通滤波器滤除缓慢变化的环境噪声， 最终得到的解调信号为 得到（8）式： 22 12s12s (C)J (C) (t)(C)J (C) Dcos(t)(t)VB GHJB GHJ （8） 相位噪声项 (t) 通常情况下为缓变信号，将V通过高通滤波器滤除相位噪声，就 可以得到待测信号，实现传感信号

5、的解调（9）式。 2 12s (C)J (C)Dcos(t) out VB GHJ (9) 由式（9）可以看出，最后的解调输出信号与待测信号成线性关系，因此与后面将 要讨论的反正切算法相比，产生的非线性失真要小的多。但是由于输出信号中的干涉幅 度B = A ，而 又与光传输中偏振态的变化有关，A 是与光源光功率的稳定度、光路 中各环节光功率的衰减、光纤干涉仪输入的光强等因素有关的量。因此，解调信号幅度 受调制深度、光强、光路损耗、耦合器分光比、偏振态等诸多因素的影响。 1 2.22.2 反正切算法反正切算法 反正切算法的原理如图 2 所示： 图 2 反正切算法原理图 反正切算法前端部分均与 D

6、CM 法相同，只是该法在两路信号分别通过低通滤波器 后，I(t) 、Q(t)进行相除得到（10）式 3： 1 2 (C)(t) tan(t) (t)(C) s GJI QHJ （10） 对式（10）进行反正切运算，得到（11）式： 2 s 1 (t)*(C) (t)Dcos(t)(t)arctan() (t)*(C) s IHJ QGJ （11） 噪声项 (t) 通常情况下为缓变信号,通过高通滤波器滤除环境噪声(t) 即可得到 待测信号Vout=D cos st。 与 DCM 算法相比，反正切算法通过除法运算，消除了 B 对解调结果的影响，而且， 如果令 G=H,那么G、H对解调结果的影响也会

7、被消除。同时，反正切算法比 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

8、(C)相关的。当 C=2.63 时，J1(C）/J2(C)=1，故对于传统的 PGC-Arctan 算法来讲，C=2.63 是最佳调制度。无 论是 PGC-DCM 算法还是 PGC-Atctan 算法， C 选取的原则是选取一个恰当的值使得解调结 果随着与 C 相关的贝塞尔函数想的变化趋于稳定。 4 根据文献4,低通滤波器的截止频率应该满足： 0/2pass ff。D 的取值范围应该 满足 0 (1)/2 s D。 三、三、 MATLABMATLAB 仿真仿真 3.13.1 PGCPGC 调制调制 设置采样频率为 4MHz，采样点数为 10K。 调制信号：调制信号： C=2.37;%调制度 f

9、0=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.23.2 微分交叉相乘（微分交叉相乘（DCMDCM）算法）算法

10、解调信号解调信号 分别乘以基频信号和分别乘以基频信号和 2 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)

11、; L=length(x1_out); DCMDCM x_d=(x2_out(1:L-1).*x11-(x1_out(1:L-1).*x22; %积分 for j=1:N-1 020004000600080001000012000 -5 0 5 (a)调 制 信 号 s1 020004000600080001000012000 -2 0 2 (b)待 测 信 号 s2 020004000600080001000012000 -2 0 2 (c)信 号 s3 50010001500200025003000 0 1 2 (d)调 制 后 的 信 号 I X(j)=0; end for i=2:L-

12、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.33.3 反正切算法反正切算法解调信号解调信号 载波信号：载波信号： C=2.63;%调制度 f0=50000;%载波频率 50KHz s1=C*cos(2*pi*f0/Fs*t);%调制信号 s1，s2，s3，I 对应的波形如图 5： 020004000600080001000012000 -1.5 -1 -0.5 0 0.5 1 1.5 解 调

13、 后 的 信 号 s3 图 5 信号波形 解调前后的波形对比如图 6 所示，由于使用卷积运算，信号的前后有失真，长度为滤波 器长度的一半 Number/2（Number=512） 。 图 6 解调前后的波形信号 020004000600080001000012000 -5 0 5 (a)调 制 信 号 s1 020004000600080001000012000 -2 0 2 (b)待 测 信 号 s2 020004000600080001000012000 -2 0 2 (c)信 号 s3 50010001500200025003000 0 1 2 (d)调 制 后 的 信 号 I 0200

14、04000600080001000012000 -1.5 -1 -0.5 0 0.5 1 1.5 解 调 后 的 信 号 s3 参考文献: 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 Genera

15、ted 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 1 1 sin(zcos )2( 1) J(z)cos(2n 1) n 02n 1 cos(zsin )J (z)2J(z)cos(2n ) n 2n 1 1 sin(zsin )2J(z)sin(2n 1) n 2. PGC 调制+ DCM

16、 算法解调信号(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);%待测信号 %调制

17、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

18、(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; %

19、积分 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;%采样点数 1

20、0K 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

21、(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 低

22、通 % 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 ；生成滤波

23、器 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

24、. % 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