功率譜估計MATLAB實現(xiàn)-圖文

功率譜估計MATLAB實現(xiàn)-圖文一、經典功率譜估計分類簡介1．間接法1958年Blackman和Turkey具體實現(xiàn)，即先由N，估計出自相關函數(shù)的估計。，求自相關函數(shù)傅里葉變換，以此變換結果作為對功2．直接法直接法功率譜估計是間接法功率譜估計的一個特例，又稱為周期圖法，它是把隨機信號的N個觀察值對功率譜直接進行傅里葉變換，得到的估計。，然后取其幅值的平方，再除以N，作為3．改進的周期圖法將NLM，然后計算L據(jù)段的周期圖的平均，作為功率譜的估計，以此來改善用N點觀察數(shù)據(jù)直接計算的周期圖特性。根據(jù)分段方法的不同，又可以分為Welch法和Bartlett法。的方差WelchBartlett所分的數(shù)據(jù)段互不重疊，選用的數(shù)據(jù)窗是矩形窗。二、經典功率譜估計的性能比較仿真結果為了比較經典功率譜估計的性能，本文采用的信號是高斯白噪聲加兩個正弦信號，采樣率F=1000Hz，兩個正弦信號的頻率分別為f1=200Hz，f2=210Hz。所用數(shù)據(jù)長度N=400.仿真結果如下：(a)(b)(c)(d)(e)Figure1經典功率譜估計的仿真結果Figure1(a)示出了待估計信號的時域波形；2(f)Figure2(b)示出了用該數(shù)據(jù)段直接求出的周期圖，所用的數(shù)據(jù)窗為矩形窗；Figure2(c)是用BT法（間接法）求出的功率譜曲線，對自相關函數(shù)用的平滑窗為矩形窗，長度M=128，數(shù)據(jù)沒有加窗；Figure2(d)是用BT法（間接法）求出的功率譜曲線，對自相關函數(shù)用的平滑窗為Hamming窗，長度M=64，數(shù)據(jù)沒有加窗；Figure2(e)是用Welch平均法求出的功率譜曲線，每段數(shù)據(jù)的長度為64點，重疊32點，使用的Hamming窗；Figure2(f)是用Welch平均法求出的功率譜曲線，每段數(shù)據(jù)的長度為100點，重疊48點，使用的Hamming窗；性能比較1)直接法得到的功率譜分辨率最高，但是方差性能最差，功率譜起伏劇烈，容易出現(xiàn)虛假譜峰；2)間接法由于使用了平滑窗對直接法估計的功率譜進行了平滑，因此方差性能比直接法好，功率譜比直接法估計的要平滑，但其分辨率比直接法低。3)Welch平均周期圖法是三種經典功率譜估計方法中方差性能最好的，估計的功率譜也最為平滑，但這是以分辨率的下降及偏差的增大為代價的。3．關于經典功率譜估計的總結功率譜估計，不論是直接法還是間接法都可以用FFT物理概念明確，因而仍是目前較常用的譜估計方法。2π/N，N方差性能不好，不是真實功率譜的一致估計，且N4)平均主要是用來改善周期圖的方差性能，但往往又減小了分辨率和增加了偏差，沒有一個窗函數(shù)能使估計的功率譜在方差、偏差和分辨率各個方面都得到改善，因此使用窗函數(shù)只是改進估計質量的一個技巧問題，并不能從根本上解決問題。三、AR模型功率譜估計ARAR模型功率譜估計是現(xiàn)代譜估計中最常用的一種方法，這是因為AR模型參數(shù)的精確估計可以用解一組線性方程（Yule-Walker方程）的方法求得。其核心思想是：將信號看成是一個p階AR過程，通過建立Yule-Walker方程求解AR模型的參數(shù)，從而得到功率譜的估計。由于已知的僅僅是長度有限的觀測數(shù)據(jù)，因此AR模型參數(shù)的求得，通常是首先通過某種算法求得自相關函數(shù)的估計值，進而求得AR模型參數(shù)的估計值。常用的幾種AR模型參數(shù)提取方法有：1)自相關法0向預測誤差最小。3由此估計的自相關矩陣式正定的，且具有Toeplitz性，可以用Levion-Durbin算法求解。2)協(xié)方差法0據(jù)的前向預測誤差最小。由此估計的自相關矩陣式半正定的，且不具有ToeplitzAR3)修正的協(xié)方差法0據(jù)的前向預測誤差具有ToeplitzARAR4)Burg在約束ARLevion-DurbinAR仍然用前面的仿真信號，取AR模型的階數(shù)p=48，用上述各種AR型參數(shù)提取方法估計的功率譜如figure2Figure2AR模型功率譜估計的仿真結果AR分別采用直接法和AR模型功率譜估計中的自相關法得到的上述信號的功率譜估計如figure3所示：4Figure3AR模型功率譜估計與經典功率譜估計比較小結：由于AR滑。AR的分辨率。關于ARFigure4ARAR功率譜估計中的BurgFigure4AR小結：階數(shù)越高，得到的譜的分辨率也越高，但方差也越大，將會產生更多的虛假譜峰；5四、附錄本文所用來仿真的MATLAB程序如下：%ReearchOnClaicPSDEtimateMethod%Author:ChenFeiqiang%Date:2022-12-14%%GenerateSignalclearall;cloeall;clc;randn('tate',0);F=1000;%amplefrequencyt=0:1/F:.3;igma=1;%noievariance某=co(2pit200)+co(2pit210)+igmarandn(ize(t));%generateignal:%coine+noiefigure;plot(t,某),某label('t'),ylabel('某');title('SignalInTimeDomain');gridon;%%EtimatePSDuingperiodogrammethodwindow=rectwin(length(某));%windowfunctionued某n=某'.某window;inde=netpow2(length(某));N=pow2(inde某);某w=fft(某n,N);%doN-pointFFTPw=某w.某conjw)/N;%CalculatePSDk=0:floor(N/2-1);figure;plot(k某F/N,10某log10(Pw(k+1)/ma(Pw)));title('PeriodogramPSDEtimate');某label('Frequency(Hz)');ylabel('NormalizedPSD(dB)'),gridon;holdon%%EtimatePSDuingBTmethodwindow_a=rectwin(length));%windowfunctionfordata某n=某'.某window_a;R某=某corr(某n);%auto-correlationfunctionetimateN=length(R某);M=floor(N/4);%thelengthofmoothwindow%M=100;window_v=rectwin(M);%moothwindowforBTmethodR某Win=R(1:M).某window_v;%moothwindowmultiplyauto-correlationfunction6Pw=ab(fft(R某Win,N));k=0:floor(N/2-1);figure;plot(kF/N,10log10(Pw(k+1)/ma(Pw)),'r');title('BTMethodPSDEtimate');某label('Frequency(Hz)');ylabel('NormalizedPSD(dB)'),gridon;%%EtimatePSDuingWelchmethodwindow=32;%thelengthofeachegmentnoverlap=8;%overlapnumberforeachegmentnfft=pow2(netpow2(length)));%nfft-pointFFTforeachegment[P某,f]=pwelch(某,window