試驗(yàn)六功率譜估計(jì)

實(shí)驗(yàn)六功率譜估計(jì)1、實(shí)驗(yàn)要求：理解功率譜估計(jì)的基本概念；掌握經(jīng)典功率譜估計(jì)方法——直接法和間接法；掌握改進(jìn)的經(jīng)典功率譜估計(jì)方法，例如Welch法。2、實(shí)驗(yàn)原理：功率譜估計(jì)就是基于有限的數(shù)據(jù)尋找信號(hào)、隨機(jī)過(guò)程或系統(tǒng)的頻率成分。它表示隨機(jī)信號(hào)頻率域的統(tǒng)計(jì)特性。隨機(jī)信號(hào)是無(wú)始無(wú)終具有無(wú)限能量的，所以其傅立葉變換并不存在，因?yàn)樗粷M(mǎn)足絕對(duì)可積的條件。因此需要研究其在頻率域上的功率分布情況，即功率譜密度或功率譜。經(jīng)典功率譜分析分為直接法和間接法。直接法，即周期圖法，是由傅立葉變換得到的：將隨機(jī)信號(hào)x(n)的N點(diǎn)樣本值x(n)看作能量有限信號(hào)，取其傅立葉N變換，得到X(ejw)；然后再取其幅值的平方，并除以N作為x(n)的真實(shí)功率譜NP(ejw)的估計(jì)，即1P(ejw)=IX(w)I2NN間接法，又稱(chēng)為自相關(guān)法或BT法,是由隨機(jī)信號(hào)N個(gè)觀(guān)察值x(0),...,x(N-1)，估計(jì)出自相關(guān)函數(shù)R(m)，然后再求R(m)的傅立葉變換作為功率譜的估計(jì)：NNS(ejw)=丈R(m)e-jWm IMI<N—1Nm=—M直接法和間接法的方差性能很差，而且當(dāng)數(shù)據(jù)長(zhǎng)度太大時(shí)，譜曲線(xiàn)起伏加??；若數(shù)據(jù)長(zhǎng)度太小，則譜的分辨率又不好，所以需要改進(jìn)。改進(jìn)的直接譜估計(jì)方法由Bartlett法和Welch法。Bartlett法將采樣數(shù)據(jù)x(n)分成L段,每段的長(zhǎng)度都是M,即N=LM,對(duì)每N段數(shù)據(jù)加矩形窗，再計(jì)算其各自的功率譜Pi(w)，把Pi(w)對(duì)應(yīng)相加，再取平PER PER均,得到平均周期圖P_(w)。Welch法是對(duì)Bartlett法的改進(jìn)：一,在對(duì)x(n)分PER N段時(shí)，可允許每段數(shù)據(jù)有部分重疊；二，每段數(shù)據(jù)窗口可以不是矩形窗口，例如使用漢寧窗或哈明窗，記為d(n)。然后按Bartlett法求每一段的功率譜，記為Pi(w)=Pi(w)=PER1MUI￡xi(n)d(n)e-jwnI2,Nn=0其中U=M遲d2(n。平均后的功率譜為:n=0(段數(shù)L(段數(shù)L=N-M/2)M/2P(w)=1￡P(guān)i(w)PER L PERi=13、實(shí)驗(yàn)內(nèi)容與結(jié)果(1)、若輸入信號(hào)為n=0?2兀的x=cos(50n)+cos(80n)余弦序列，用改進(jìn)的經(jīng)典功率譜估計(jì)方法Welch算法估計(jì)其功率譜，圖形如下：(2)、若輸入信號(hào)為n=0?2兀的cos(50n)+cos(80n)與均值為零、方差為1白噪聲組成的序列，用改進(jìn)的經(jīng)典功率譜估計(jì)方法Welch算法估計(jì)其功率譜，圖形如下：

以上兩幅圖為用Welch方法求出的平均周期圖，輸入數(shù)據(jù)為630點(diǎn)，每段30點(diǎn)，疊合15點(diǎn)，使用了漢明窗，不但功率譜變得平滑，而且分辨率也較好。4、Matlab實(shí)現(xiàn)程序(1)、輸入信號(hào)為兩個(gè)余弦函數(shù)之和的序列clear;h=1;forn=0:0.01:6.29x(h)=cos(50*n)+cos(80*n);h=h+1;end %輸入信號(hào)為兩個(gè)余弦函數(shù)之和的序列%用Welch平均估計(jì)試驗(yàn)數(shù)據(jù)的功率譜；p=welchlxb(x,30);mmax=max(p);p=p/mmax;w=0:1/629:1;plot(w,p);gridon;axis([0,0.5,0,1])xlabel('w/(2*pi)');ylabel('p(w)');(2)、輸入信號(hào)為余弦和噪聲之和的序列clear;h=1;forn=0:0.01:6.29x(h)=cos(50*n)+cos(80*n);h=h+1;endu=randn(1,630);x=x+u; %輸入信號(hào)為余弦和噪聲之和的序列%用Welch平均估計(jì)試驗(yàn)數(shù)據(jù)的功率譜；p=welchlxb(x,30);mmax=max(p);p=p/mmax;w=0:1/629:1;plot(w,p);gridon;axis([0,0.5,0,1])xlabel('w/(2*pi)');ylabel('p(w)');(3)、用Welch方法編寫(xiě)的函數(shù)functionp=welchlxb(x,M)N=length(x); %觀(guān)測(cè)數(shù)據(jù)的長(zhǎng)度L=(2*N-M)/M; %分段后的段數(shù),M為每段長(zhǎng)度d=hamming(M); %采用hamming窗U=0;forn=0:M-1U=(U+d(n+l)人2)/M;end%估計(jì)功率譜p=zeros(630,1); %在[0,2*pi]內(nèi)按100Hz取樣psd=zeros(L,l);k=1;forw=0:0.01:6.29forj=1:Listart=(j-1)*M/2; %第j段數(shù)據(jù)的第1個(gè)數(shù)據(jù)在原數(shù)據(jù)列中的位置forn=0:(M-1)psd(j)=psd(j)