濾波與功率譜估計(jì)

1、清 華 大 學(xué)數(shù)字信號(hào)處理期末作業(yè)2013 年 1 月第一題 掌握去噪的方法1.1 題目描述MATLAB 中的數(shù)據(jù)文件 noisdopp 含有噪聲，該數(shù)據(jù)的抽樣頻率未知。調(diào)出該數(shù)據(jù)，用你學(xué)過的濾波方法和奇異值分解的方法對(duì)其去噪。要求： 1盡可能多地去除噪聲，而又不損害原信號(hào)； 2給出你去噪的原理與方法；給出說明去噪效果的方法或指標(biāo)； 3形成報(bào)告時(shí)應(yīng)包含上述內(nèi)容及必要的圖形，并附上原程序。1.2 信號(hào)特性分析MATLAB所給noisdopp信號(hào)極其頻域特征如圖1.1、圖1.2。圖1.1 含有噪聲的noisdopp信號(hào)圖1.2 noisdopp信號(hào)頻域特性其中橫坐標(biāo)f采用歸一化頻率，即未知抽樣頻率

2、Fs對(duì)應(yīng)2（與濾波器設(shè)計(jì)時(shí)參數(shù)一致）。信號(hào)基本特性是一個(gè)幅值和頻率逐漸增加的正弦信號(hào)疊加噪聲，噪聲為均勻的近似白噪聲，沒有周期等特點(diǎn)。因?yàn)樵肼曅盘?hào)能量在全頻帶均勻分布，濾波器截止頻率過低則信號(hào)損失大，過高則噪聲抑制小，認(rèn)為頻譜中含有毛刺較多的部分即為信噪比較小的部分，濾除這部分可以達(dá)到較好的濾波效果。先給定去噪效果的評(píng)定指標(biāo)。信號(hào)開始階段頻率較高（如圖1.3，紅圈為信號(hào)值），一周期內(nèi)采樣點(diǎn)45個(gè)，即信號(hào)歸一化頻率達(dá)到0.40.5（Fs=2），難以從頻域?qū)⑦@部分信號(hào)同噪聲分離，濾波后信號(hào)損失較大，故對(duì)前128點(diǎn)用信噪比考察其濾波效果，定義：其中，為原nosidopp信號(hào)，為濾波后信號(hào)。SNR越大

3、表示濾除部分能力越小，可以反映濾波后信號(hào)對(duì)原信號(hào)的跟蹤能力，對(duì)前128點(diǎn)主要考察SNR，越大濾波器性能越好。對(duì)128點(diǎn)以后的點(diǎn)，信號(hào)能量集中在低頻部分，濾波后效果不僅僅需考察SNR，同時(shí)考察濾波后信號(hào)的平滑程度。定義平滑指標(biāo)1r越小表示濾波后信號(hào)相比原信號(hào)更加平滑，濾波效果更好。圖1.3 信號(hào)前100點(diǎn)本文采用整體濾波、分時(shí)段濾波和奇異值分解的方法分別對(duì)信號(hào)進(jìn)行去噪，并采用信噪比和平滑指數(shù)兩個(gè)指標(biāo)評(píng)定去噪效果。1.3 濾波器設(shè)計(jì)根據(jù)信號(hào)的頻譜，設(shè)計(jì)IIR濾波器和FIR濾波器對(duì)信號(hào)去噪，用重疊相加法對(duì)信號(hào)進(jìn)行卷積運(yùn)算（每128點(diǎn)截為16段）。根據(jù)信號(hào)的頻譜，均采用低通濾波方法去噪，將濾波器截止頻

4、率設(shè)為0.2。1.3.1 IIR濾波器濾波設(shè)計(jì)巴特沃斯低通濾波器，給定參數(shù)，由buttord得到階數(shù)。用MATLAB函數(shù)butter設(shè)計(jì)濾波器，并得到70點(diǎn)頻率響應(yīng)如圖1.4，濾波器沖擊響應(yīng)如圖1.5。圖1.4 巴特沃斯低通濾波器幅頻特性圖1.5 巴特沃斯低通濾波器沖擊響應(yīng)采用重疊相加法計(jì)算nosidopp信號(hào)于巴特沃斯低通濾波器卷積結(jié)果。nosidopp信號(hào)每128點(diǎn)分為一段，共16段，得到濾波后信號(hào)有1093點(diǎn)。由圖1.5可以得到，設(shè)計(jì)得到的巴特沃斯濾波器存在30點(diǎn)的延時(shí)，取卷積后信號(hào)31點(diǎn)開始的1024點(diǎn)為濾波后信號(hào)y。得到濾波后信號(hào)如圖1.6。與原信號(hào)比較如圖1.7。圖1.6 濾波后信

5、號(hào)圖1.7 濾波后信號(hào)與原信號(hào)直觀看濾波后信號(hào)去除了大部分噪聲，但開始高頻部分信號(hào)損失較多（如圖1.7），后尾部信號(hào)頻率較低，濾波后噪聲含量仍比較多。圖1.8 信號(hào)起始部分前128點(diǎn)信噪比SNR為3.429，后896點(diǎn)信號(hào)SNR為12.715，全信號(hào)平滑指標(biāo)為0.0642。本題中設(shè)計(jì)的IIR濾波器的問題在于不是線性相位，只能近似得到濾波后信號(hào)的延時(shí)，濾波后信號(hào)難以同原信號(hào)“對(duì)齊”沖擊響應(yīng)無限長(zhǎng)，有限長(zhǎng)序列通過IIR系統(tǒng)得到的輸出也是無限長(zhǎng)，存在截?cái)嗾`差分時(shí)段設(shè)計(jì)濾波器時(shí)計(jì)算較慢。1.3.2 窗函數(shù)法設(shè)計(jì)FIR濾波器濾波選擇性能較優(yōu)的漢寧窗設(shè)計(jì)低通濾波器，截止頻率0.2，N取34，得到35點(diǎn)長(zhǎng)度

6、的濾波器。得到濾波器幅頻響應(yīng)和沖擊響應(yīng)如圖1.9，圖1.10圖1.9 低通濾波器幅頻特性圖1.10 低通濾波器沖擊響應(yīng)該低通濾波器具有很好的線性相位，輸出信號(hào)存在N/2=17點(diǎn)的延時(shí)，得到的濾波后信號(hào)如圖1.11。與原信號(hào)比較如圖1.12。圖1.11 濾波后信號(hào)圖1.12 濾波后信號(hào)與原信號(hào)直觀感覺比巴特沃斯濾波器的濾波結(jié)果在信號(hào)尾部紋波減小，信號(hào)開始部分損失也比較大。得到前128點(diǎn)信噪比SNR為6.448，后896點(diǎn)信號(hào)SNR為14.418，全信號(hào)平滑指標(biāo)為0.0643。由計(jì)算得到的SNR和平滑指標(biāo)看，窗函數(shù)法得到的濾波器濾波效果比巴特沃斯濾波器要好。采用矩形窗設(shè)計(jì)同樣階數(shù)的FIR濾波器，得

7、到的濾波效果為：前128點(diǎn)信噪比SNR為10.303，后896點(diǎn)信號(hào)SNR為15.738，全信號(hào)平滑指標(biāo)為0.184，如圖1.13。圖1.13 濾波后信號(hào)與原信號(hào)顯然，全信號(hào)采用同樣的濾波器參數(shù)濾波，信號(hào)高頻分量保持較好，則會(huì)保留較多的噪聲。若在頻域?qū)π盘?hào)直接加窗，相當(dāng)于N=的窗函數(shù)濾波器，恢復(fù)信號(hào)有明顯的吉布斯現(xiàn)象，如圖1.14。圖1.14 濾波后信號(hào)的吉布斯現(xiàn)象1.3.3 一致逼近法設(shè)計(jì)FIR濾波器濾波選擇性切比雪夫一致逼近法設(shè)計(jì)低通濾波器，截止頻率0.2，N取34，通帶、阻帶加權(quán)均取1，得到35點(diǎn)長(zhǎng)度的濾波器。濾波器幅頻響應(yīng)如圖1.15。圖1.15 低通濾波器幅頻特性具有很好的線性相位，

8、同樣，輸出信號(hào)存在N/2=17點(diǎn)的延時(shí)，得到的濾波后信號(hào)如圖1.16。與原信號(hào)比較如圖1.17。圖1.16 濾波后信號(hào)圖1.17 濾波后信號(hào)與原信號(hào)前128點(diǎn)信噪比SNR為6.167，后896點(diǎn)信號(hào)SNR為13.987，全信號(hào)平滑指標(biāo)為0.0655，與窗函數(shù)法得到的濾波效果類似。對(duì)全信號(hào)采用同樣的濾波器參數(shù)濾波，因?yàn)闆]有考慮信號(hào)頻率隨時(shí)間的變化，不能得到更好的濾波效果。1.4 分時(shí)段濾波用短時(shí)傅里葉變化（）估計(jì)信號(hào)的時(shí)頻特性大致如圖1.18。圖1.18 Nosidopp信號(hào)的時(shí)頻特性同樣將信號(hào)分為16段，每段64點(diǎn)，對(duì)每段信號(hào)做離散傅里葉變化求其頻譜，得到幅值最大點(diǎn)頻率為，每?jī)牲c(diǎn)間，因?yàn)樾盘?hào)的

9、包絡(luò)為低頻分量，需要保留，故設(shè)計(jì)為截止頻率的低通濾波器。采用切比雪夫一致逼近法，N取34，通帶、阻帶加權(quán)均取1，輸出信號(hào)存在N/2=17點(diǎn)的延時(shí)。得到的濾波后信號(hào)如圖1.19。與原信號(hào)比較如圖1.20。圖1.19 濾波后信號(hào)圖1.20 濾波后信號(hào)與原信號(hào)圖1.21 濾波后信號(hào)與原信號(hào)前100點(diǎn)前128點(diǎn)信噪比SNR為8.144，后896點(diǎn)信號(hào)SNR為13.734，全信號(hào)平滑指標(biāo)為0.0773，得到了兼顧信號(hào)高頻分量和去噪的效果，濾波效果比不分時(shí)段要好。1.5 奇異值分解方法濾波奇異值分解濾波的基本原理是通過矩陣分解方法，得到信號(hào)的特征值，特征值中較大的k個(gè)分量反應(yīng)了信號(hào)低頻成分，而較小的分量反

10、應(yīng)高頻噪聲，將較小的特征值置零再恢復(fù)信號(hào)，得到去噪的信號(hào)。奇異值分解方法濾波需要考慮兩個(gè)問題一維信號(hào)構(gòu)成用于奇異值分解的矩陣的方法如何去除特征值。本文采用教材式（9.5.18）給定的方式構(gòu)造矩陣，對(duì)長(zhǎng)度N的信號(hào)x(n)，定義矩陣X1為：得到的矩陣為方陣可最大程度利用其特征值，對(duì)nosidopp信號(hào)，N=1024，故選擇L=513，M=512構(gòu)造矩陣。利用svd分解得到特征值，如圖1.21。同時(shí)得到矩陣U，V。圖1.22 矩陣特征值可以看出特征值存在明顯的轉(zhuǎn)折點(diǎn)，將轉(zhuǎn)折點(diǎn)（第30點(diǎn)）以后的特征值值為0?；謴?fù)信號(hào)矩陣：此時(shí)得到的矩陣將不再對(duì)稱，按上式中對(duì)應(yīng)位置得到濾波后信號(hào)為：得到的濾波后信號(hào)如圖

11、1.22。與原信號(hào)比較如圖1.23。圖1.23 濾波后信號(hào)圖1.24 濾波后信號(hào)與原信號(hào)前128點(diǎn)信噪比SNR為6.145，后896點(diǎn)信號(hào)SNR為13.827，全信號(hào)平滑指標(biāo)為0.0708。直觀觀察濾波后信號(hào)質(zhì)量好與濾波器濾波后的信號(hào)。特別對(duì)信號(hào)前100點(diǎn)，如圖1.24，直觀觀察要好于同樣信噪比下濾波器去噪的效果。圖1.25 濾波后信號(hào)與原信號(hào)前100點(diǎn)若采用將特征值中小于均值的數(shù)置零的方式濾波，得到濾波后信號(hào)如圖1.26。圖1.26 濾波后信號(hào)前128點(diǎn)信噪比SNR為10.83，后896點(diǎn)信號(hào)SNR為15.357，全信號(hào)平滑指標(biāo)為0.315。對(duì)噪聲去除減少，對(duì)原信號(hào)特征保持較好。構(gòu)造矩陣時(shí)，

12、減小矩陣的行數(shù)、增加矩陣列數(shù)，或采用文獻(xiàn)2中提出的構(gòu)造矩陣的方法：不能使得到的特征值拐點(diǎn)給清晰。1.6 小結(jié)在以信噪比SNR和輸出信號(hào)平滑指標(biāo)兩者作為濾波效果評(píng)定標(biāo)準(zhǔn)下，采用切比雪夫一致逼近法分時(shí)段設(shè)計(jì)低通濾波器得到了最好的濾波效果，最優(yōu)指標(biāo)為：前128點(diǎn)信噪比SNR為8.144，后896點(diǎn)信號(hào)SNR為13.734，全信號(hào)平滑指標(biāo)為0.0773。直觀觀察奇異值分解去噪后濾波效果最佳。本文去噪方式存在的不足有給定的評(píng)定濾波效果的標(biāo)準(zhǔn)不夠合理，缺乏標(biāo)準(zhǔn)信號(hào)的情況下，信噪比作為標(biāo)準(zhǔn)難以區(qū)分估計(jì)誤差和噪聲，給定的標(biāo)準(zhǔn)不能完全反應(yīng)直觀觀察出的濾波效果缺乏對(duì)于信號(hào)所分段數(shù)L的研究奇異值分解方法中，對(duì)矩陣構(gòu)

13、造方式對(duì)濾波效果的研究不足。第一題濾波器去噪部分MATLAB程序見附錄I，奇異值分解去噪部分MATLAB程序見附錄II。第二題 掌握功率譜估計(jì)的方法第2章2.1 試驗(yàn)數(shù)據(jù)的產(chǎn)生2.1.1 產(chǎn)生x(n)第一步，產(chǎn)生均值為0，功率為0.12的白噪聲序列。利用教材所給出的式(1.10.2)，可以得到給定功率的白噪聲序列，即rand函數(shù)產(chǎn)生的500點(diǎn)均勻分布的白噪聲序列，再乘以常數(shù)得到。第二步，通過MATLAB函數(shù)sos2tf得到5個(gè)FIR子系統(tǒng)級(jí)聯(lián)組成的FIR系統(tǒng)傳遞函數(shù)的系數(shù)firb和fira，firb即為FIR系統(tǒng)的沖擊響應(yīng)h(n)。將第一步得到的白噪聲信號(hào)分別與h(n)卷積，即得到通過FIR系

14、統(tǒng)后的輸出，并截取其中256點(diǎn)數(shù)據(jù)得到信號(hào)。第三步，產(chǎn)生四個(gè)復(fù)正弦信號(hào)，相加，再與v(n)相加，得到x(n)。四個(gè)復(fù)正弦信號(hào)疊加得到的信號(hào)幅值和相位圖如圖2.1，x(n) 幅值和相位圖如圖2.2。四個(gè)復(fù)正弦信號(hào)疊加信號(hào)是周期信號(hào)，周期為276。圖2.1 復(fù)正弦信號(hào)波形圖2.2 x(n)波形2.1.2 求信號(hào)真是功率譜第一步，計(jì)算白噪聲信號(hào)功率譜。第二步，計(jì)算FIR系統(tǒng)頻率響應(yīng)。第三步，計(jì)算v(n)功率譜。由教材式(10.3.2)給出的隨機(jī)信號(hào)通過線性系統(tǒng)理論，得到功率譜為求v(n)功率譜，先計(jì)算自相關(guān)函數(shù)在計(jì)算功率譜，得到第四步，求復(fù)正弦信號(hào)功率譜。設(shè)S(n)為4個(gè)復(fù)正弦信號(hào)的和。S(n)是確

15、定信號(hào)，且是功率信號(hào)，相求自相關(guān)函數(shù)再求功率譜按教材式(10.4.4)，用時(shí)域傅里葉變換模平方再除以時(shí)間長(zhǎng)度得到的功率譜為S(n)是周期信號(hào)，故上式化為一周期N點(diǎn)內(nèi)平均，取02等分取4096點(diǎn)轉(zhuǎn)換為DFT若信號(hào)S(n)存在截?cái)啵盘?hào)長(zhǎng)度N，傅里葉變換長(zhǎng)度M，此時(shí)S(n)功率譜為得到信號(hào)x(n)的功率譜為取02等分取4096點(diǎn)，離散化的功率譜取上式中N=256，M=4096，得到真實(shí)功率譜如圖2.3。圖2.3 信號(hào)x真是功率譜2.1.3 計(jì)算給定頻率點(diǎn)信噪比求出復(fù)正弦信號(hào)在四個(gè)頻率點(diǎn)處的功率。根據(jù)教材式(3.2.32)，幅度和頻率分別為的復(fù)正弦信號(hào)的功率為用離散的功率譜計(jì)算得到的結(jié)果相同，即復(fù)正

16、弦信號(hào)在功率集中在處，幅度為。噪聲信號(hào)在處功率為。故在頻率處的信噪比為信噪比都比較低。2.2 對(duì)x(n)做功率譜估計(jì)2.2.1 周期圖法估計(jì)x(n)的功率譜周期圖法估計(jì)隨機(jī)信號(hào)功率譜的基本原理是根據(jù)的N點(diǎn)觀察值做傅里葉變換，得到，然后取其平方，再除以N，作為對(duì)其真是功率譜的估計(jì)，周期圖譜法估計(jì)的功率譜為令在02等分取4096點(diǎn)，得到離散化的功率譜對(duì)256點(diǎn)做4096點(diǎn)DFT，得到估計(jì)的曲線如圖2.4。圖2.4 周期圖法估計(jì)x功率譜得到比較好的功率譜估計(jì)值，256點(diǎn)實(shí)際分辨率為，故四個(gè)頻率點(diǎn)的譜線均可分開，能量較小處也有所反應(yīng)。2.2.2 Welch方法改進(jìn)周期圖W點(diǎn)漢明窗表達(dá)式為每長(zhǎng)度M=64

17、分一段，得到Bartlett法估計(jì)的每段功率譜表達(dá)式為其中平均后的功率譜作為最后估計(jì)的功率譜令在02等分取4096點(diǎn)，用DFT代替DTFT，可得到離散化的功率譜估計(jì)。取M=64，每段重疊32點(diǎn)，分7段，用上述Welch法得到的功率譜估計(jì)曲線如圖2.5。圖2.5 Welch法估計(jì)x功率譜功率譜較周期譜法更光滑，但是分辨率下降，不能區(qū)分0.23、0.24兩處頻率值的信號(hào)。64點(diǎn)漢明窗的主瓣寬度為，故不能分辨兩個(gè)頻率值，若到達(dá)0.01的分辨率，至少需400點(diǎn)的窗函數(shù)，超過信號(hào)長(zhǎng)度，等同于周期譜法，故改進(jìn)的周期譜法不能達(dá)到更好的功率譜估計(jì)效果。2.2.3 自相關(guān)法估計(jì)功率譜自相關(guān)法估計(jì)功率譜的基本原理

18、是，根據(jù)信號(hào)的N點(diǎn)觀察值，計(jì)算自相關(guān)函數(shù)的估計(jì)值，再對(duì)做傅里葉變化得到估計(jì)的功率譜。復(fù)信號(hào)的自相關(guān)函數(shù)計(jì)算方法為再得到功率譜估計(jì)為根據(jù)上述公式，取M=63，令在02等分取4096點(diǎn)，用DFT代替DTFT，得到的離散的功率譜計(jì)算方式為得到估計(jì)的功率譜曲線如圖2.6（已舍去負(fù)值）。圖2.6 自相關(guān)法估計(jì)x功率譜譜線比周期譜法光滑，但不能分辨0.23、0.24頻率處的兩個(gè)譜線，減小M或使用窗函數(shù)可使譜線光滑，但將更加降低頻率分辨率。2.2.4 AR模型估計(jì)功率譜的自相關(guān)法自相關(guān)法求解AR模型的基本原理是用信號(hào)的N點(diǎn)觀察值計(jì)算得到的自相關(guān)函數(shù)估計(jì)值代替Yule-Walker方程中的，得到AR模型參數(shù)。

19、采用Levinson-Durbin遞推算法計(jì)算。初始條件遞推關(guān)系為用MATLAB函數(shù)levinson可以求解上述方程參數(shù)，得到p階AR模型參數(shù)。自相關(guān)求解AR模型的最小預(yù)測(cè)誤差隨p階數(shù)增加而減小，p為6時(shí)，最小預(yù)測(cè)誤差77.15，p為15時(shí)最小預(yù)測(cè)誤差為43.3，取p為15，得到AR模型系數(shù)，和最小預(yù)測(cè)誤差，從而得到估計(jì)的功率譜N=4096離散得到得到估計(jì)的功率譜曲線如圖2.7。圖2.7 AR模型自相關(guān)算法估計(jì)x功率譜，p=15基本估計(jì)出四處頻率點(diǎn)的復(fù)正弦信號(hào)，但不能分辨0.23和0.24兩處頻率的譜線，信噪比較低的-0.16頻率點(diǎn)處的信號(hào)估計(jì)不佳，若降低階數(shù)p，則根本估計(jì)不出該處的譜線，p取

20、10時(shí)得到的功率譜如圖2.8。圖2.8 AR模型自相關(guān)算法估計(jì)x功率譜，p=102.2.5 AR模型估計(jì)功率譜的Burg算法根據(jù)教所給的Burg算法計(jì)算。設(shè)定初始條件，由得到。由得到m=1時(shí)的參數(shù)。再由求得，由表達(dá)式求得。由求得m=2時(shí)參數(shù)。重復(fù)上述步驟，可得到m=p時(shí)參數(shù)。根據(jù)上述算法，在MATLAB中編寫程序?qū)崿F(xiàn)上述算法。得到p=18時(shí)的參數(shù)，并利用MATLAB函數(shù)arburg驗(yàn)證了計(jì)算結(jié)果的正確性。得到估計(jì)的功率譜如圖2.9。圖2.9 AR模型Burg算法估計(jì)x功率譜，p=18比自相關(guān)法有更好的頻譜，可以分辨0.23和0.24頻率處的譜線，-0.16處的信號(hào)也更加明顯。但Burg算法在p

21、較小時(shí)（如10），得到的功率譜不僅不能分辨較近的頻譜和幅度小的信號(hào)，且譜線的相對(duì)幅度與信號(hào)幅度對(duì)應(yīng)發(fā)生變化，如圖2.10。圖2.10 AR模型Burg算法估計(jì)x功率譜，p=102.2.6 AR模型估計(jì)功率譜的Burg算法根據(jù)信息論準(zhǔn)則利用自相關(guān)算法，k由1開始增加，通過levinson函數(shù)求解最小估計(jì)誤差，計(jì)算，得到最小值，即為最合適的階次p，此時(shí)，最小估計(jì)誤差49.38。p為10時(shí)利用自相關(guān)算法的AR模型的功率譜估計(jì)結(jié)果如前2.8圖。2.3 小結(jié)對(duì)本題中信號(hào)的功率譜估計(jì)不足在于產(chǎn)生的實(shí)驗(yàn)信號(hào)中的白噪聲不是理想白噪聲，不同次實(shí)驗(yàn)產(chǎn)生的功率不穩(wěn)定，兩者的獨(dú)立性也不佳缺乏對(duì)功率譜估計(jì)效果的定量指標(biāo)

22、。2.12.2.4，2.2.6節(jié)相關(guān)MATLAB程序見附錄III，2.2.5節(jié)相關(guān)MATLAB程序見附錄VI。29參考文獻(xiàn) 1陳強(qiáng), 黃聲享, 王韋. 小波去噪效果評(píng)價(jià)的另一指標(biāo)J. 測(cè)繪信息與工程, 2008(5):13-14. 2周江, 楊浩. 分段濾波在心率變異性功率譜分析中的應(yīng)用J. 重慶大學(xué)學(xué)報(bào)(自然科學(xué)版), 2001(4):95-97. 3海蘭萍, 姜占才, 李振起. 一種改進(jìn)的奇異值分解語音降噪算法J. 科技信息, 2012(6):156-158. 4遲華山, 王紅星, 趙培洪, 等. 基于S變換的線性調(diào)頻信號(hào)時(shí)頻濾波J. 無線電通信技術(shù), 2012(1):21-24. 5DE

23、LL. 常用圖像去噪比較與分析Z. 20101. 6奇異值分解降噪的改進(jìn)方法-張磊彭偉才原春暉劉彥Z. 7胡廣書. 數(shù)字信號(hào)處理M. 第二版. 北京: 清華大學(xué)出版社, 2008. 8自適應(yīng)線性調(diào)頻信號(hào)時(shí)變?yōu)V波器在分?jǐn)?shù)傅里葉域Z. 9姚文俊. 自相關(guān)法和Burg法在AR模型功率譜估計(jì)中的仿真研究J. 計(jì)算機(jī)與數(shù)字工程, 2007(10):32-34.10白噪聲背景下正弦信號(hào)頻率和功率估計(jì)的有效算法Z.11李軍, 王凱. 三種現(xiàn)代功率譜估計(jì)方法性能研究J. 科技信息, 2010(30):554-555.附錄I% 期末大作業(yè)第一題clear all;clc;close allload noisdo

24、pp %讀入信號(hào)%信號(hào)特性x=noisdopp;nx=length(x);plot(x);axis tight;xlabel('n');ylabel('noisdopp');grid onplot(1:100,x(1:100),1:100,x(1:100),'ro');axis tight;xlabel('n');ylabel('noisdopp');grid onx_fft=fft(x);figuresubplot(2,1,1);plot(1:nx/2*2./nx,abs(x_fft(1:nx/2);axis t

25、ight;xlabel('f');ylabel('|X(k)|');subplot(2,1,2);plot(1:nx/2*2./nx,unwrap(angle(x_fft(1:nx/2);axis tight;xlabel('f');ylabel('phase X(k)');% 分段處理信號(hào)；每M點(diǎn)分一段M=128/2;% 設(shè)計(jì)濾波器 %FIR，窗函數(shù)法階數(shù)為L(zhǎng)，有L+1點(diǎn)，有L/2點(diǎn)的延時(shí) L=34;fc=0.2; %Fs=2 b1=fir1(L,fc,hanning(L+1);figure;freqz(b1);figure;s

26、tem(b1);xlabel('n');ylabel('h(n)');grid on %FIR，切比學(xué)夫一致逼近法 F=0,fc-0.03,fc+0.03,1; A=1,1,0,0; W=1,1;b2=remez(L,F,A,W);freqz(b2) %將H(1)置為一% b2_mo=max(abs(fft(b2);% b2=remez(L,F,A,W)./b2_mo; %IIR，巴特沃斯濾波器 Niir,wniir=buttord(fc-0.03,fc+0.03,0.1,60);bb,ab=butter(Niir,fc); %bb,ab=bilinear(bb

27、s,abs,2); L_b3=2*L+1;b3=impz(bb,ab,L_b3);freqz(b3);figure;stem(b3);xlabel('n');ylabel('h(n)');grid on %FIR，窗函數(shù)法,吉布斯現(xiàn)象 nb51=1:L/2-1;nb52=L/2+1:L+1;wc_b5=fc*2*pi; b5=sin(nb51-L/2)*wc_b5)./(pi*(nb51-L/2),wc_b5./pi,sin(nb52-L/2)*wc_b5)./(pi*(nb52-L/2);% 用于重建信號(hào)y1=zeros(1,nx+length(b1)-1);

28、y2=zeros(1,nx+length(b2)-1);y3=zeros(1,nx+length(b3)-1);y4=zeros(1,nx+length(b1)-1);y5=zeros(1,nx+length(b5)-1);ylw=zeros(1,nx);for n=1:(nx-M)/(M-n_over)+1 xw=x(n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M); %直接頻域加窗 xw_fft=fft(xw); xl=xw_fft(1:5),zeros(1,length(xw)-10),xw_fft(length(xw)-4:length(xw); yl=iff

29、t(xl,'symmetric'); ylw(n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M)=ylw(n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M)+yl; %分段調(diào)整濾波頻率,取頻率最大值左右一定帶寬的信號(hào) max_fft,n_max_fft=max(xw_fft(1:M/2); %通帶中心頻率 f0=n_max_fft./M*2; %設(shè)置帶寬 fw=3;df=1./M*2; if n_max_fft-fw<=0 %即為低通濾波器，F(xiàn)IR，切比學(xué)夫一致逼近法 F=0,f0,f0+fw*df,1; b4=remez(

30、L,F,A,W); else %F=0,f0-fw*df,f0-fw*df+df,f0+fw*df-df,f0+fw*df,1; %即為帶通濾波器 F=0,f0,f0+fw*df,1;b4=remez(L,F,A,W); end %信號(hào)還原 yw5=conv(xw,b5); y5(n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M+length(b5)-1)=y5(n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M+length(b5)-1)+yw5; yw4=conv(xw,b4); y4(n-1)*(M-n_over)+1:(n-1)*(M-n_

31、over)+M+length(b4)-1)=y4(n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M+length(b4)-1)+yw4; yw3=conv(xw,b3); y3(n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M+length(b3)-1)=y3(n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M+length(b3)-1)+yw3.' yw2=conv(xw,b2); y2(n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M+length(b2)-1)=y2(n-1)*(M-n_

32、over)+1:(n-1)*(M-n_over)+M+length(b2)-1)+yw2; yw1=conv(xw,b1); %信號(hào)幅值不對(duì) y1(n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M+length(b1)-1)=y1(n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M+length(b1)-1)+yw1; endfigure my1,ny_y1=max(b1);y1_yx=y1(ny_y1:length(y1);plot(y1_yx);axis tight;xlabel('n');ylabel('y')

33、;grid onfigure;plot(1:nx,x,1:nx,y1_yx(1:nx),'r-');axis tight;xlabel('n');ylabel('x,y');grid onmy2,ny_y2=max(b2);y2_yx=y2(ny_y2:length(y2);plot(y2_yx);axis tight;xlabel('n');ylabel('y');grid onfigure;plot(1:nx,x,1:nx,y2_yx(1:nx),'r-');axis tight;xlabel(

34、'n');ylabel('x,y');grid on%figure;%plot(ylw) % 直接頻域加窗，等于時(shí)域用sinc函數(shù)截?cái)?，效果不好，吉布斯現(xiàn)象my3,ny_y3=max(b3);y3_yx=y3(ny_y3:length(y3);plot(y3_yx);axis tight;xlabel('n');ylabel('y');grid onplot(1:nx,x,1:nx,y3_yx(1:nx),'r-');axis tight;xlabel('n');ylabel('x,y

35、9;);grid onmy4,ny_y4=max(b4);y4_yx=y4(ny_y4:length(y4);plot(y4_yx);axis tight;xlabel('n');ylabel('y');grid onfigure;plot(1:nx,x,1:nx,y4_yx(1:nx),'r-');axis tight;xlabel('n');ylabel('x,y');grid onmy5,ny_y5=max(b5);y5_yx=y5(ny_y5:length(y5);% plot(y5_yx);axis ti

36、ght;xlabel('n');ylabel('y');grid on% plot(1:nx,x,1:nx,y5_yx(1:nx),'r-');axis tight;xlabel('n');ylabel('x,y');grid on%計(jì)算峰值信噪比% y_MSN=y3_yx(129:nx);% y_MSN128=y3_yx(1:128);% y_MSN=y1_yx(129:nx);% y_MSN128=y1_yx(1:128);% y_MSN=y5_yx(129:nx);% y_MSN128=y5_yx(1:128

37、);% y_MSN=y2_yx(129:nx);% y_MSN128=y2_yx(1:128);y_MSN=y4_yx(129:nx);y_MSN128=y4_yx(1:128);x_MSN128=x(1:128);x_MSN=x(129:nx);MSN=(y_MSN-x_MSN)*(y_MSN-x_MSN)'PSNR=10*log10(x_MSN*x_MSN'/MSN);MSN128=(y_MSN128-x_MSN128)*(y_MSN128-x_MSN128)'PSNR128=10*log10(x_MSN128*x_MSN128'/MSN128);%平滑尺度

38、% y_r1=y3_yx(2:nx);y_r2=y3_yx(1:nx-1);% y_r1=y1_yx(2:nx);y_r2=y1_yx(1:nx-1);% y_r1=y5_yx(2:nx);y_r2=y5_yx(1:nx-1);% y_r1=y2_yx(2:nx);y_r2=y2_yx(1:nx-1);y_r1=y4_yx(2:nx);y_r2=y4_yx(1:nx-1);x_r1=x(2:nx);x_r2=x(1:nx-1);r_phx=(x_r1-x_r2)*(x_r1-x_r2)'r_phy=(y_r1-y_r2)*(y_r1-y_r2)'r_ph=r_phy/r_phx

39、;% 用短時(shí)傅里葉變化估計(jì)信號(hào)時(shí)頻特性% spectrogram(x,128,127,128,2); %計(jì)算自相關(guān)函數(shù)% rx=xcorr(x,x);subplot(3,1,1);plot(rx);xlabel('n');ylabel('r')% ry=xcorr(y,y);subplot(3,1,2);plot(ry);% ry1=xcorr(y1,y1);subplot(3,1,3);plot(ry1);%對(duì)不足一周期正弦信號(hào)FFT分析效果% nt=0:1/300:127/300;% yt=sin(2*pi*nt);% Eb=(b4*b4');% y

40、tl=conv(yt,b4);% plot(yt',ytl(1:length(yt)')%plot(abs(fft(yt)附錄II% 奇異值分解方法降噪clear all;clc;close all;load noisdopp %讀入信號(hào)%信號(hào)特性x=noisdopp;%plot(x);nx=length(x);% LxL=nx方法M=nx/2;X1=zeros(M+1,M);for i=1:M+1 for j=1:M X1(i,j)=x(i+j-1); endendU,S,V=svd(X1);plot(abs(diag(S);axis tight;xlabel('k&

41、#39;);grid on%濾波S1=zeros(i,j);% S1(1:30,1:30)=S(1:30,1:30); %衰減很明顯s_u=sum(diag(S)./length(diag(S);s_n=find(diag(S)-s_u>0);s_nq=max(s_n);S1(1:s_nq,1:s_nq)=S(1:s_nq,1:s_nq); %均值下去除%重建,奇異值分解降噪的改進(jìn)方法A=U*S1*V.'y=0*x;for i=1:nx ytemp=0; if i<=M for j=1:i ytemp=ytemp+A(j,i-j+1); end y(i)=ytemp/i;

42、else for j=1:nx-i+1 ytemp=ytemp+A(M+1-j+1,i-M+j-1); end y(i)=ytemp/(nx-i+1); endendplot(y);axis tight;xlabel('n');ylabel('y');grid on;figureplot(1:nx,x,1:nx,y,'r-');axis tight;xlabel('n');ylabel('x,y');grid on%計(jì)算峰值信噪比y_MSN=y(129:nx);y_MSN128=y(1:128);x_MSN128=

43、x(1:128);x_MSN=x(129:nx);MSN=(y_MSN-x_MSN)*(y_MSN-x_MSN)'PSNR=10*log10(x_MSN*x_MSN'/MSN);MSN128=(y_MSN128-x_MSN128)*(y_MSN128-x_MSN128)'PSNR128=10*log10(x_MSN128*x_MSN128'/MSN128);%平滑尺度y_r1=y(2:nx);y_r2=y(1:nx-1);x_r1=x(2:nx);x_r2=x(1:nx-1);r_phx=(x_r1-x_r2)*(x_r1-x_r2)'r_phy=(y_

44、r1-y_r2)*(y_r1-y_r2)'r_ph=r_phy/r_phx;附錄III% 第二題clear all;clc;close allN=500; %信號(hào)長(zhǎng)度%功率為0.12的白噪聲產(chǎn)生方式(1.10.2)P=0.12;a=sqrt(12*P);%產(chǎn)生噪聲序列u1=rand(1,N)*a;u2=rand(1,N)*a;%5個(gè)FIR子系統(tǒng)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;B=B1;B2;B3;B4;B5;sos=B,ones(5,1),zeros(5,

45、2);%得到級(jí)聯(lián)系統(tǒng)傳遞函數(shù)firb,fira=sos2tf(sos);v1=conv(u1,firb);v2=conv(u2,firb);nx=256;v=v1(100:100+nx-1)+1i*v2(100:100+nx-1);%產(chǎn)生復(fù)正弦信號(hào)exp_a=6,12,12,2;exp_f=0.12,0.23,0.24,-0.16;for i=1:length(exp_a) exp_n=0:nx-1; exp_x(i,1:nx)=exp_a(i)*exp(1i*2*pi*exp_f(i)*exp_n);endexp_s=sum(exp_x);x=exp_s+v;save('x_sv.m

46、at','x') %保存用于burg算法subplot(2,1,1);plot(abs(exp_s);axis tight;xlabel('n');ylabel('|exp|');grid onsubplot(2,1,2);plot(angle(exp_s);axis tight;xlabel('n');ylabel('phase exp/rad');grid on%畫x波形圖figuresubplot(2,1,1);plot(abs(x);axis tight;xlabel('n');yl

47、abel('|x|');grid onsubplot(2,1,2);plot(angle(x);axis tight;xlabel('n');ylabel('phase x/rad');grid onM=4096;nw=-M/2:M/2-1;fw=nw./M;w=nw*(2*pi./M);bt=zeros(1,length(w);for k=0:10 bt=firb(k+1).*exp(-1i*w*k)+bt;endPv=abs(bt).2*2*P;%Ps=1./length(w)*(abs(fft(exp_s,length(w).2);Ps=z

48、eros(1,length(w);for i=1:4 expndt=find(2*pi*exp_f(i)-w<0); expd=min(expndt); Ps(expd)=length(x).2./length(w)*exp_a(i).2;end% plot(Ps)Px=Pv+Ps;figure;plot(fw,Px);grid on;axis tight;xlabel('f');ylabel('Px');%計(jì)算信噪比SNR=10*log10(exp_a.2./2/P);%周期圖法功率譜估計(jì)%求信號(hào)FFTXper=fftshift(fft(x,M);Ppe

49、r=abs(Xper).2./M;figureplot(fw,Pper);grid on;axis tight;xlabel('f');ylabel('Pper');%用Welch法估計(jì)mw=64;nxw=length(x);lw=(nxw-mw/2)/(mw/2);%得到漢明窗wnw=0:mw-1;wn_h=0.54-0.46*cos(2*pi*wnw/mw);Uw=1./mw*(wn_h*wn_h');%將x分段Pper_w=zeros(1,M);for i=1:lw xni=(x(i-1)*mw/2+1:(i+1)*mw/2).*wn_h; xni

50、_fft=fft(xni,M); Pwt=1./(M*Uw)*(abs(xni_fft).2); Pper_w=Pper_w+Pwt;endPper_w=1./lw*fftshift(Pper_w)*16;figure;plot(fw,Pper_w);grid on;axis tight;xlabel('f');ylabel('Pper Welch');%自相關(guān)法估計(jì)m_bt=63;n_bt=length(x);%計(jì)算自相關(guān)函數(shù)r_bt=zeros(1,m_bt);for i=0:m_bt-1 r_bt(i+1)=1./n_bt*sum(x(1:n_bt-i)'.'.*x(1+i:n_bt);endPbt=(2*r