窄带信号包络及相位检波

1、随机实验窄带信号及包络和相位检波分析孙浩01101365郑博年01101379李彪2012年11月23日窄带信号及包络和相位检波分析I：摘要 当窄带系统(接收机)的输入噪声(如热噪声)的功率谱分布在足够宽的频带(相对于接收机带宽)上时，系统的输出即为窄带过程。对于一个窄带信号，通过包络检波器之后，在检波器的输出端可以得到包络A（t）。当窄带信号通过一个相位检波器之后，可以得到有关相位的信息。论文通过用matlab软件对窄带信号的包络和相位，以及窄带信号的数字特征、概率密度、功率谱密度等进行了画图分析，进一步研究了窄带信号的包络和相位的特性。II：实验背景与目的 通过实验掌握窄带随机信号的特点，

2、关键在于包络和相位检波分析。分析并了解了解窄带信号的特性，包括均值、均方值、方差、相关函数、概率密度、频谱及功率谱密度等。熟悉运用常用的信号处理仿真软件平台：matlab软件。III：窄带信号的产生与验证一、窄带随机过程一个实平稳随机过程X(t)，若它的功率谱密度Sx（）具有下述性质：Sx=Sx c- 0c+ 00 others 中心频率为c，带宽为=20，当<<c 时，就可认为满足窄带条件。若随机过程的功率谱满足该条件则称为窄带随机过程。若带通滤波器的传输函数满足该条件则称为窄带滤波器。随机过程通过窄带滤波器传输之后变成窄带随机过程。图1 为典型窄带随机过程的功率谱密度图。若用一

3、示波器来观测次波形，则可看到，它接近于一个正弦波，但此正弦波的幅度和相位都在缓慢地随机变化，图2所示为窄带随机过程的一个样本函数。图1 典型窄带随机过程的功率谱密度图图2 窄带随机过程的一个样本函数二、窄带随机过程的数学表示1、用包络和相位的变化表示由窄带条件可知,窄带过程是功率谱限制在c附近的很窄范围内的一个随机过程,从示波器观察(或由理论上可以推知):这个过程中的一个样本函数(一个实现)的波形是一个频率为c且幅度和相位都做缓慢变化的余弦波。写成包络函数和随机相位函数的形式：X(t)=A(t)*cosct+ (t)其中:A(t)称作X(t)的包络函数; (t)称作X(t)的随机相位函数。包络

4、随时间做缓慢变化，看起来比较直观，相位的变化，则看不出来。2、莱斯（Rice）表示式任何一个实平稳随机过程X(t)都可以表示为：X(t)=Ac(t) cosct-AS(t) sinct其中同相分量：Ac(t)= X(t) cost= X(t) cosct Xtsinct=LPX(t) *2cosct正交分量：AS(t) = X(t)sint=Xt cosct X(t) sinct= LP-X(t) *2sinct（Xt为Xt的希尔伯特变换。LPA表示取A的低频部分）。Ac(t)和AS(t)都是实随机过程，均值为0，方差等于X(t)的方差。三、窄带随机过程仿真建模要求1、根据要求得到X(t)的表

5、达式：x= (l+a) .*cos (2*pi*4000*t+2) +noisy/10;其中：noisy为高斯白噪声，由wgn函数生成，a=cos (2*pi*l000*t)，均值：Ex=mean (x)，方差：Dx=var (x)，计算可得：X(t)的均值为0.0019，X(t)的方差为0.7590。如图4所示，其中蓝色线为X(t)一个样本的时域波形，红色点连成的线为X(t)的均值，绿色点连成的线为X(t)的方差。图4 窄带随机信号时域波形2、求X(t)的概率密度，方法是将最大最小区间分成14等份，然后分别计算各个区间的个数，如图02中柱形条所示，利用曲线拟合， 得到合适的概率密度函数。 为

6、了得到光滑的曲线，利用了多项式拟合，经过测试，9次拟合曲线比较符合要求，获得的曲线如图5中曲线所示：图5 X(t)的概率分布密度函数3、对X(t)进行频谱分析，在Matlab中，利用fft函数可以很方便得求得X(t)的频谱，然后用abs和angle函数求得幅值和相位，画出图像如图6-1所示，图6-2为无噪声叠加时X(t)的频谱图：图6 -1X(t)的频谱图（有噪声叠加）图6 -2X(t)的频谱图（无噪声叠加） 4、求X(t)的自相关函数，用xcorr函数求出自相关序列，得到X(t)自相关函数的时域波形，如图7-1所示，图7-2为无噪声叠加时X(t)的自相关函数图。图7 -1X(t)自相关函数的

7、时域波形（有噪声叠加）图7-2 X(t)自相关函数的时域波形（无噪声叠加）5、对X(t)自相关函数进行fft变换，得到X(t)的功率谱密度，如图8-1所示。图8-2为无噪声叠加时X(t)的功率谱密度图：图8 -1X(t)的功率谱密度（有噪声叠加）图8 -2X(t)的功率谱密度（无噪声叠加）6、建立滤波器，建立一个巴特沃思滤波器，对产生的x(t)进行检测。滤波器的幅度谱和相位谱所示：图9 地通滤波器的幅度谱和相位谱7、求Ac(t)的统计特性，Ac(t)为X(t) *2cosct通过低通滤波器的信号，Ac(t)的均值Eh = -0.4075 4（带有直流分量），Ac(t)的均方值是E2h =0.2

8、458Ac(t)的方差Dh = 0.0798Ac(t)的波形如图10、图11所示：图10 -1Ac(t)的时域波形图和频谱图（有噪声叠加）图10-2 Ac(t)的时域波形图和频谱图（无噪声叠加）图11-1 Ac(t)的自相关函数的时域波形图和Ac(t)的功率谱密度（有噪声叠加）图11 -2Ac(t)的自相关函数的时域波形图和Ac(t)的功率谱密度（无噪声叠加）8、求AS(t)的统计特性，AS(t)为X(t) *2cosct通过低通滤波器的信号，AS(t)的均值Eh =0.8972（带有直流分量），AS(t)的均方值是E2h = 1.1565AS(t)的方差Dh = 0.3518 AS(t)的波

9、形如图13、图14所示：图13-1 AS(t)的时域波形图和频谱图（有噪声叠加）图13 -2AS(t)的时域波形图和频谱图（无噪声叠加）图14-1 AS(t)的自相关函数的时域波形图和AS(t)的功率谱密度（有噪声叠加）图14 -2AS(t)的自相关函数的时域波形图和AS(t)的功率谱密度（无噪声叠加）9、求出Y(t)的统计特性，Y (t)=Ac(t) cosct-AS(t) sinct，其统计特性如下输出信号Y(t)的均值Eh = -4.4011e-004s输出信号Y(t)的均方值E2h = 3.0280输出信号Y(t)的方差Dh = 3.0303Y(t)的仿真图形如图15、图16所示。图1

10、5-1 Y(t)的时域波形图和频谱图（有噪声叠加）图15-2 Y(t)的时域波形图和频谱图（无噪声叠加）图16 -1Y(t)的自相关函数的时域波形图和Y(t)的功率谱密度（有噪声叠加）图16 -2Y(t)的自相关函数的时域波形图和Y(t)的功率谱密度（无噪声叠加）IV：调幅调相波的相干解调在实际应用中，常常需要检测出包络和的信息。若将窄带随机过程X(t)送入包络检波器，则在检波器的输出端可得到包络；若将窄带随机过程X(t)送入一个相位检波器，便可检测出相位信息。如图10所示： 窄带信号及包络和相位检波器图中，在相位检波器之前加入一个理想限幅器，其作用是消除包络起伏对相位检波器的影响。一、 带通

11、滤波器：此带通滤波器的作用是滤除已调信号频带以外的噪声。带通滤波器的参数如下：通带截止频率上下限：2800HZ,5000HZ阻带截止频率上下线：2500HZ,5200HZ阻带衰减：30DB通带衰减：0.1DB该带通滤波器的幅频和相频特性曲线如图9所示：二、 信号通过带通滤波器后的幅频特性：附件：clcfs=16000; %设定采样频率N=1300;n=0:N-1; %取的样本点数t=n/fs; %获得以1/16000为时间间隔采样序列noisy=wgn(1,N,0); %产生高斯白噪声w1=2*pi*1000;%低频信号角频率w2=2*pi*4000;%载波角频率a=cos(2*pi*1000

12、*t); %获取A(t)的采样点x=(1+a).*cos(2*pi*4000*t+2)+noisy/10; %获取x(t)的采样点%以t为横坐标画出x(t)的时域图型figure(1); subplot(2,1,1); plot(n,x); axis(0 140 -3 3);xlabel('采样点');ylabel('X(t)/V');title('窄带随机信号波形');grid on;%求X(t)的统计特性 并画出来disp('X(t)的均值为'); Ex=mean(x); disp(Ex);%求X(t)均值hold on; p

13、lot(n,Ex,'r.'); disp('X(t)的方差为');Dx=var(x); disp(Dx);%求x(t)方差hold on; plot(n,Dx,'g.');%画出X(t)的概率分布函数each=linspace(min(x),max(x),14); %将最大最小区间分成14等份,然后分别计算各个区间的个数nr=hist(x,each); %计算各个区间的个数nr=nr/length(x); %计算各个区间的个数归一化subplot(2,1,2); p=polyfit(each,nr,9); %画出概率分布直方图bar(each,n

14、r); %多项式拟合hold on; plot(each,nr,'g')eachi=-2:0.1:2;nri=polyval(p,eachi);plot(eachi,nri,'r')axis tight;title('X(t)概率密度分布');xlabel('X(t)');ylabel('P(x)');grid on;%对X(t)进行频谱分析Fx=fft(x,N); %对x(t)进行fft变换，在016000区间内得到2N-1个频率值magn=abs(Fx); %求x(t)幅值xangle=angle(Fx); %

15、求X(t)相位labelang=(0:length(x)-1)*16000/length(x); %在016000区间内求横坐标刻度figure(2); plot(labelang,magn*10); %在016000区间内做频谱和相位图axis(0 16000 -0.5 600); xlabel('频率/Hz');ylabel('幅值');title('X(t)频谱图');grid on;%求X(t)的自相关函数c,lags=xcorr(x,'coeff'); %求出自相关序列figure(3); subplot(2,1,1);

16、 plot(lags/fs,c); %在时域内画自相关函数axis tight; xlabel('T');ylabel('Rx(T)');title('X(t)的自相关函数');grid on;%求X(t)的功率谱密度long=length(c); Sx=fft(c,long); labelx=(0:long-1)*2*pi;plot_magn=10*log10(abs(Sx);subplot(2,1,2); plot(labelx,plot_magn); %画功率谱密度axis tight;xlabel('w');ylabel(

17、'Sx(w)');title('X(t)的功率谱密度');grid on; %窄带系统检测z1=2.*cos(2*pi*4000*t); z2=-2.*sin(2*pi*4000*t); Ac=z1.*x; %滤波后生成Ac(t)As=z2.*x; %滤波后生成As(t)y=Ac.*cos(2*pi*4000*t)-As.*sin(2*pi*4000*t);%滤波器设计f_p=1000;f_s=1600;R_p=1;R_s=35; %设定滤波器参数; 通、阻带截止频率,通、阻带衰减Ws=2*f_s/fs;Wp=2*f_p/fs; %频率归一化n,Wn=butto

18、rd(Wp,Ws,R_p,R_s); %采用巴特沃思滤波器b,a=butter(n,Wn); %求得滤波器传输函数的多项式系数figure(4);H,W=freqz(b,a); %求得滤波器传输函数的幅频特性subplot(2,1,1); plot(W*fs/(2*pi),abs(H); %在02pi区间内作幅度谱title('低通滤波器幅度谱'); grid on;subplot(2,1,2); plot(W*fs/(2*pi),angle(H); %在02pi区间内作相位谱title('低通滤波器相位谱'); grid on;%求Ac(t)滤波后的统计特性m

19、c=filter(b,a,Ac); %上支路通过滤波器 Ac(t)disp('Ac(t)的均值');Eh=mean(mc) %求Ac(t)的均值disp('Ac(t)的均方值是');E2h=mc*mc'/N %求Ac(t)的均方值disp('Ac(t)的方差');Dh=var(mc) %求Ac(t)的方差%画Ac(t)的时域波形figure(6); subplot(2,1,1); n=0:N-1; plot(n,mc); axis(0 300 -1 1);xlabel('采样点');ylabel('幅值')

20、;title('Ac(t)的时域波形');grid on;%画Ac(t)的频谱图yc=fft(mc,length(mc); %对Ac(t)进行fft变换longc=length(yc); %求傅里叶变换后的序列长度labelx=(0:longc-1)*16000/longc;magnl=abs(yc); %求Ac(t)的幅值subplot(2,1,2); plot(labelx,magnl); %画Ac(t)的频谱图axis tight; xlabel('频率(Hz)'); ylabel('幅值'); title('Ac(t)频谱图

21、9;); grid on;%求Ac(t)的自相关函数c1,lags1=xcorr(mc,'coeff'); %求出Ac(t)的自相关序列figure(7); subplot(2,1,1); plot(lags1/fs,c1); %在时域内画Ac(t)的自相关函数xlabel('T');ylabel('Rx(T)');axis tight;title('Ac(t)的自相关函数');grid on;%求Ac(t)的双边功率谱Sac=fft(c1,length(c1); %对Ac(t)的自相关函数进行傅里叶变换magnc=abs(Sac

22、); %求Ac(t)的双边功率谱幅值long=length(Sac); %求傅里叶变换后的序列长度labelc=(0:long-1)*16000/long;subplot(2,1,2); plot(labelc,10*log10(magnc); %画Ac(t)的自相关函数频谱 即为Ac(t)的双边功率谱xlabel('频率(Hz)');ylabel('功率谱(dbW)');axis tight;title('Ac(t)的双边功率谱');grid on;%求得As(t)的统计特性ms=filter(b,a,As); %对下支路信号进行滤波得As(t

23、)disp('As(t)的均值'); Eh=mean(ms) %求As(t)的均值disp('As(t)的均方值是'); E2h=ms*ms'/N %求As(t)的均方值disp('As(t)的方差'); Dh=var(ms) %求As(t)的方差 %作As(t)的时域波形figure(8);subplot(2,1,1); n=0:N-1;plot(n,ms); %画出As(t)的时域波形axis(0 300 -0.5 2); xlabel('采样点');ylabel('幅值');title('As

24、(t)的时域波形');grid on;%对As(t)进行FFT变换并做频谱图ys=fft(ms,length(ms); %对As(t)进行fft变换longs=length(ys); %求傅里叶变换后的序列长度labelx=(0:longs-1)*16000/longs;magn2=abs(ys); %求As(t)的幅值subplot(2,1,2); plot(labelx,magn2); %画出As(t)的频谱图axis tight;xlabel('频率(Hz)');ylabel('幅值');title('As(t)的频谱图');gri

25、d on;%求As(t)的自相关函数c2,lags2=xcorr(ms,'coeff'); %求出As(t)的自相关序列figure(9);subplot(2,1,1);plot(lags2/fs,c2); %画出As(t)自相关函数的时域波形xlabel('T');ylabel('Rx(T)');axis tight;title('As(t)的的自相关函数');grid on;%求As(t)的双边功率谱Sas=fft(c2,length(c2); %对As(t)的自相关函数进行傅里叶变换magnc=abs(Sac); %求As(

26、t)的双边功率谱幅值long=length(Sas); %求傅里叶变换后的序列长度labels=(0:long-1)*16000/long;subplot(2,1,2); plot(labelc,10*log10(magnc); %画As(t)的自相关函数频谱xlabel('频率(Hz)');ylabel('功率谱(dbW)');axis tight;title('As(t)的双边功率谱'); % 求y(t)的统计特性disp('输出信号Y(t)的均值');Eh=mean(y) %求输出信号Y(t)的均值 disp('输出

27、信号Y(t)的均方值');E2h=y*y'/N %求输出信号Y(t)的均方值disp('输出信号Y(t)的方差');Dh=var(y) %求输出信号Y(t)的方差 %作输出信号Y(t)的时域波形figure(10); subplot(2,1,1);n=0:N-1;plot(n,y); axis(0 150 -2 2);xlabel('采样点');ylabel('幅值');title('Y(t)的时域波形');grid on;%进行FFT变换并做频谱图yy=fft(y,length(y); %对相加后的信号进行fft

28、变换longy=length(yy); %Y(t)傅里叶变换后的序列长度labelx=(0:longy-1)*16000/longy;magn3=abs(yy); %求Y(t)的幅值subplot(2,1,2); plot(labelx,magn3); %做Y(t)的频谱图axis tight;xlabel('频率(Hz)');ylabel('幅值');title('Y(t)的频谱图');grid on;%求输出信号Y(t)的自相关函数c3,lags3=xcorr(y,'coeff'); %求出Y(t)的自相关序列figure(11); subplot(2,1,1); plot(lags3/fs,c3); %画Y(t)自相关函数的时域波形xlabel('T');ylabel('Rx(T)');axis tight;title('Y(t)的的自相关函数');grid on;%求输出信号Y(t)的双边功率谱Sy=fft(c3,length(c3); %对Y(t)的自相关函数进行傅里叶变换magny=abs(Sy); %求Y(t)双边功率谱幅值long=length(Sy);labely=(0:long-1)*16000/long;