現(xiàn)代數(shù)字信號(hào)處理-Levinson-Durbin法和Burg法

1、現(xiàn)代數(shù)字信號(hào)處理實(shí)驗(yàn)報(bào)告Levinson-Durbin 算法與 Burg 算法目錄Levinson-Durbin 算法 1一、實(shí)驗(yàn)要求: 1二、實(shí)驗(yàn)原理:1三、編程實(shí)現(xiàn)2四、仿真結(jié)果及分析3Burg 法 5一、實(shí)驗(yàn)要求:5二、實(shí)驗(yàn)原理:5三、編程實(shí)現(xiàn):6四、仿真結(jié)果及分析7附錄（完整代碼）9Levinson-Durbin 算法： 9Burg 法： 10Levinson-Durbin 算法:使用burg 法計(jì)算 x(n)的功率譜:80 階估計(jì)，數(shù)據(jù)采樣分別取128 和 512。u(n)是均值為其中wn 是高斯白噪聲，階數(shù)選用二、實(shí)驗(yàn)原理:實(shí)際應(yīng)用的隨機(jī)過程大多數(shù)可以用有理傳輸函數(shù)模型很好逼近，輸

2、入激勵(lì)零、方差為 2的白噪聲序列，線性系統(tǒng)傳輸函數(shù)為：其中bk 是前饋(或動(dòng)平均)支路的系數(shù)，ak 是反饋(或自回歸)支路的系數(shù)，當(dāng)除了b0=1 之外其余的MA 系數(shù)都為0，則p 階自回歸模型，即AR 模型 :此時(shí)傳輸函數(shù)為：傳輸模型的輸出功率譜為：Wold 分解定理認(rèn)為任何廣義平穩(wěn)隨機(jī)過程都可分解為一個(gè)完全隨機(jī)的部分和一個(gè)確定的部分，如果功率譜完全連續(xù)，那么任何ARMA 過程或者M(jìn)A 過程都可以用一個(gè)無限階AR 過程表示。而Yule-Walker 方程則把模型p階系數(shù)和方差2的關(guān)系：也就是只要估計(jì)出P+1 個(gè)自相關(guān)函數(shù)值，就可以由該式解出P+1 個(gè)模型參數(shù)。Levinson-Durbin 算

3、法的運(yùn)算數(shù)量級(jí)為p2.這是一種按階次進(jìn)行遞推的算法，即首先以AR(0) 和 AR(1) 模型參數(shù)作為初始條件，計(jì)算AR(2) ，如此類推。由k 階到 k+1 階的計(jì)算公式為 :11三、編程實(shí)現(xiàn)仿真工具使用MATLAB2015a ，根據(jù)上述的算法原理，代碼可以分為信號(hào)采樣模塊、自相關(guān)系數(shù)計(jì)算模塊、初始值模塊、算法迭代模塊、功率譜計(jì)算模塊和畫圖模塊，接下來進(jìn)行依次分析。信號(hào)采樣模塊如下，其中N 是采樣數(shù)據(jù)個(gè)數(shù)，利用wgn 函數(shù)直接產(chǎn)生高斯白噪聲，wgn(1,N,0) 表示隨機(jī)產(chǎn)生一個(gè)強(qiáng)度為0dB 的 1 行 N 列的高斯白噪聲矩陣，矩陣n 記錄了 N個(gè)時(shí)間值，矩陣x 則記錄了N 個(gè)輸入信號(hào)采樣值。

4、n0: N 1 ;wwgn 1,N ,0 ;x cos 0.3 * pi * n cos 0.32 * pi * n w;自相關(guān)函數(shù)計(jì)算模塊如下，因?yàn)殡A數(shù)為80，所以需要計(jì)算R(0)R(80) 的值，此處直接采用 xcorr 函數(shù)進(jìn)行自相關(guān)系數(shù)的求解，b, c = xcorr(x, p, 'biased') 表示計(jì)算x 序列時(shí)間差為-p, p 的自相關(guān)系數(shù)值，并依次返回時(shí)間差值到矩陣b 中， 而返回相應(yīng)相關(guān)系數(shù)到矩陣c中，取 0, p 所對(duì)應(yīng)的相關(guān)系數(shù)值，即得相關(guān)系數(shù)矩陣r。p80;b,c xcorr x, p,' biased ' ;r;for i 1: le

5、ngth bif c i 0r rbi;endend初始值模塊如下，此模塊用于模型參數(shù)的初始化，a 表示 AR 參數(shù)， s 表示2， d 表示D，g 表示 。且02=R(0)， D0=R(1)， a1,0=1，所以= D0/ 02=R(1)/R(0) ，進(jìn)而a1,1=- 。此處需要注意的是，MATLAB 矩陣標(biāo)號(hào)按列從1 開始，但是2、 D、 a的下標(biāo)均從0 開始。a zeros p, p 1 ;szeros 1, p1;dzeros 1, p;gzeros 1,p;a 1,11;s1r1;d1r2;g 1 r 2/r 1;a 1,2 g 1 ;算法迭代模塊如下，首先確定k-1 階的a、 2、

6、 D 和 ，然后根據(jù)原理部分的公式依次計(jì)算 k 階的2、 D 和 ， 再根據(jù) k 階的 計(jì)算 k 階的 a 參數(shù)， 最后再根據(jù)p 階 2計(jì)算p+1 階 2。for i 2 : ps(i) (1 g(i 1)2)* s(i 1);d i0;for k1: id i d i a i 1,k * r i 2 k ;endg idi /si ;a i,11;for m1: i 1a i,m 1 a i 1,m 1 g i * a i 1,i m 1 ;enda i,i 1 g i ;ends(p 1) (1 g(p)2)* s( p);功率譜計(jì)算模塊如下，在-pi, pi范圍內(nèi)的2000 個(gè)等間隔頻率

7、點(diǎn)上均勻采樣，直接用linspace 函數(shù)生成，接著利用書中的式4-27 計(jì)算功率譜。W linspace pi, pi,2000 ;S zeros 1,length W ;for i 1: 2000spot W(i);Ssum 1;for k 1: pSsum Ssum a( p, k 1)* exp( j * spot * k);endS(i) s(p 1)/ (abs(Ssum) 2);end以 512 點(diǎn)估計(jì)為例，畫圖模塊如下，此處畫出輸入信號(hào)以及利用Levinson-Durbin 算法求得的功率譜，將兩幅圖顯示在一個(gè)figure 中。subplot 1,2,1 ;plot n,x ;

8、title ('512點(diǎn)信號(hào)');subplot 1,2,2 ;plot W,S ;title ('512點(diǎn) 80階 Levinson估 計(jì) ');由圖可知，大致在0.9 與 0.12 附近有峰值，這也與原題目中的0.3pi 與 0.32pi 相對(duì)應(yīng)，由于噪聲的存在，使得除了峰值點(diǎn)的其余頻率的功率譜不為0，但接近于0，這是因?yàn)榇a中設(shè)置的高斯白噪聲的強(qiáng)度為0dB 比較小。 另外可以看到N=128 和 N=512 對(duì)應(yīng)的峰值不一樣，這也是由于編程方式使得噪聲在取了不同的n 時(shí)噪聲信號(hào)又重置了。urg 法:使用 burg 法計(jì)算x(n)的功率譜:其中wn 是高斯白噪

9、聲，分別使用50 階和 80 階估計(jì)。二、實(shí)驗(yàn)原理:在實(shí)際應(yīng)用中，常需根據(jù)信號(hào)的有限個(gè)取樣值來估計(jì)AR 模型的參數(shù)，應(yīng)用較多的是Yule-Walker 法或自相關(guān)法，協(xié)方差法，Burg 法。與自相關(guān)法的計(jì)算效率很高，而保證預(yù)測誤差濾波器是最小相位的，但數(shù)據(jù)兩端要附加零取樣值，實(shí)質(zhì)上等于數(shù)據(jù)加窗，這會(huì)使參數(shù)估計(jì)的精度下降，當(dāng)數(shù)據(jù)段很短時(shí)，加窗效應(yīng)更嚴(yán)重，是一種直接估計(jì)AR 參數(shù)的方法。Burg 法則一方面希望利用已知數(shù)據(jù)兩端以外的未知數(shù)據(jù)(但它相對(duì)于這些未知數(shù)據(jù)不做主觀臆測)，另一方面有設(shè)法保證預(yù)測誤差濾波器是最小相位的，Burg 法與自相關(guān)法和協(xié)方差法不同，它不直接估計(jì)AR 參數(shù)，二是先估計(jì)反

10、射系數(shù)，然后利用Levinson 遞推算法由反射系數(shù)求得AR 參數(shù)。Burg 法首先要估計(jì)反射系數(shù)，所使用的準(zhǔn)則是前向和后向預(yù)測誤差功率估計(jì)的平均值最小準(zhǔn)則。在這里，預(yù)測誤差功率估計(jì)仍然用時(shí)間平均來代替集合平均。因此，Burg 法估計(jì)反射系數(shù)的準(zhǔn)則為：前向和后向預(yù)測誤差濾波器的工作都是在數(shù)據(jù)段上進(jìn)行的數(shù)據(jù)段兩端不需要補(bǔ)充零。一般情況下，求出k 后，即可使用Levinson 遞推算法中的公式4.25 由 k-1 階 AR 參數(shù)計(jì)算出k 階 AR 參數(shù)。 Burg 法的公式可歸納如下：Burg 法估計(jì) AR( p) 模型參數(shù)的具體計(jì)算方法如下：1) 確定初始條件22) 確定 k-1 階 AR 參數(shù)

11、(迭代計(jì)算時(shí)，k 的值從 1 開始選取)Ak-1(z)，2k-1,k-1nN -1;3) 計(jì)算 k;4) 計(jì)算A k(z);5) 計(jì)算和 ， knN -1;6) 計(jì)算k 階均方誤差7) 回到第二部，進(jìn)行下一次迭代。一般來說，如果處理的數(shù)據(jù)來自AR 過程，那么采用Burg 法可以獲得精確地AR 譜估計(jì)。三、編程實(shí)現(xiàn):仿真工具同樣使用MATLAB2015a ，根據(jù)原理，代碼也可以分為信號(hào)采樣模塊、初始值模塊、 算法迭代模塊、功率譜計(jì)算模塊和畫圖模塊，其中信號(hào)采樣模塊、功率譜計(jì)算模塊以及畫圖模塊與Levinson-Durbin 算法代碼大致相同，代碼中各字母含義也與之前類似，此處不再贅述，下面將著重

12、分析初始值模塊、算法迭代模塊。要計(jì)算的初始參數(shù)包括e0+(n)、 e0-(n)、 02、，其中e0+(n) 和 e0-(n)均等于采樣的信號(hào)矩陣，用矩陣erP 和 erN 表示，02則等于采樣矩陣的平方均值，根據(jù)e0+(n) 和 e0-(n)可以計(jì)算出 ，根據(jù)和a1,0=1 進(jìn)一步可計(jì)算出a1,1，具體代碼如下。此處同樣需要注意的是，MATLAB 矩陣標(biāo)號(hào)按列從1 開始，但是2、 a 、 e0+(n) 和 e0-(n)的下標(biāo)均從0 開始。erP zeros( p 1,N); erN zeros( p 1,N);erP(1,:) x; erN(1,:) x;g zeros(1, p); s ze

13、ros(1, p 1);for i 1: Ns(1) s(1) x(i)2;ends(1) s(1)/ N;a zeros( p, p 1);a(1,1) 1;temp1 0;temp2 0;for i 2 : Ntemp1 temp1 2* x(i)* x(i 1);temp2 temp2 x(i)2 x(i 1)2;endg(1) temp1/ temp2;a(1,2)g(1);算法迭代模塊如下，首先確定k-1 階的a、 2、 ek-1+(n)、 ek-1-(n)， 然后根據(jù)k-1 階的ek-1+(n)、ek-1-(n)計(jì)算k 階的ek+(n)、 ek-(n)和，再根據(jù)k 階的 計(jì)算 k

14、階的 a 參數(shù)和2，最后再根據(jù)p階 2 計(jì)算p+1 階 2。for i 2 : perP(i,1) x(1);erN (i,1) x(1);for k 2 : NerP(i,k) erP(i 1,k) g(i 1)* erN (i 1,k 1);erN(i,k) erN(i 1,k 1) g(i 1)* erP(i 1,k); end 2s(i) (1 g(i 1)2)* s(i 1);temp1 0; temp2 0;for m (i 1): N22temp1 temp1 erP (i , m) erN (i , m 1) ;temp2 temp2 2* erP(i, m)* erN (i,

15、 m 1);endg (i) temp2 / temp1;a(i,1) 1;for m 2 : ia(i,m) a(i 1,m) g(i)* a(i 1,i m 2);enda(i,i 1)g(i);ends( p 1) (1 g(p)2)* s( p);由圖可知，大致在 0.9與 0.12 附近有峰值，但由于相離較近難以辨別，這也與原題目中的 0.3pi 與 0.32pi 相一致，一高一低現(xiàn)象是隨機(jī)造成的，因?yàn)閱未蔚男蛄羞M(jìn)行計(jì)算會(huì)有隨機(jī)性，但總體是會(huì)趨向于一致。由于是按照書上4.8.3 節(jié)標(biāo)準(zhǔn)公式所寫，放大可以發(fā)現(xiàn)，AR 譜估計(jì)中存在兩個(gè)靠的很近的譜峰，即有譜線分裂現(xiàn)象，Burg 算法得到的

16、譜估計(jì)，其譜峰位置移動(dòng)有可能達(dá)到原位置的16%，如果要進(jìn)行修正，那么可以選擇相應(yīng)的算法進(jìn)行改進(jìn)，例如復(fù)數(shù)據(jù)和反射系數(shù)修正等等。而AR 模型階數(shù)的確定，可以用最終預(yù)測誤差（FPE）準(zhǔn)則或者AIC 信息準(zhǔn)則來確定，階數(shù)選的太低，功率譜顯得平滑，階數(shù)太高，會(huì)出現(xiàn)虛假峰，因此需要選擇合適的階數(shù)，此試驗(yàn)中階數(shù)已經(jīng)確定了。Levinson-Durbin 算法：% 清除工作空間和變量空間 clear ;close all;clc;% 信號(hào)采樣%采樣數(shù)據(jù)個(gè)數(shù)N = 512;n = 0:(N-1);%高斯白噪聲w = wgn(1,N,0);%輸入信號(hào)x(n)x = cos(0.3 * pi * n) + cos

17、(0.32 * pi * n) + w;% 計(jì)算自相關(guān)系數(shù)矩陣，R( 0) -R(80)p = 80;b,c = xcorr(x,p,'biased');r = ;for i = 1:length(b)if c(i) >= 0r = r b(i);endend% 算法迭代a = zeros(p,p+1);% a 參數(shù)2s = zeros(1,p+1);% d = zeros(1,p);% Dkg = zeros(1,p);%初始化參數(shù)a(1,1) = 1;s(1) = r(1);d(1) = r(2);g(1) = r(2)/r(1);a(1,2) = -g(1);%迭代

18、for i = 2:ps(i) = (1-g(i-1)2) * s(i-1);for k = 1:id(i) = d(i) + a(i-1,k) * r(i+2-k);endg(i) = d(i)/s(i);a(i,1) = 1;for m = 1:(i-1)a(i,m+1) = a(i-1,m+1) - g(i) * a(i-1,i-m+1);enda(i,i+1) = -g(i);ends(p+1) = (1 - g(p)2) * s(p);% 功率譜W = linspace(-pi,pi,2000);S = zeros(1,length(W);for i = 1:2000spot = W

19、(i);Ssum = 1;for k = 1:pSsum = Ssum + a(p,k+1) * exp(-j * spot * k); endS(i) = s(p+1)/(abs(Ssum)2);End% 畫圖模塊subplot(1,2,1);plot(n,x);title('512 點(diǎn)信號(hào) ');subplot(1,2,2);plot(W,S);title('512 點(diǎn) 80 階 Levinson 估計(jì) ');Burg 法：% 清除工作空間和變量空間clear ;close all;clc;% 信號(hào)產(chǎn)生N = 512;n = 0:(N-1);w = wgn(

20、1,N,0);x = cos(0.3 * pi * n) + cos(0.32 * pi * n) + w;% 初始化參數(shù)p = 80; %設(shè)置階數(shù)erP = zeros(p+1,N); % ek+(n)erN = zeros(p+1,N); % ek-(n)erP(1,:) = x;erN(1,:) = x;g = zeros(1,p); %2s = zeros(1,p+1); % for i = 1:Ns(1) = s(1) + x(i)2;ends(1) = s(1)/N;a = zeros(p,p+1);a(1,1) = 1;temp1 = 0;temp2 = 0;for i = 2:Ntemp1 = temp1 + 2 * x(i) * x(i-1);temp2 = temp2 + x(i)2 + x(i-1)2;endg(1) = temp1 / temp2;a(1,2) = -g(1);% 迭代for i = 2:perP(i,1) = x(1);erN(i,1) = x