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

以上兩幅圖為用Welch方法求出的平均周期圖，輸入數(shù)據(jù)為630點，每段30點，疊合15點，使用了漢明窗，不但功率譜變得平滑，而且分辨率也較好。4、Matlab實現(xiàn)程序(1)、輸入信號為兩個余弦函數(shù)之和的序列clear;h=1;forn=0:0.01:6.29x(h)=cos(50*n)+cos(80*n);h=h+1;end %輸入信號為兩個余弦函數(shù)之和的序列%用Welch平均估計試驗數(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)、輸入信號為余弦和噪聲之和的序列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; %輸入信號為余弦和噪聲之和的序列%用Welch平均估計試驗數(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方法編寫的函數(shù)functionp=welchlxb(x,M)N=length(x); %觀測數(shù)據(jù)的長度L=(2*N-M)/M; %分段后的段數(shù),M為每段長度d=hamming(M); %采用hamming窗U=0;forn=0:M-1U=(U+d(n+l)人2)/M;end%估計功率譜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個數(shù)據(jù)在原數(shù)據(jù)列中的位置forn=0:(M-1)psd(j)=psd(j)