实验二典型时间序列的功率谱估计.doc

《实验二典型时间序列的功率谱估计.doc》由会员分享，可在线阅读，更多相关《实验二典型时间序列的功率谱估计.doc（13页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、【精品文档】如有侵权，请联系网站删除，仅供学习与交流实验二典型时间序列的功率谱估计.精品文档.实验二. 典型时间序列的功率谱估计一、实验内容与目标：了解有限长数据对谱估计的影响，重点研究周期图法和改进型周期图法的谱估计方法，并分析噪声对随机信号谱估计结果的影响。使学生了解随机信号的功率谱分析与主要估计方法。二、 实验任务（一） 对AR(2)模型产生的序列进行分析，并估计其数字特征。理论知识在随机信号分析课程的第五章时间序列模型中，我们对AR、MA及ARMA模型进行了分析。对于AR(2)模型：其解为：其中和是由初始条件确定的待定系数，；而根据其自相关函数，有：；功率谱为：。已知条件设有AR(2)

2、模型为：即是高斯白噪声，均值为0，方差为4；任务一产生指定AR(P)模型的典型时间序列X(n)；分别画出X(n)的500、2000点和10000个观测点的波形，并估计他们的均值与方差；试根据上述理论知识推导并计算出该模型理论的均值和方差。程序及结果程序：b=1;a=1 0.9 0.1;noise=normrnd(0,2,1,500);x=filter(b,a,noise); subplot(211);plot(x);title(AR(2)随机序列500点);m=mean(x);sigma2=var(x);msigma2noise=normrnd(0,2,1,2000);x=filter(b,a

3、,noise);subplot(211);plot(x);title(AR(2)随机序列2000点);m=mean(x);sigma2=var(x);msigma2noise=normrnd(0,2,1,10000);x=filter(b,a,noise);subplot(211);plot(x);title(AR(2)随机序列10000点);m=mean(x);sigma2=var(x);msigma2m=0.0404 sigma2=14.3090m=-0.0486 sigma2=14.5322m=0.0058 sigma2=15.0434理论的均值和方差均值为0，因为是一个线性系统，w(n

4、)为白噪声，将其看成输入，输出x(n)也为0，所以方差就等于R(0) m(x)=0 2=r(0)=15.77分析随着点数的增加（500,2000,10000），波形越来越密集，随着点数的增加波形越能体现出时间序列模型的特点。任务二求出该AR(2)模型理论功率谱，在上等距选取K=500个点，画出其理论功率谱曲线；程序及结果程序：fs=1000;%采样频率w=0:pi/2000:pi;G=4*(abs(1./(1+0.9*exp(-j*w)+0.1*exp(-2*j*w).2);%AR模型系统函数G=G/max(G);%归一化f=w*fs/(2*pi);plot(f,G,r)title(理论功率谱

5、密度曲线)xlabel(f)ylabel(幅值)任务三估计500、2000点和10000个观测点的典型时间序列X(n)的自相关函数与功率谱，并与理论的功率谱曲线比较；程序、结果及分析程序：(功率谱)b=1;a=1 0.9 0.1;noise=normrnd(0,4,1,500);x=filter(b,a,noise);window=hamming(20); % 采用hanmming 窗，长度为20noverlap=10;Nfft=512;fs=1000;Px,f=pwelch(x,window,noverlap,Nfft,fs, onesided); % 估计功率谱密度f=-fliplr(f)

6、 f(1:end); %对称频率(反转)Px=fliplr(Px) Px(1:end); % 对称谱Px/max(Px)subplot(221);plot(f,Px,b);hold on;fs=1000;w=-pi:1/fs:pi;G=4*(abs(1./(1+0.9*exp(-j*w)+0.1*exp(-2*j*w).2);G=G/max(G);f=w*fs/(2*pi);subplot(221);plot(f,G,-r)legend(实际功率谱曲线,理论功率谱曲线)title(实际功率谱与理论功率谱曲线比较图（500点）)b=1;a=1 0.9 0.1;noise=normrnd(0,4,

7、1,2000);x=filter(b,a,noise);window=hamming(20); % 采用hanmming 窗，长度为20noverlap=10;Nfft=512;fs=1000;Px,f=pwelch(x,window,noverlap,Nfft,fs, onesided); % 估计功率谱密度f=-fliplr(f) f(1:end); %对称频率(反转)Px=fliplr(Px) Px(1:end); % 对称谱Px/max(Px);subplot(222);plot(f,Px,b);hold on;fs=1000;w=-pi:1/fs:pi;G=4*(abs(1./(1+

8、0.9*exp(-j*w)+0.1*exp(-2*j*w).2);G=G/max(G);f=w*fs/(2*pi);subplot(222);plot(f,G,-r);legend(实际功率谱曲线,理论功率谱曲线)title(实际功率谱与理论功率谱曲线比较图（2000点）)b=1;a=1 0.9 0.1;noise=normrnd(0,4,1,10000);x=filter(b,a,noise);window=hamming(20); % 采用hanmming 窗，长度为20noverlap=10;Nfft=512;fs=1000;Px,f=pwelch(x,window,noverlap,N

9、fft,fs, onesided); % 估计功率谱密度f=-fliplr(f) f(1:end); %对称频率(反转)Px=fliplr(Px) Px(1:end); % 对称谱Px/max(Px)subplot(223);plot(f,Px,b);hold on;fs=1000;w=-pi:1/fs:pi;G=4*(abs(1./(1+0.9*exp(-j*w)+0.1*exp(-2*j*w).2);G=G/max(G);f=w*fs/(2*pi);subplot(223);plot(f,G,-r);legend(实际功率谱曲线,理论功率谱曲线)title(实际功率谱与理论功率谱曲线比较图

10、（10000点）)由图可知：（1）实际功率谱曲线相对于理论功率谱曲线有点误差，但其包络还是比较一致的。（2）有限长的数据和噪声都会影响谱估计结果，数据越长，功率谱估计越准确程序：(自相关)b=1;a=1 0.9 0.1;noise=normrnd(0,4,1,500);x=filter(b,a,noise);R=xcorr(x,coeff);subplot(211);plot(R)title(自相关函数估计（500点）)noise=normrnd(0,4,1,2000);x=filter(b,a,noise);R=xcorr(x,coeff);subplot(212);plot(R)title

11、(自相关函数估计（2000点）)noise=normrnd(0,4,1,10000);x=filter(b,a,noise);R=xcorr(x,coeff);subplot(213);plot(R)title(自相关函数估计（10000点）)（二） 不同信噪比下的功率谱估计与比较。任务对于给定的时间序列Y(n)，设计2种不同白噪声W(n)，分别画出10000点的波形Y(n)+W(n)，估计功率谱，比较和分析估计结果；给定时间序列,该时间序列信号离散化采样频率为1KHz。程序：Fs=1000;N=1024;Nfft=10000;n=0:N-1;%第一种白噪声yn=sin(0.3*pi*n)+c

12、os(0.6*pi*n)+0.2*randn(1,N);Px1=10*log10(abs(fft(yn,Nfft).2)/N);f=(0:length(Px1)-1)*Fs/length(Px1);subplot(211);plot(f,Px1);xlabel(频率/Hz);ylabel(功率谱/dB);title(第一种白噪声功率谱);%第二种白噪声yn=sin(0.3*pi*n)+cos(0.6*pi*n)+0.8*randn(1,N);Px2=10*log10(abs(fft(yn,Nfft).2)/N);f=(0:length(Px2)-1)*Fs/length(Px2);subplo

13、t(212);plot(f,Px2);xlabel(频率/Hz);ylabel(功率谱/dB);title(第二种白噪声功率谱);（三） 设计2种不同特性的线性滤波器，考察分布为N(0,4)的高斯白噪声通过不同线性滤波器后的波形变化情况，画出两种滤波器下输出前后的波形及系统的功率谱估计，并对波形变化情况加以讨论。线性滤波器一程序及结果程序：N=1000;x=normrnd(0,1,N,1);b=1;a=1,-0.1,-0.8;y=filter(b,a,x);Psd=abs(fft(y).2/N;subplot(2,2,1);plot(x);title(经过滤波器前的x的波形);xlabel(n

14、);ylabel(x);subplot(2,2,2);plot(y);title(经过滤波器后的y的波形);xlabel(n);ylabel(y);subplot(2,2,3);plot(Psd);title(系统功率谱);xlabel(n);ylabel(Psd);分析用r,p,k=residue(1,0,1,a1,a2)函数求出H(z)的极点求的其极点p1=0.9458 p2=-0.8458，可得极点一个为正，一个为负，所以功率谱在w=0和出现峰值，w=0时极点0.9458离单位圆距离比w=时极点-0.8458离单位圆的距离小，则此时正极点作用大，w=0时|H(ejw)|最大。即|H(ej

15、w)|在低频时分量较大，可以看作其具有低通性质。线性滤波器二程序及结果N=1000;x=normrnd(0,1,N,1);b=1;a=1,0.1,-0.8;y=filter(b,a,x);Psd=abs(fft(y).2/N;subplot(2,2,1);plot(x);title(经过滤波器前的x的波形);xlabel(n);ylabel(x);subplot(2,2,2);plot(y);title(经过滤波器后的y的波形);xlabel(n);ylabel(y);subplot(2,2,3);plot(Psd);title(系统功率谱);xlabel(n);ylabel(Psd);分析用r,p,k=residue(1,0,1,a1,a2)函数求出H(z)的极点求的其极点p1=-0.9458 p2=0.8458，可得极点一个为正，一个为负，所以功率谱在w=0和出现峰值，w=0时极点-0.9458离单位圆距离比w=时极点0.8458离单位圆的距离小，则此时负极点作用大，w=时|H(ejw)|最大。即|H(ejw)|在高频时分量较大，可以看作其具有高通性质。三、 体会与收获