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

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

3、者共同特點是由信號觀測數據求模型系數時采用信號預測誤差最小的原則。對于長記錄數據，這些方法的估計質量是相似的，但對于短記錄數據，不同方法之間存在差別。2.2 自相關法列文森（Levenson）遞推法自相關法的出發點是選擇AR模型參數使預測誤差功率最小，預測誤差功率為假設信號的數據區在范圍，有個預測系數，個數據經過沖激響應為的濾波器，輸出預測誤差的長度為，因此應用下式計算：的長度長于數據的長度，上式中數據的兩端需補充零點，相當于對無窮長的信號加窗處理，得到長度為N的數據。上式對系數的實部和虛部求微分使預測誤差功率最小，得到 （1）式中自相關函數采用有偏自相關估計，即對比上式，可知式（1）即為已推

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

5、算這種方法計算簡單，但需要預先估計出信號自相關函數，實際中只能按照信號的有限個觀測數據估計自相關函數。當觀測數據長度較短時，估計誤差較大，會出現譜峰頻率偏移和譜線分裂（在信號譜峰附近產生虛假譜線）；如數據很長，估計自相關函數較準確，但計算量大，應適當選擇數據長度。2.3 伯格（Burg）遞推法Levinson-Durbin遞推法需要由觀測數據估計自相關函數，這是它的缺點。而伯格遞推法則由信號觀測數據直接計算AR模型參數。伯格遞推法利用Levinson-Durbin遞推公式，導出前向預測誤差與后向預測誤差，并按照使它們最小的原則求出，從而實現不用估計自相關函數，直接用觀測數據得出結果。Burg遞

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

7、（1） 編寫程序產生題目要求的信號和噪聲（2） 然后分別用兩種方法的遞推流程進行譜估計（3） 改變題目中要求的變量參數，分析結果的變化4. 代碼Levensonclc;clear all;fs=100; %采樣頻率Ts=1/fs;N=27; %數據長度p1=20; %階數 f1=0.2*fs;f2=0.25*fs; %設置信號頻率pha1=0;pha2=0;%初始相位SNR=2; %設置信噪比%產生信號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)

8、;xx=x1+x2;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;%計算自相關函數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=

9、1:p-1 sum=sum+a(p-1,k)*R(p-k+1); end a(p,p)=-(R(p+1)+sum)/cov(p-1);%計算反射系數 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 %計算功率譜，功率譜以2*pi為周期，又信號為實信號，只需輸出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)

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

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= numer

12、ator+ (-2)*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

13、=2:p1 for 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); %角頻率坐標換算成頻率坐標figure(1)plot(F,pxx);xlabel(Frequency(Hz);ylabel(Power Spectral Density);title(伯格（Burg）遞推法的功率譜估計);grid;% figure(2)% p=1:p1;% plo

14、t(p,cov(p),b);% xlabel(模型階次);% ylabel(預測誤差功率大小);% title(預測誤差功率的變化趨勢);% grid;5. 實驗結果及分析5.1 產生信號5.2 Levenson遞推法1. 數據長度的影響（階數設置為20階）N=32N=64N=128N=1024N=8192N=163842. 模型階數的影響（數據長度設置為N=1024）預測誤差功率的變化（從1階到50階）不同階數（p1表示）功率譜的圖像如下：P1=5P1=10P1=20P1=30P1=503. 相位的影響（設置N=128 p1=20）0Pi/6Pi/3Pi/2Pi4. 頻率的影響（N=128

15、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. 數據長度的影響N=32N=64N=256N=1024N=8192N=16384 2. 模型階數的影響（N=128）預測誤差功率的變化趨勢不同階數的功率譜P1=5P1=10P1=20P1=503. 相位的影響（N=256 p1=20）0Pi/