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

实验二.典型时间序列的功率谱估计一、实验内容与目标：了解有限长数据对谱估计的影响，重点研究周期图法和改进型周期图法的谱估计方法，并分析噪声对随机信号谱估计结果的影响。使学生了解随机信号的功率谱分析与主要估计方法。二、实验任务（一）对AR(2)模型产生的序列进行分析，并估计其数字特征。理论知识在《随机信号分析》课程的第五章时间序列模型中，我们对AR、MA及ARMA模型进行了分析。对于AR(2)模型：)()2()1()(21nnxanxanx其解为：02112112211)()(kkknnknzzzzzAzAnx其中1A和2A是由初始条件确定的待定系数，22112,1421aaaz；而根据其自相关函数，有：])1)[(1()1()0(122222aaaaRnx；功率谱为：,1)(22212jjnxeaeaS。已知条件设有AR(2)模型为：)()2(1.0)1(9.0)(nnxnxnx即)()2(1.0)1(9.0)(nnxnxnx)(n是高斯白噪声，均值为0，方差为4；任务一产生指定AR(P)模型的典型时间序列X(n)；分别画出X(n)的500、2000点和10000个观测点的波形，并估计他们的均值与方差；试根据上述理论知识推导并计算出该模型理论的均值和方差。程序及结果程序：b=1;a=[10.90.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('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.0404sigma2=14.3090m=-0.0486sigma2=14.5322m=0.0058sigma2=15.0434理论的均值和方差均值为0，因为是一个线性系统，w(n)为白噪声，将其看成输入，输出x(n)也为0，所以方差就等于R(0)m(x)=0δ2=r(0)=15.77分析随着点数的增加（500,2000,10000），波形越来越密集，随着点数的增加波形越能体现出时间序列模型的特点。任务二求出该AR(2)模型理论功率谱，在0上等距选取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('理论功率谱密度曲线')xlabel('f')ylabel('幅值')任务三估计500、2000点和10000个观测点的典型时间序列X(n)的自相关函数与功率谱，并与理论的功率谱曲线比较；程序、结果及分析程序：(功率谱)b=1;a=[10.90.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)f(1:end)];%对称频率(反转)Px=[fliplr(Px)Px(1:end)];%对称谱Px/max(Px)subplot(221);plot(f,Px,'b');holdon;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=[10.90.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:end)];%对称谱Px/max(Px);subplot(222);plot(f,Px,'b');holdon;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=[10.90.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,Nfft,fs,'onesided');%估计功率谱密度f=[-fliplr(f)f(1:end)];%对称频率(反转)Px=[fliplr(Px)Px(1:end)];%对称谱Px/max(Px)subplot(223);plot(f,Px,'b');holdon;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('实际功率谱与理论功率谱曲线比较图（10000点）')由图可知：（1）实际功率谱曲线相对于理论功率谱曲线有点误差，但其包络还是比较一致的。（2）有限长的数据和噪声都会影响谱估计结果，数据越长，功率谱估计越准确程序：(自相关)b=1;a=[10.90.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,noise);R=xcorr(x,'coeff');subplot(213);plot(R)title('自相关函数估计（10000点）')（二）不同信噪比下的功率谱估计与比较。任务对于给定的时间序列Y(n)，设计2种不同白噪声W(n)，分别画出10000点的波形Y(n)+W(n)，估计功率谱，比较和分析估计结果；给定时间序列)6.0cos()3.0sin()(nnnY,该时间序列信号离散化采样频率为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*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);subplot(212);plot(f,Px2);xlabel('频率/Hz');ylabel('功率谱/dB');title('第二种白噪声功率谱');（三）设计2种不同特性的线性滤波器，考察分布为N(0,4)的高斯白噪声通过不同线性滤波器后的波形变化情况，画出两种滤波器下输出前后的波形及系统的功率谱估计，并对波形变化情况加以讨论。线性滤波器一)()2(8.0-)1(1.0-)(nnxnxnx程序及结果程序：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.9458p2=-0.8458，可得极点一个为正，一个为负，所以功率谱在w=0和π出现峰值，w=0时极点0.9458离单位圆距离比w=π时极点-0.8458离单位圆的距离小，则此时正极点作用大，w=0时|H(ejw)|最大。即|H(ejw)|在低频时分量较大，可以看作其具有低通性质。线性滤波器二)()2(8.0-)1(1.0)(nnxnxnx程序及结果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.9458p2=0.8458，可得极点一个为正，一个为负，所以功率谱在w=0和π出现峰值，w=0时极点-0.945