實驗二典型時間序列的功率譜估計

1、實驗二. 典型時間序列的功率譜估計一、實驗內(nèi)容與目標：了解有限長數(shù)據(jù)對譜估計的影響，重點研究周期圖法和改進型周期圖法的譜估計方法，并分析噪聲對隨機信號譜估計結(jié)果的影響。使學生了解隨機信號的功率譜分析與主要估計方法。2、 實驗任務(wù)（1） 對AR(2)模型產(chǎn)生的序列進行分析，并估計其數(shù)字特征。理論知識在隨機信號分析課程的第五章時間序列模型中，我們對AR、MA及ARMA模型進行了分析。對于AR(2)模型：其解為：其中和是由初始條件確定的待定系數(shù)，；而根據(jù)其自相關(guān)函數(shù)，有：；功率譜為：。已知條件設(shè)有AR(2)模型為：即是高斯白噪聲，均值為0，方差為4；任務(wù)一產(chǎn)生指定AR(P)模型的典型時間序列X(n)

2、；分別畫出X(n)的500、2000點和10000個觀測點的波形，并估計他們的均值與方差；試根據(jù)上述理論知識推導并計算出該模型理論的均值和方差。1 / 13程序及結(jié)果程序：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);ti

3、tle('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，因為是一個線性系統(tǒng)，w(n)為白噪聲，將其看成輸

4、入，輸出x(n)也為0，所以方差就等于R(0) m(x)=0 2=r(0)=15.77分析隨著點數(shù)的增加（500,2000,10000），波形越來越密集，隨著點數(shù)的增加波形越能體現(xiàn)出時間序列模型的特點。任務(wù)二求出該AR(2)模型理論功率譜，在上等距選取K=500個點，畫出其理論功率譜曲線；程序及結(jié)果程序：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模型系統(tǒng)函數(shù)G=G/max(G);%歸一化f=w*fs/(2*pi);plot(f,G,'r')title('理

5、論功率譜密度曲線')xlabel('f')ylabel('幅值')任務(wù)三估計500、2000點和10000個觀測點的典型時間序列X(n)的自相關(guān)函數(shù)與功率譜，并與理論的功率譜曲線比較；程序、結(jié)果及分析程序：(功率譜)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, 

6、9;onesided'); % 估計功率譜密度f=-fliplr(f) f(1:end); %對稱頻率(反轉(zhuǎn))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、39;理論功率譜曲線')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); %對稱頻率(反轉(zhuǎn))Px=fliplr(Px)

8、 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+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.

9、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); %對稱頻率(反轉(zhuǎn))Px=fliplr(Px) Px(1:end); % 對稱譜Px/max(Px)subplot(223);plot(f,Px,'b');hold

10、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('實際功率譜與理論功率譜曲線比較圖（10000點）')由圖可知：（1）實際功率譜曲線相對于理論功率譜曲線有點誤差，但其包絡(luò)還是比較一致的。（2）有限長的數(shù)據(jù)和噪聲都會影響譜估計結(jié)果，數(shù)據(jù)越長，功率譜估計越準確程序：(自

11、相關(guān))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('自相關(guān)函數(shù)估計（500點）')noise=normrnd(0,4,1,2000);x=filter(b,a,noise);R=xcorr(x,'coeff');subplot(212);plot(R)title('自相關(guān)函數(shù)估計（2000點）')noise=normrnd(0,4,1,10000);x=filter

12、(b,a,noise);R=xcorr(x,'coeff');subplot(213);plot(R)title('自相關(guān)函數(shù)估計（10000點）')（2） 不同信噪比下的功率譜估計與比較。任務(wù)對于給定的時間序列Y(n)，設(shè)計2種不同白噪聲W(n)，分別畫出10000點的波形Y(n)+W(n)，估計功率譜，比較和分析估計結(jié)果；給定時間序列,該時間序列信號離散化采樣頻率為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

13、=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

14、(212);plot(f,Px2);xlabel('頻率/Hz');ylabel('功率譜/dB');title('第二種白噪聲功率譜');（3） 設(shè)計2種不同特性的線性濾波器，考察分布為N(0,4)的高斯白噪聲通過不同線性濾波器后的波形變化情況，畫出兩種濾波器下輸出前后的波形及系統(tǒng)的功率譜估計，并對波形變化情況加以討論。線性濾波器一程序及結(jié)果程序：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

15、);title('經(jīng)過濾波器前的x的波形');xlabel('n');ylabel('x');subplot(2,2,2);plot(y);title('經(jīng)過濾波器后的y的波形');xlabel('n');ylabel('y');subplot(2,2,3);plot(Psd);title('系統(tǒng)功率譜');xlabel('n');ylabel('Psd');分析用r,p,k=residue(1,0,1,a1,a2)函數(shù)求出H(z)的極點求的其極點p1

16、=0.9458 p2=-0.8458，可得極點一個為正，一個為負，所以功率譜在w=0和出現(xiàn)峰值，w=0時極點0.9458離單位圓距離比w=時極點-0.8458離單位圓的距離小，則此時正極點作用大，w=0時|H(ejw)|最大。即|H(ejw)|在低頻時分量較大，可以看作其具有低通性質(zhì)。線性濾波器二程序及結(jié)果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('經(jīng)過濾波器前的x的波形');xlabel('n');ylabel('x');subplot(2,2,2);plot(y);title('經(jīng)過濾波器后的y的波形');xlabel('n');ylabel('y');subplot(2,2,3);plot(Psd);title('系統(tǒng)功率譜');xlabel('n');ylabel(