功率譜估計(jì)Levinson 遞推法和 Burg 法

1、數(shù)字信號處理實(shí)驗(yàn)報(bào)告姓名： 學(xué)號： 日期：1. 實(shí)驗(yàn)任務(wù)信號為兩個正弦信號加高斯白噪聲，各正弦信號的信噪比均為10dB，長度為，信號頻率分別為和，初始相位，取，取不同的數(shù)值：0.3，0.25。為采樣率。（1） 分別用 Levinson 遞推法和 Burg 法進(jìn)行功率譜估計(jì)，并分析改變數(shù)據(jù)長度、模型階數(shù)對譜估計(jì)結(jié)果的影響。（2） 當(dāng)正弦信號相位、頻率、信噪比改變后，上述譜估計(jì)的結(jié)果有何變化？并作分析說明。2. 原理分析2.1 現(xiàn)代譜估計(jì)中的參數(shù)建模根據(jù)參數(shù)模型來描述隨機(jī)信號的方法，我們可以知道，如果能確定信號的信號模型，根據(jù)信號觀測數(shù)據(jù)求出模型參數(shù)，系統(tǒng)函數(shù)用表示，模型輸入白噪聲，其方差為，信號

2、的功率譜用下式求出：按照這種求功率譜的思路，功率譜估計(jì)可分為三個步驟：（1） 選擇合適的信號模型；（2） 根據(jù)有限的觀測數(shù)據(jù)，或者它的有限個自相關(guān)函數(shù)的估計(jì)值，估計(jì)模型的參數(shù)；（3） 計(jì)算墨香的輸出功率譜。其中以（1）、（2）兩步最為關(guān)鍵。按照模型的不同，譜估計(jì)的方法有許多種，它們共同的特點(diǎn)是對信號觀測區(qū)以外的數(shù)據(jù)不假設(shè)為0，而先根據(jù)信號觀測數(shù)據(jù)估計(jì)模型參數(shù)，按照求模型輸出功率的方法估計(jì)信號功率譜，回避了數(shù)據(jù)觀測區(qū)以外的數(shù)據(jù)假設(shè)問題。下面分析AR譜估計(jì)的兩種方法：自相關(guān)法列文森（Levenson）遞推法和伯格（Burg）遞推法。這兩種方法均為已知信號觀測數(shù)據(jù)，估計(jì)功率譜，兩者共同特點(diǎn)是由信號觀

3、測數(shù)據(jù)求模型系數(shù)時(shí)采用信號預(yù)測誤差最小的原則。對于長記錄數(shù)據(jù)，這些方法的估計(jì)質(zhì)量是相似的，但對于短記錄數(shù)據(jù)，不同方法之間存在差別。2.2 自相關(guān)法列文森（Levenson）遞推法自相關(guān)法的出發(fā)點(diǎn)是選擇AR模型參數(shù)使預(yù)測誤差功率最小，預(yù)測誤差功率為假設(shè)信號的數(shù)據(jù)區(qū)在范圍，有個預(yù)測系數(shù)，個數(shù)據(jù)經(jīng)過沖激響應(yīng)為的濾波器，輸出預(yù)測誤差的長度為，因此應(yīng)用下式計(jì)算：的長度長于數(shù)據(jù)的長度，上式中數(shù)據(jù)的兩端需補(bǔ)充零點(diǎn)，相當(dāng)于對無窮長的信號加窗處理，得到長度為N的數(shù)據(jù)。上式對系數(shù)的實(shí)部和虛部求微分使預(yù)測誤差功率最小，得到 （1）式中自相關(guān)函數(shù)采用有偏自相關(guān)估計(jì)，即對比上式，可知式（1）即為已推導(dǎo)出的Yule-Wa

4、lker 方程，因此自相關(guān)法也是基于解Yule-Walker 方程的一種方法。但是直接解該方程，需要計(jì)算逆矩陣，不方便，因此，基于Yule-Walker 方程中自相關(guān)矩陣的性質(zhì)，導(dǎo)出Levinson-Durbin遞推法，這是一種高效的解方程的方法。Levinson-Durbin算法首先由一階AR模型開始：一階AR模型的Yule-Walker方程為由該方程解出然后令，以此類推，可以得到一般遞推公式如下：稱為反射系數(shù)，。，隨著階數(shù)增加，預(yù)測誤差功率將減少或不變。由k=1開始遞推，遞推到k=p，依次得到各階模型參數(shù)，AR模型的各個系數(shù)及模型輸入白噪聲方差求出后，信號功率譜用下式計(jì)算這種方法計(jì)算簡單，

5、但需要預(yù)先估計(jì)出信號自相關(guān)函數(shù)，實(shí)際中只能按照信號的有限個觀測數(shù)據(jù)估計(jì)自相關(guān)函數(shù)。當(dāng)觀測數(shù)據(jù)長度較短時(shí)，估計(jì)誤差較大，會出現(xiàn)譜峰頻率偏移和譜線分裂（在信號譜峰附近產(chǎn)生虛假譜線）；如數(shù)據(jù)很長，估計(jì)自相關(guān)函數(shù)較準(zhǔn)確，但計(jì)算量大，應(yīng)適當(dāng)選擇數(shù)據(jù)長度。2.3 伯格（Burg）遞推法Levinson-Durbin遞推法需要由觀測數(shù)據(jù)估計(jì)自相關(guān)函數(shù)，這是它的缺點(diǎn)。而伯格遞推法則由信號觀測數(shù)據(jù)直接計(jì)算AR模型參數(shù)。伯格遞推法利用Levinson-Durbin遞推公式，導(dǎo)出前向預(yù)測誤差與后向預(yù)測誤差，并按照使它們最小的原則求出，從而實(shí)現(xiàn)不用估計(jì)自相關(guān)函數(shù)，直接用觀測數(shù)據(jù)得出結(jié)果。Burg遞推法思想：借助格型預(yù)

6、測誤差濾波器，求前向、后向預(yù)測誤差平均功率，選擇使其最小，求出。之后，再利用Levinson-Durbin遞推法求模型參數(shù)和輸入噪聲方差。設(shè)信號的觀測數(shù)據(jù)區(qū)間：,前向、后向預(yù)測誤差功率分別用和表示，預(yù)測誤差平均功率用表示，公式分別為前向、后向觀測誤差公式分別為 上式中，信號項(xiàng)的自變量最大的是n,最小的是n-p，為了保證計(jì)算范圍不超出給定的數(shù)據(jù)范圍，在和計(jì)算公式中，選擇求和范圍為： 。為求預(yù)測誤差平均功率最小時(shí)的反射系數(shù)，令,將前、后向預(yù)測誤差的遞推公式代入得Burg遞推法求AR模型參數(shù)的遞推公式總結(jié)如下：(1) (2) (3) (4) (5) (6) (7) 3. 編程思想（1） 編寫程序產(chǎn)生

7、題目要求的信號和噪聲（2） 然后分別用兩種方法的遞推流程進(jìn)行譜估計(jì)（3） 改變題目中要求的變量參數(shù)，分析結(jié)果的變化4. 代碼Levensonclc;clear all;fs=100; %采樣頻率Ts=1/fs;N=27; %數(shù)據(jù)長度p1=20; %階數(shù) f1=0.2*fs;f2=0.25*fs; %設(shè)置信號頻率pha1=0;pha2=0;%初始相位SNR=2; %設(shè)置信噪比%產(chǎn)生信號w=randn(1,N);Am=sqrt(2*10(SNR/10);k=1:N;x1=Am*sin(2*pi*f1*k*Ts+pha1);x2=Am*sin(2*pi*f2*k*Ts+pha2);xx=x1+x2;

8、x=w+x1+x2;figure(1)subplot(2,1,1);plot(xx);ylim(-20,20);title(兩個正弦信號波形);grid on;subplot(2,1,2);plot(x);ylim(-20,20);title(加噪信號波形);grid on;%計(jì)算自相關(guān)函數(shù)R=;for m=1:N s=0; for n=1:N-(m-1) s=s+x(n)*x(n+m-1); end R(m)=s/N;end %列文森遞推法a(1,1)=-R(2)/R(1); cov(1)=R(1)+a(1,1)*R(2);for p=2:p1 sum=0; for k=1:p-1 sum=

9、sum+a(p-1,k)*R(p-k+1); end a(p,p)=-(R(p+1)+sum)/cov(p-1);%計(jì)算反射系數(shù) for k=1:p-1 a(p,k)=a(p-1,k)+a(p,p)*a(p-1,p-k); end cov(p)=(1-(abs(a(k,k)2)*cov(p-1);end %計(jì)算功率譜，功率譜以2*pi為周期，又信號為實(shí)信號，只需輸出0到P1即可；W=0.01:0.01:pi;Z=0;for k=1:p1 Z=Z+(a(p1,k).*exp(-j*k*W);endPSD=1./(abs(1+Z).2); %算出功率譜F=W*fs/(pi*2); %角頻率坐標(biāo)換算

10、成頻率坐標(biāo)figure(2)plot(F,abs(PSD);xlabel(頻率(Hz);ylabel(功率譜);title(自相關(guān)法-列文森Levenson遞推法的功率譜估計(jì));grid;figure(3)p=1:p1;plot(p,cov(p),b);xlabel(模型階數(shù));ylabel(預(yù)測誤差功率大小);title(預(yù)測誤差功率變化趨勢);grid;Burgclc;clear all;fs=100; %采樣頻率Ts=1/fs;N=28; %數(shù)據(jù)長度p1=20; %階數(shù) f1=0.2*fs;f2=0.25*fs; %設(shè)置信號頻率pha1=0;pha2=0;SNR=2; %設(shè)置信噪比%產(chǎn)生

11、信號w=randn(1,N); % 噪聲為高斯白噪聲，功率為1.Am=sqrt(2*10(SNR/10);k=1:N;x1=Am*sin(2*pi*f1*k*Ts+pha1);x2=Am*sin(2*pi*f2*k*Ts+pha2);x=w+x1+x2;%遞推ef=zeros(p1,N);eb=zeros(p1,N);ef(1,:)=x; eb(1,:)=x; cov(1)=x*x/N; k(1)=0; a=zeros(p1+1,p1+1);for p=2:p1+1 numerator=0; denominator=0; for n=p:N numerator= numerator+ (-2)

12、*ef(p-1,n)*eb(p-1,n-1); denominator=denominator+(ef(p-1,n)2+(eb(p-1,n-1)2; end k(p)=numerator./(denominator+0.0001); cov(p)=(1-abs(k(p)2).*cov(p-1); a(p,p)=k(p); for n=p:N ef(p,n)=ef(p-1,n)+k(p).*eb(p-1,n-1); eb(p,n)=eb(p-1,n-1)+k(p).*ef(p-1,n); endendk=k(1,2:p1+1);a=a(2:p1+1,2:p1+1);for p=2:p1 for

13、i=1:p-1 a(p,i)=a(p-1,i)+k(p)*a(p-1,p-i); endend%功率譜Z=0;W=1:0.01:pi;for i=1:p Z=Z+(a(p,i)*exp(-j*W*i); end;pxx=cov(p1+1)./(abs(1+Z).2);F=W*fs/(pi*2); %角頻率坐標(biāo)換算成頻率坐標(biāo)figure(1)plot(F,pxx);xlabel(Frequency(Hz);ylabel(Power Spectral Density);title(伯格（Burg）遞推法的功率譜估計(jì));grid;% figure(2)% p=1:p1;% plot(p,cov(p)

14、,b);% xlabel(模型階次);% ylabel(預(yù)測誤差功率大小);% title(預(yù)測誤差功率的變化趨勢);% grid;5. 實(shí)驗(yàn)結(jié)果及分析5.1 產(chǎn)生信號5.2 Levenson遞推法1. 數(shù)據(jù)長度的影響（階數(shù)設(shè)置為20階）N=32N=64N=128N=1024N=8192N=163842. 模型階數(shù)的影響（數(shù)據(jù)長度設(shè)置為N=1024）預(yù)測誤差功率的變化（從1階到50階）不同階數(shù)（p1表示）功率譜的圖像如下：P1=5P1=10P1=20P1=30P1=503. 相位的影響（設(shè)置N=128 p1=20）0Pi/6Pi/3Pi/2Pi4. 頻率的影響（N=128 p1=20 初始相位為0）Fs=100, f1=0.2*fs f2=0.25*fsFs=100, f1=0.2*fs f2=0.3*fsFs=1000,f1=0.2*fs f2=0.25*fsFs=1000,f1=0.2*fs f2=0.3*fs5. 信噪比的影響SNR =20dB SNR =10dBSNR =6dBSNR=2dB5.3 Burg遞推法1. 數(shù)據(jù)長度的影響N=32N=64N=256N=1024N=8192N=16384 2. 模型階數(shù)的影響（N=128）預(yù)測誤差功率的變化趨勢不同階數(shù)的功率譜P1=5P1=10P1=20P1=503. 相位的影響（N=256 p1=20）0Pi/