【正文】 figure(2) plot(f,plot_Pxx1)。 南京郵電大學(xué) 通達(dá)學(xué)院 2020 屆本科生畢業(yè)設(shè)計(jì)（論文） figure(1) plot(f,plot_Pxx)。 plot_Pxx1=10*log10(Pxx1)。 [Pxx2,f]=pwelch(xn,window2,noverlap,nfft,Fs,range)。 %頻率間隔為 [0 Fs/2]，只計(jì)算一半的頻率 [Pxx,f]=pwelch(xn,window,noverlap,nfft,Fs,range)。 %數(shù)據(jù)無重疊 range=39。%漢明窗 window2=blackman(100)。 window=boxcar(100)。 xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n))。 Fs=1000。 南京郵電大學(xué) 通達(dá)學(xué)院 2020 屆本科生畢業(yè)設(shè)計(jì)（論文） 0 50 100 150 200 250 300 350 400 450 50005101520253035自相關(guān)仿真圖 自相關(guān)法仿真圖 以下是利用 matlab 對(duì)加窗周期圖法進(jìn)行仿真，并得到仿真圖。plot_Pxx=10*log10(Pxx(index+1))。%求絕對(duì)值 ； index=0:round(nfft/21)。%計(jì)算序列的自相關(guān)函數(shù) CXk=fft(cxn,nfft) %對(duì)自相關(guān)函數(shù)進(jìn)行傅里葉變換 。unbiased39。 nfft=1024。 n=0:1/Fs:1。程序如下： clear。 figure(2) plot(k,[plot_Pxx plot_Pxxplot_Pxxc plot_Pxx+plot_Pxxc])。 南京郵電大學(xué) 通達(dá)學(xué)院 2020 屆本科生畢業(yè)設(shè)計(jì)（論文） figure(1) plot(k,plot_Pxx)。 plot_Pxx=10*log10(Pxx(index+1))。 index=0:round(nfft/21)。 %數(shù)據(jù)無重疊 p=。window=boxcar(length(n))。 xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n))。 Fs=1000。 0 100 200 300 400 500 600 700 800 900 1000 6 0 5 0 4 0 3 0 2 0 1 0010 f Pxx周期圖 周期圖法仿真圖 以下是利用 matlab 對(duì)平均周期圖法進(jìn)行仿真，并得到仿真圖。%加矩形窗 [Pxx,f]=periodgram(xn,window,nfft,Fs)。 nfft=1024。%采樣頻率 n=0:1/Fs:1。程序如下： clear?，F(xiàn)代譜估計(jì)主要是針對(duì)經(jīng)典譜估計(jì)的分辨率低和方差性能不好的問題而提出的，其內(nèi)容極其豐富，涉 及的學(xué)科和領(lǐng)域也相當(dāng)廣泛。 AR模型的 Burg 法也存在問題，比如計(jì)算量大；信號(hào)起始相位變動(dòng)可導(dǎo)致譜線偏移和分裂 ；低信噪比可導(dǎo)致譜分辨率下降等等。 利用求得的 AR系數(shù)先得到一個(gè) FIR系統(tǒng)為 ?? ???pkkk zazA1 ?1)(? 序列 )(nx 經(jīng)此 FIR 系統(tǒng)濾波，得到一個(gè)輸出序列 ?? ???pk k knxanxny 1 )(?)()( ， ARMA(p,q)模型與FIR系統(tǒng) )(?zA 級(jí)聯(lián)，近似于模型 )(zB 。因此 p 和 q 的不正確指定有可能導(dǎo)致自相關(guān)陣出現(xiàn)奇異。 ②式中階次 p和 q都是未知的，需要事先指定。從第 q+1 個(gè)方程開始是線性的，可以解出 AR部分的系數(shù)，將上式中的第二 個(gè)方程寫成如下展開形式： ?????????????????????????????????????????????????????)()2()1()()2()1()2()()1()1()1()(21pqrqrqraaaqrpqrpqrpqrqrqrpqrqrqrxxxpxxxxxxxxx???????? 上式雖然可解出 AR 部分的系數(shù)，但存在以下兩個(gè)問題： ①由于式中的真實(shí)自相關(guān)函數(shù) )(mrx 是未知的，因此只能使用估計(jì)值 )(?mrx 來代替，且要用到大延遲的估計(jì)值（最大延遲是 p+q），而對(duì)于給定的信號(hào)長(zhǎng)度，這將造成 )(?mrx 估計(jì)很不準(zhǔn)確。 ARMA 模型譜估計(jì) ARMA(p,q)模型的差分方程 ?? ?? ?????qr rpk k rnwbknxanx 01 )()()( (49) 式中 10?b 。 因?yàn)?? ?)()()( nxkmnwEkmrxw ???? ?????? ???? ??? 0 )()()( i inwihkmnwE )()()( 202 mkhikmihi ????? ??? ??? 對(duì) MA(q)模 型，由式 (2) 式得 qibih i ,1,0)( ??? 。粗略的講，AR譜的方差反比于數(shù)據(jù) )(nxN 的長(zhǎng)度 N 和信噪比 SNR。因?yàn)?1)( )(21 ??? ???? ?? deP eP jARjx，整個(gè)積分值為 1，所以又可得到譜匹配的局部性質(zhì)，也就是說 )( ?jAReP 是 )( ?jx eP 的包絡(luò)的一個(gè)好的近似。 3 .AR 模型譜的匹配性質(zhì) 由維納 辛欽定理及 AR譜對(duì)應(yīng)的一個(gè)無限長(zhǎng)的自相關(guān)函數(shù)可得 2)()()(???jjejx eA ePeP ? (46) 2m in)()( ???jjAR eAeP ? 南京郵電大學(xué) 通達(dá)學(xué)院 2020 屆本科生畢業(yè)設(shè)計(jì)（論文） 若用 AR譜去匹配信號(hào)的譜，則誤差系列的譜 應(yīng)由常數(shù)譜來匹配，體現(xiàn) )(zA 的白化性質(zhì)， 從整體上看 )( ?jAReP 將均勻地和 )( ?jx eP 相跟隨。 2 .AR 譜的分辨率 AR 譜估計(jì)的頻率分辨率，要優(yōu)于經(jīng)典譜估計(jì)方法。 AR 模型譜估計(jì) 譜估計(jì)的步驟 我們可以按如下步驟估計(jì) )(nx 的功率譜： ①由 )(nxN 估計(jì) )(nx 的自相關(guān)函數(shù)，得 ?xr ， m=0,1,....p； ②用 ?xr 代替上述地推算法中的 )(mrx ，重新求解 YuleWalker 方程，這時(shí)求出的 AR 模型參數(shù)是真實(shí)參數(shù)的估計(jì)值，即 ?)1(a 、 ?)2(a , ?)(pa 、 p?? ； 將這些參數(shù)代入（ 3）式，得到 )(nx 的功率譜 )( jwx eP 的估計(jì)。同 Burg 算法一樣，改進(jìn)協(xié)方差算法進(jìn)行功率譜估計(jì)時(shí)令南京郵電大學(xué) 通達(dá)學(xué)院 2020 屆本科生畢業(yè)設(shè)計(jì)（論文） 前后向預(yù)測(cè)誤差功率之和最小，即對(duì) )(neb 、 )(nef 前后都不加窗，但得到的協(xié)方差矩陣不是 Toeplitz 矩陣，因此正則方程不能 用 Levinson 遞推算法求解。 (3)計(jì)算步驟如下：①由初始條件 )()()( 00 nxnene bf ?? ，再由式 ? ?pmnenenenek NmnbmNmnfmNmnbmfmm ,2,1)1()()1()(2?1 211 21111????????????????????? (46) 求出 1?k ； ②由 ????102)(1)0(? Nnx nxNr 得 1?m 時(shí)的參數(shù) ，11,1 ?? ka ? ；)0(?)?1(? 211 xrk??? ③求出 )(1 nef 和 )(1neb ，再 mk? 由估計(jì) 2?k ； ④依照 Levinson 遞推關(guān)系，求 2?m 時(shí)的參數(shù) 1,2?a 、 2,2?a 及 2?? ； ⑤重復(fù)上述過程，直到 pm? ，求出所有階次時(shí)的 AR 參數(shù)。用 Burg 算法進(jìn)行功率譜估計(jì)時(shí)令前后向預(yù)測(cè)誤差功率之和最小，即對(duì) )(neb 、 )(nef 前后都不加窗，使用 Levinson— Durbin 遞推可快速的求解AR系數(shù)。因此自相關(guān)法也是已知所有 AR 系數(shù)求解方法中簡(jiǎn)單的一種，但譜分辨率相對(duì)較差。模型參數(shù)算法就是基于上述最小均方誤差時(shí)由模型參數(shù)估計(jì)信號(hào)功率的方法，主要有以下幾種經(jīng)典算法： )()()( 1 knxkanx pk f ??? ??? )()()( nxnxne f ??? 南京郵電大學(xué) 通達(dá)學(xué)院 2020 屆本科生畢業(yè)設(shè)計(jì)（論文） ? ?2)(neEp ff ? )()()( 1 knxkanx pk b ??? ??? )()()( pnxpnxne b ???? ? ? ?2)(neEp fbb ? (45) (1)自相關(guān)法 (BT 法 )。 根據(jù)線性預(yù)測(cè)理論知：一個(gè) P 階 AR 模型的 p+1 個(gè)參數(shù)同樣可用來構(gòu)成 P階的最佳線性預(yù)測(cè)器，其預(yù)測(cè)的最小均方誤差等于 AR 模型激勵(lì)白噪聲的能量，即AR模型是在最小方差意義上對(duì) 數(shù)據(jù)的擬合。 AR 模型參數(shù)求解的典型算法 用線性方程組的常用解法 (例如高斯消元法 )解 Yule— Walker 方程，需要的運(yùn)算量數(shù)量級(jí)為 3p ，但若利用系數(shù)矩陣的對(duì)稱性和 Toeplitz 性質(zhì)，則可構(gòu)成一些高效算法， Levinson— Durbin 算法是其中最著名、應(yīng)用最廣泛的一種，這種算法的運(yùn)算量數(shù)量級(jí)為 2p 。 將方程兩邊同乘以 x(n+m)，并求得均值，最后得到 )()()( 1 mrkmramr xupk xkx ???? ?? (42) 又因?yàn)? ? ? ??? ???? 0)0( 00)()(2 mhmmnxnuE ? (43) 由 Z變換的定義 ,在 式中，當(dāng) ??z 時(shí)，有 h(0)=1 綜合兩式，有 南京郵電大學(xué) 通達(dá)學(xué)院 2020 屆本科生畢業(yè)設(shè)計(jì)（論文） ??????????????????)0()()1)(()(121mkramkmramr pkxkpkxkx? (44) 在上面的推導(dǎo)中，應(yīng)用了自相關(guān)函數(shù)的偶對(duì)稱性。 AR 模型又稱為自回歸模型，它是一個(gè)全極 點(diǎn)模型，其當(dāng)前輸出是現(xiàn)在輸入和過去輸入的加權(quán)和，表示如下 (其中 u(n)為白噪聲序列； p 為 AR 模型的階數(shù) : )()()( 1 nuknxanXpk k ???? ?? (41) ?? ???? pk kk zazAzH 111)(12)( 由隨機(jī)信號(hào)通過線性系統(tǒng)理論知輸出序列的功率譜 2121)(????? pkjw kkjwxeaep ? 其中 2? 為白噪聲序列的方差，因此進(jìn)行功率譜估計(jì)，必需求得 AR模型的參數(shù) ak(k=l， 2? p)及 2? 。本章將會(huì)詳細(xì)的討論AR模型參數(shù)的計(jì)算、譜的性能及其他算法（如線性預(yù)測(cè)、最大熵譜估計(jì)等）的關(guān)系，最后簡(jiǎn)要給出 MA模型及 ARMA 模型譜估計(jì)算法。不難想象， AR模型易于反映譜中的峰值，而 MA模型易于反映譜中的谷值。工程實(shí)際中所遇到的功率譜大體分為三種，一種是“平滑”，即白噪聲的譜，另一種是“線譜”，這是由一個(gè)或多個(gè)純正弦信號(hào)所組成的信號(hào)的功率譜，這兩種是極端的情況；介于二者之間的是既有峰值又有谷值，這種譜稱為 ARMA 譜。輸出序列 )(nx 可以是平穩(wěn)的隨機(jī)序列，也可以是確定性的時(shí)間序列。 （ 3）由 )(zH 的參