实验： 典型时间序列的功率谱估计

1、实验二. 典型时间序列的功率谱估计一、实验内容与目标：了解有限长数据对谱估计的影响，重点研究周期图法和改进型周期图法的谱估计方法，并分析噪声对随机信号谱估计结果的影响。使学生了解随机信号的功率谱分析与主要估计方法。2、 实验任务（1） 对AR(2)模型产生的序列进行分析，并估计其数字特征。理论知识在随机信号分析课程的第五章时间序列模型中，我们对AR、MA及ARMA模型进行了分析。对于AR(2)模型：其解为：其中和是由初始条件确定的待定系数，；而根据其自相关函数，有：；功率谱为：。已知条件设有AR(2)模型为：即是高斯白噪声，均值为0，方差为4；任务一产生指定AR(P)模型的典型时间序列X(n)

2、；分别画出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,noise);subplot(211);plot(x);title(&#

3、39;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)为白噪声，将其看成输入，输出x(

4、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, 'ones

6、ided'); % 估计功率谱密度f=-fliplr(f) 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('实际功率谱曲线','理论功

7、率谱曲线')title('实际功率谱与理论功率谱曲线比较图（500点）')b=1;a=1 0.9 0.1;noise=normrnd(0,4,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:

8、end); % 对称谱Px/max(Px);subplot(222);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(222);plot(f,G,'-r');legend('实际功率谱曲线','理论功率谱曲线')title('实际功率谱与理论功率谱曲线比较图（2000点）')b=1;a=1 0.9 0.1;nois

9、e=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,Nfft,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=

10、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('实际功率谱与理论功率谱曲线比较图（10000点）')由图可知：（1）实际功率谱曲线相对于理论功率谱曲线有点误差，但其包络还是比较一致的。（2）有限长的数据和噪声都会影响谱估计结果，数据越长，功率谱估计越准确程序：(自相关)b=1

11、;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('自相关函数估计（2000点）')noise=normrnd(0,4,1,10000);x=filter(b,a,n

12、oise);R=xcorr(x,'coeff');subplot(213);plot(R)title('自相关函数估计（10000点）')（2） 不同信噪比下的功率谱估计与比较。任务对于给定的时间序列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)+cos(0.6*pi*n)+0.2*randn(1,N);Px1=10*lo

13、g10(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);subplot(212);

14、plot(f,Px2);xlabel('频率/Hz');ylabel('功率谱/dB');title('第二种白噪声功率谱');（3） 设计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);titl

15、e('经过滤波器前的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.945

16、8 p2=-0.8458，可得极点一个为正，一个为负，所以功率谱在w=0和出现峰值，w=0时极点0.9458离单位圆距离比w=时极点-0.8458离单位圆的距离小，则此时正极点作用大，w=0时|H(ejw)|最大。即|H(ejw)|在低频时分量较大，可以看作其具有低通性质。线性滤波器二程序及结果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);tit