窄带随机信号分析与仿真

1、随机信号实验实验报告题 目 窄带随机信号的仿真与分析实验十三 窄带随机信号的仿真与分析实验目的利用计算机仿真窄带随机信号，考察其数字特征，以此加深对窄带随机信号的特点及分析方法。熟悉常用的信号处理仿真软件平台：matlab、c/c+。 实验原理如果一个随机过程的功率谱密度，只分布在高频载波附近的一个窄频率范围内，在此范围之外全为零，则称之为窄带过程。窄带过程是在信息传输系统，特别是接收机中经常遇到的随机信号。当窄带系统(接收机)的输入噪声(如热噪声)的功率谱分布在足够宽的频带(相对于接收机带宽)上时，系统的输出即为窄带过程。窄带信号的确切定义如下：一个实平稳随机x(t)，如果它的功率谱密度sx

2、()具有下述性质 而且带宽满足，则称此过程为窄带平稳随机过程。窄带平稳随机过程的功率谱密度函数如图24所示。 图24 窄带随机过程的功率谱密度从示波器观看窄带随机过程的一个样本函数，可看到如图25所示的波形，从这个波形可以看出，窄带随机过程可以表示成具有角频率以及慢变幅度与相位的正弦振荡，这就是说可以把它写成 z(t)=b(t)cost + (t) 式中，b(t)是随机过程的慢变幅度，(t)是过程的慢变相位。称之为准正弦振荡，这就是窄带过程的数学模型。 图25 窄带过程的某个样本函数 实验任务与要求 用matlab或c/c+语言编写和仿真程序。 产生满足下列条件的窄带随机信号： 其中a(t)包

3、络频率为1khz,幅值为1v。载波频率为：4khz，幅值为1v，是一个固定相位，n(t)为高斯白噪声，采样频率设为16khz。实际上，这是一个带有载波的双边带调制信号，如下图：计算窄带随机信号的均值、均方值、方差、概率密度、频谱及功率谱密度、相关函数，用图示法来表示。 窄带系统检测框图如图26所示： 图26 窄带系统检测 x(t)：窄带随机信号、：载波 低通滤波器设计：低通滤波器技术要求：通带截止频率1khz阻带截止频率2khz。过渡带:1khz阻带衰减：35db通带衰减：1db采样频率：44.1khz 计算a点、b点、ac(t)、as(t)、y(t)的均值、均方值、方差、频谱及功率谱密度、相

4、关函数，用图示法来表示。下图为ac(t)的时域、频域波形： 注意观察、在系统中起到什么作用，调制信号的频谱结构有没有变化。 y(t)还原了什么信号。 按要求写实验报告。实验步骤：1. 根据实验要求得到x(t)的表达式：x=(1+a).*cos(2*pi*4000*t+2)+noisy/10;其中：noisy为高斯白噪声，由wgn函数生成a=cos(2*pi*1000*t)均值：ex=mean(x);方差：dx=var(x);可得：x(t)的均值为3.1440e-004x(t)的方差为0.7667x(t)的时域波形，均值还有方差如图01所示：图01-x(t)的波形（蓝色），方差（上方的直线），均

5、值（下方的直线）2. 求x(t)的概率密度方法是将最大最小区间分成14等份,然后分别计算各个区间的个数，如图02中柱形条所示，利用曲线拟合，得到合适的概率密度函数。为了得到光滑的曲线，我们利用了多项式拟合，经过测试，9次拟合曲线比较符合我们的要求，获得的曲线如图02中曲线所示：图02-x(t)概率密度分布3. 对x(t)进行频谱分析在matlab中，利用fft函数可以很方便得求得x(t)的频谱，然后用abs和angle函数求得幅值和相位，画出图像如图03所示：图03-x(t)的频谱分析图4. 求x(t)的自相关函数我们的方法是用xcorr函数求出自相关序列，得到自相关函数的图像如图04所示：图

6、04-x(t)的自相关函数图像5. 求x(t)功率谱密度有了自相关函数，对它进行fft变换，可得x(t)的功率谱密度，如图05所示：图05-x(t)的功率谱密度6. 建立滤波器下面我们建立一个巴特沃思滤波器对产生的x(t)进行检测。滤波器的幅度谱和相位谱如图06，图07所示：图06-巴特沃思滤波器幅度谱图07-巴特沃思滤波器相位谱7. 求ac的统计特性ac为x(t)*2*cos(2*4000*t)通过步骤6中建立的低通滤波器的信号，其各项统计特性如下：ac输出信号均值eh=-0.0024ac的均方值8.8873e-004ac的方差8.8383e-004各图像如下列图所示：图08-ac(t)的时

7、域图图09-ac(t)的频谱图图10-ac(t)的自相关函数图11-ac(t)的功率谱密度8. 求as的统计特性as为x(t)*2*sin(2*4000*t)通过步骤6中建立的低通滤波器的信号，其各项统计特性如下：as输出信号均值eh= -0.0018as的均方值0.0039as的方差0.0039各图像如下列图所示：图12-as(t)的时域图图13-as(t)的频谱图图14-as(t)的自相关函数图15-as(t)的功率谱密度9. 求输出y(t)的统计特性最后，我们来求y(t)的统计特性,y(t)的表达式如下：y(t)=ac(t)*cos(2*4000*t)+as(t)*sin(2*4000*

8、t)所用的方法和前述一样，其统计特性如下y(t)的均值：0.1325y(t)的均方值：0.0720y(t)的方差：0.0544各图像如下列图所示：图16-y(t)的时域图图17-y(t)的频谱图图18-y(t)的自相关函数图19-y(t)的功率谱密度附件1：心得体会终于在小组的努力下完成了这个我们大学真正独立完成的随机信号实验，辛勤的汗水换回那一张一张图片是那么的让人赏心悦目。虽然这不是我们的第一次合作，但合作的重要性又一次让我们明白其价值的不菲。对于我们这些只知道念书读书的大学生来说，这次实验的意义远不止完成这个实验，更重要是让我们明白了知识与实践的桥梁不是我们一直认为的那么简单，只有在那一

9、次次的调试我们自己的程序，才会明白一切的一切都是需要我们去重新学习。两年的时光足以忘记很多东西，大一学习的matlab再次回来自己的手边是那么陌生，自以为学会编程的我们面对那一个一个的错误才会明白知识的无止境。实验中存在的问题：1 对于matlab太陌生，没有系统的学习过那些真正用于信号分析的那些工具箱。2 对于随机信号的知识书本上的那些知识，没有学宽，学出去。3 对某些小的知识点了解不够问题解决方法：1. 合理的分工，大大提高了工作效率。2. 合理的利用可以使用的一切资源，对于遇到的问题集体讨论并上网寻求最佳的答案。3. 对于无法解决的问题，进一步咨询了老师，争取做到答案的完美。在三人的共同

10、努力下，我们完成了这次实验，希望老师验收并加以指正。附件2：程序% 说明：此程序是随机实验十三 窄带随机信号的仿真与分析的程序 % 班级：010852班% 作者：邓晓光 刘春燕 谭正中 %- 生成x(t)并求得x(t)的统计特性 -%clcfs=16000; %设定采样频率n=1300; %取的样本点数n=0:n-1; %取样点序列length(n)t=n/fs; %获得以1/16000为时间间隔采样序列noisy=wgn(1,n,0); %产生高斯白噪声a=cos(2*pi*1000*t); %获取a(t)的采样点x=(1+a).*cos(2*pi*4000*t+2)+noisy/10; %

11、获取x(t)的采样点figure(1); %画图subplot(2,1,1);plot(n,x);%以t为横坐标画出x(t)的时域图型axis(0 140 -3 3);%axis tight; xlabel(采样点);ylabel(x(t)/v);title(图1-窄带随机信号波形);grid on;%求x(t)均值disp(x(t)的均值为);ex=mean(x);hold on;plot(n,ex,r.);disp(ex);%求x(t)方差disp(x(t)的方差为)dx=var(x);hold on;plot(n,dx,g.)disp(dx);each=linspace(min(x),m

12、ax(x),14);%将最大最小区间分成19等份,然后分别计算各个区间的个数nr=hist(x,each);%计算各个区间的个数nr=nr/length(x);%计算各个区间的个数归一化subplot(2,1,2);%多项式拟合bar(each,nr);hold on;plot(each,nr,g)hold on;p=polyfit(each,nr,9);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)

13、进行频谱分析fx=fft(x,n); %对x(t)进行fft变换,在016000区间内得到2n-1个频率值magn=abs(fx); %求x(t)幅值xangle=angle(fx);%求x(t)相位labelang=(0:length(fx)-1)*16000/length(x);%在016000区间内求横坐标刻度figure(2); %subplot(2,2,3);plot(labelang,magn*10,labelang,xangle);%在016000区间内做频谱和相位图%plot(labelang,magn);%在02区间内做频谱图axis(0 16000 -0.5 600) ;

14、xlabel(频率/hz);ylabel(幅值 );title(x(t)频谱图);grid on;%求x(t)的自相关函数c,lags=xcorr(x,coeff);%求出自相关序列figure(3); subplot(2,1,1);plot(lags/fs,c); %在时域内画自相关函数%axis(-0.003 0.003 -0.5 1);axis tight; xlabel();ylabel( rx();title(x(t)的自相关函数);grid on;%求x(t)的功率谱密度long=length(c);sx=fft(c,long); %求出自相关序列 index=0:round(n/

15、2-1);labelx=(0:long-1)*2*pi/long;plot_magn=10*log10(abs(sx);subplot(2,1,2);plot(labelx,plot_magn); %画功率谱密度axis tight; xlabel(w);ylabel( s(w);title(x(t)的功率谱密度);grid on;%窄带系统检测z1=x.*cos(2*pi*4000*t);z2=x.*sin(2*pi*4000*t);ac=z1.*x;as=z2.*x;y=ac.*cos(2*pi*4000*t)+as.*sin(2*pi*4000*t);f_p=1000;f_s=1600;

16、r_p=1;r_s=35;%设定滤波器参数:通、阻带截止频率，通、阻带衰减ws=2*f_s/fs;wp=2*f_p/fs; %频率归一化n,wn=buttord(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);%在 02 区间内作幅度谱title(低通滤波器幅度谱)grid on;subplot(2,1,2);plot(w*fs/(2*pi),angle(h);%在 02

17、 区间内作相位谱title(低通滤波器相位谱);grid on;%- 生成ac(t)并求得ac(t)的统计特性 -%信号通过滤波器mc=filter(b,a,ac);%对信号进行滤波%求输出信号均值disp(ac输出信号均值);eh=mean(mc)%求输出信号的均方值disp(ac的均方值是);e2h=mc*mc/n%求输出信号方差disp(ac输出信号方差)dh=var(mc)%作输出信号时域波形figure(6)subplot(2,1,1);n=0:n-1;plot(n,mc);%作信号时域波形axis(0 300 -0.1 0.1);% axis tight; xlabel(采样点 )

18、;ylabel(幅值) title(ac时域谱);grid on;%进行 fft 变换并做频谱图yc=fft(mc,length(mc);%对滤波后信号进行 fft 变换longc=length(yc);%求傅里叶变换后的序列长度labelx=(0:longc-1)*16000/longc;magn1=abs(yc);% 求滤波后信号幅值subplot(2,1,2);plot(labelx,magn1);%做频谱图axis tight; xlabel(频率(hz); ylabel(幅值 ); title(ac频谱图);grid on;%求输出信号的自相关函数c1,lags1=xcorr(mc,

19、coeff);%求出自相关序列figure(7);subplot(2,1,1);plot(lags1/fs,c1);%在时域内画自相关函数xlabel(); ylabel(rx(); axis tight;title(ac的自相关函数);grid on;%求输出信号双边功率谱sac=fft(c1,length(c1);%对自相关函数进行傅里叶变换magnc=abs(sac);%求双边功率谱幅值subplot(2,1,2);long=length(sac);%求傅里叶变换后的序列长度labelc=(0:long-1)*16000/long;plot(labelc,10*log10(magnc);

20、%画自相关函数谱xlabel(频率(hz);ylabel(功率谱(dbw);axis tight;title(ac双边功率谱);grid on;%- 生成as(t)并求得as(t)的统计特性 -%信号通过滤波器ms=filter(b,a,as);%对信号进行滤波%求输出信号均值disp(as输出信号均值);eh=mean(ms)%求输出信号的均方值disp(as的均方值是);e2h=ms*ms/n%求输出信号方差disp(as输出信号方差)dh=var(ms)%作输出信号时域波形figure(8)subplot(2,1,1);n=0:n-1;plot(n,ms);%作信号时域波形axis(0

21、300 -0.1 0.1);% axis tight; xlabel(采样点 );ylabel(幅值) title(as时域谱);grid on;%进行 fft 变换并做频谱图ys=fft(ms,length(ms);%对滤波后信号进行 fft 变换longs=length(ys);%求傅里叶变换后的序列长度labelx=(0:longs-1)*16000/longs;magn2=abs(ys);% 求滤波后信号幅值subplot(2,1,2);plot(labelx,magn2);%做频谱图axis tight; xlabel(频率(hz); ylabel(幅值 ); title(as频谱图

22、);grid on;%求输出信号的自相关函数c2,lags2=xcorr(ms,coeff);%求出自相关序列figure(9);subplot(2,1,1);plot(lags2/fs,c2);%在时域内画自相关函数xlabel(); ylabel(rx(); axis tight;title(as的自相关函数);grid on;%求输出信号双边功率谱sas=fft(c2,length(c2);%对自相关函数进行傅里叶变换magnc=abs(sac);%求双边功率谱幅值subplot(2,1,2);long=length(sas);%求傅里叶变换后的序列长度labels=(0:long-1)*16000/long;plot(labelc,10*log10(magnc);%画自相关函数谱xlabel(频率(hz);ylabel(功率谱(dbw);axis tight;title(as双边功率谱);grid on;%- 生成y(t)并求得y(t)的统计特性 -%求输出信号均值disp(y输出信号均值);eh=mean(y)%求输出信号的均方值disp(y的均方值是);e2h=y*y/n%求输出信号方差disp(y输出信号方差)