這兩天在做一個《數字信號處理》的作業，主要是功率譜估計方面的知識運用，這里又不得不再次檢討一下自己的磨蹭了，一個期末大作業居然做了整整兩天，簡直慢到懷疑人生。
話不多說，上題目。
一、試驗數據的產生
分別產生兩個零均值的高斯白噪聲數據u1(n)和u2(n)，其功率都為 σ2 = 0.12，讓 u1(n)和u2(n)分別通過一個FIR系統，得到輸出為v1(n)和v2(n)。該FIR系統由 5個FIR子系統級聯而成：
B1 = [1;1.98;0.98]，
B2 = [1; −1.98;0.98]，
B3 = [1; −1.8418;0.98]，
B4 = [1; −1.5;0.98]，
B5 = [1; −1.2727;0.98]，
令v(n) = v1(n) + jv2(n)。
在v(n)上再加上三個復正弦信號，它們的幅度分別為a1=12;a2 =12;a3 =6，歸一化頻率分別為f1 = 0.22; f2 = 0.21; f3 = 0.12，這樣可得到已知功率譜的試驗信號x(n) 。
1. 寫出得到 x(n) 的計算過程。令所得到的數據長度 N = 256，描繪該波形。
2. 描繪出該試驗信號的真實功率譜Px(ejw)。Px(ejw)在一個周期內取4096個點（注：以下凡是功率譜曲線，都這樣取）。
3. 求x(n)中三個正弦信號在上述三個頻率處的信噪比（注：為簡化考慮，此處信噪比僅考慮正弦對白噪聲，而不考慮白噪聲通過系統后的影響）。
二、 對x(n)做功率譜估計
1 用周期圖估計 x(n) 的功率譜，令N = 256，畫出Pper(ejw)曲線。
2. 用 Welch 方法對上述周期圖法作改進，例如，可選每段64點，重疊32點，用 Hamming 窗，畫出所求功率譜曲線。
3. 用自相關法估計x(n)的功率譜，令 M = 63，畫出PBT(ejw) 曲線
4. 用 AR 模型估計 x(n) 的功率譜，階次p可在6到15之間，由自己給定，給出估計的功率譜曲線 PAR( ejw) 。分別用自相關法和Berg算法。
5. 試利用 FPE 或 AIC 准則來確定本信號所需要的 AR 模型的階次。
題解：
1、白噪聲信號的產生
利用rand產生一組白噪聲數據。題目中要求的白噪聲數據長度是N=296。為了保證產生的數據能較好的近似白噪聲，數據的長度要盡量取得長一些，此處取原始數據長度NN=2000。將得到的白噪聲數據先減去均值，再乘以 σ ，即可得到滿足要求的白噪聲數據 。同理，再次利用上述方法可以生成另一組白噪聲 ，因為是計算機隨機產生的兩組數據，所以兩組序列不相關如果程序中取 σ2 =0.12。
2、系統級聯
利用MATLAB的conv函數，可以求得級聯系統為11階FIR系統，系統響應為
歸一化系統響應為：
3、試驗信號的疊加
令u1(n)和u2(n) 分別通過該級聯系統，分別輸出v1(n) ，v2(n) 。注意，在計算u(n) 與h(n)的卷積時可以調用零相濾波器函數filtfilt避免時延對系統的影響。
再令
v(n) = v1(n) + jv2(n)
則 是一個復值的噪聲序列，其實部和虛部也基本不相關，那么 的真實功率譜可以求出，即
其中b(k)是模型的系數。
在v(n)的基礎上加上三個復正弦，其歸一化頻率已知。f1 = 0.22; f2 = 0.21; f3 = 0.12。即
然后可以取中間的數據段截取N=256個點這樣就可以求出試驗信號xn了。
4、求出xn的真實功率譜。復正弦信號的強度是piA2 。因此，xn的真實功率譜是
下圖所示即為xn的真實功率譜，其歸一化頻率范圍是從-0.5~0.5，
5、求x(n)中三個正弦信號在上述三個頻率處的信噪比。根據信噪比的定義
其計算方法是10log(Ps/Pn)，其中Ps和Pn分別代表信號和噪聲的有效功率。
那么：
將數據代入公式得
第一個正弦信號的信噪比SNR1 = 27.7815
第二個正弦信號的信噪比SNR2 = 27.7815
第三個正弦信號的信噪比SNR3 = 21.7609
6、用周期圖法對試驗信號進行譜估計。其原理是把隨機信號xn的N點觀察數據視為能量有限信號，直接進行傅立葉變換，然后再取其幅值的平方，並除以N，具體公式如下
在MATLAB中可以直接調用periodogram函數來計算功率譜密度。
如圖所示，由於主瓣寬度過小，f1和f2不能完全分開，只是在波形的頂部能看出是兩個頻率的分量。
7、用自相關法對試驗信號進行譜估計。其原理是先估計出 xN(n)自相關函數r(m)，然后對r(m)求傅立葉變換得到xN(n)的功率譜，以此作為P(w)的估計。
下圖是自相關法估計的效果圖，顯然譜變得光滑，但是較大的邊瓣使分辨率降低了很多。
8、用Welch法估計試驗信號的功率譜，Welch方法是一種修正周期圖功率譜密度估計方法，它通過選取的窗口對數據進行加窗處理，先 xN(n)分段時允許數據重疊求功率譜之后再進行平均。
MATLAB中可以直接調用pwelch函數估計信號的功率譜，其基本的調用格式
[Pxx,f]=pwelch(x,window,noverlap,nfft,fs)
其中窗口的長度N表示每次處理的分段數據長度，Noverlap是指相鄰兩段數據之間的重疊部分長度。
根據題目要求，選每段64點，重疊32點，用 Hamming 窗，畫出所求功率譜曲線如下圖所示。譜變得更加平滑，但是分辨率相比周期圖法有所降低。
9、用AR模型估計x(n)的功率譜。其最常用的有兩種方法，自相關法和Burg法。自相關法的基本思想是利用Levinson遞推公式求解Yule-walker公式求得的AR模型的參數等效於前向預測器的參數。因此AR模型的自相關法等效於對前向預測序列efp(n)前后加窗；Burg算法的基本思想是直接從觀測的數據利用線性預測器的前向和后向預測的總均方誤差之和為最小的准則來估計反射系數，進而通過Levinson遞推公式求出 AR模型優化參數。
根據本題的數據可以得出下面的AR模型估計譜，其中自相關法估計的階數取的是15階，Burg算法取得是6階，這里可以看出無論是在譜的分辨率方面還是模型的階次方面，Burg算法都比自相關算法有明顯的優點。經過試驗，自相關算法的階次大概取30階時才能把 和 兩個譜峰分開，但是Burg算法在階次等於5的時候就能完全把兩個譜峰完全分開。
10、利用FPE准則來確定本信號所需要的AR模型的階次。因為從低階到高階，階次選得高，模型的最小預測誤差功率是遞減的，固然會提高譜估計的分辨率，但同時會產生虛假的譜峰或譜的細節。 所以必須選擇合適的AR模型的階次。這里可以根據最終預測誤差准則(FPE)來確定AR模型的階次。當階次k由1增加時，FPE(k)將在某一個k處取一個極小值，可以將此時的k定為最合適的階次p。具體的計算公式如下：
在MATLAB中可以調用aryule函數和arburg函數來確定k階預測器的最小預測誤差 ，實際計算
Optimal_FPE_yule =16；Optimal_FPE_burg=49；
但是可以畫出FPE(k)的取值趨勢變化曲線，可以看到，實際上使用Burg算法時，在階次大於5后 的變化不大，所以實際計算中可以直接取6階計算。這也與上一步的試驗分析一致。
