現代信號處理大作業(yè)

1、現代信號處理大作業(yè)姓名：潘曉丹學號：0140349045班級：A1403492作業(yè)1LD算法實現AR過程估計1.1 AR模型 p階AR模型的差分方程為：，其中是均值為0的白噪聲。AR過程的線性預測方法為：先求得觀測數據的自相關函數，然后利用Yule-Walker方程遞推求得模型參數，再根據公式求得功率譜的估計。Yule-Walker方程可寫成矩陣形式：1.2 LD算法介紹Levinson-Durbin算法可求解上述問題，其一般步驟為：1) 計算觀測值各自相關系數；i=1；2) 利用以下遞推公式運算：3) i=i+1，若i>p，則算法結束；否則，返回(2)。1.3 matlab編程實現以A

2、R模型：xn=12xn-1-12xn-2+w(n)為例，Matlab 程序代碼如下：clear; clc;var = 1;noise = var*randn(1,10000);p = 2;coefficient = 1 -0.5 0.5;x = filter(1,coefficient,noise);divide = linspace(-pi,pi,200);for ii = 1:200 w = divide(ii); S1(ii) = var/(abs(1+coefficient(2:3)*exp(-j*w*(1:2)')2;enda_p var_p=Levinson_Durbin(

3、x,p);for ii = 1:200 w = divide(ii); Sxx(ii) = var_p/(abs(1+a_p(2:p+1)*exp(-j*w*(1:p)')2;endfigure;subplot(2,2,1); plot(divide,S1,'b');grid onxlabel('w'); ylabel('功率'); title('AR 功率譜');subplot(2,2,2); plot(divide,Sxx,'r-');grid onxlabel('w'); ylabel

4、('功率'); title('L-D算法估計');subplot(2,2,3); plot(divide,S1,'b');hold onplot(divide,Sxx,'r-');hold offgrid onxlabel('w'); ylabel('功率'); title('AR功率譜和算法比較');子函數：Levinson_Durbin.mfunction a_p var_p = Levinson_Durbin(x,p)N = length(x);for ii=1:N Rxx(i

5、i) = x(1:N-ii+1)*(x(ii:N)'/N; enda(1)=1; a(2)=-Rxx(2)/Rxx(1);for k=1:p-1 % Levinson-Durbin algorithm var(k+1) = Rxx(0+1)+a(1+1:k+1)*Rxx(1+1:k+1)' reflect_coefficient(k+1+1) = -a(0+1:k+1)*(fliplr(Rxx(2:k+1+1)'/var(k+1); var(k+1+1) = (1-(reflect_coefficient(k+1+1)2)*var(k+1); a_temp(1) = 1

6、; for kk=1:k a_temp(kk+1) = a(kk+1)+reflect_coefficient(k+1+1)*a(k+1-kk+1); end a_temp(k+1+1) = reflect_coefficient(k+1+1); a = a_temp;enda_p = a; % prediction coeffecientsvar_p = var(p+1); % prediction error power1.4 仿真結果1）p=2時，仿真結果圖如下預測系數：a20,a21,a22=1,-0.5068,0.5031誤差功率：var_p=1.01942）p=20時，仿真結果圖如

7、下預測系數：a20,a21,a22,a23,a24,=1,-0.5098,0.4999,-0.0066,0.0060,-0.0179,0.0193,誤差功率：var_p=0.99983）p=50時，仿真結果圖如下預測系數：a20,a21,a22,a23,a24,=1,-0.4951,0.5178,-0.0145,0.0117,-0.0169,0.0141,誤差功率：var_p=0.99551.5 結果分析由不同階數（P值）得到的仿真結果可得：當P的階數較低時，L-D算法估計AR模型對功率譜估計的分辨率較低，有平滑的效果，從P=2的仿真結果可以看出估計得到的功率譜與原始功率譜基本吻合，且曲線平滑

8、沒有毛刺；隨著階數增大，采用L-D算法進行估計后，得到的功率譜會產生振蕩，從仿真可以看到，當階數P較高為50時，估計得到的功率譜與原始功率譜基本吻合，但估計得到的功率譜曲線不平滑，有急劇的振蕩。從LD算法得到的預測系數可得，不論階數P(>2)取何值，通過該算法得到的預測參數與原始目標函數中的一致，其余各個參數均接近0。因此，L-D算法得到可計算得到較為精確的預測值或估計值。作業(yè)二非平穩(wěn)信號由兩個高斯信號疊加而成，為其中，分別求出的WV分布及其模糊函數，畫出二者的波形圖，指出并分析其信號項和交叉項。2.1 WV分布由WV分布的定義知：若記,則 ,其中,由此，我們計算得到：信號項：交叉項：其

9、中，2.2 模糊函數由模糊函數的定義知：若記,則 ,其中，由此，我們計算得到：信號項：交叉項：其中，2.3 Matlab 編程實現取，進行matlab編程如下clear all;format long;alpha1=20; t1=10;t2=4;w1=8;w2=4;a=alpha1/2; td=t1-t2; omegad=w1-w2;tm=0.5*(t1+t2); omegam=0.5*(w1+w2); m=1;n=1;for t=0:0.1:8for omega=-6:0.1:12W_auto(m,n)=2*(exp(-a*(t-t1)2-1/a*(omega-w1)2)+exp(-a*(t

10、-t2)2-1/a*(omega-w2)2); W_cross(m,n)=4*exp(-a*(t-tm)2-1/a*(omega-omegam)2)*cos(omega-omegam)*td+omegad*t)W(m,n)=W_auto(m,n)+W_cross(m,n); n=n+1;end m=m+1; n=1;endfigure; mesh(-6:0.1:12,0:0.1:8,W);xlabel('time');ylabel('frequency');title('WV分布');figure; mesh(-6:0.1:12,0:0.1:8,

11、W_auto);xlabel('time');ylabel('frequency');title('WV分布信號項');figure;mesh(-6:0.1:12,0:0.1:8,W_cross);xlabel('time');ylabel('frequency');title('WV分布交叉項');format long; %模糊函數a=10; t1=6;t2=2;w1=6;w2=2;td=t1-t2; wd=w1-w2;tm=0.5*(t1+t2); wm=0.5*(w1+w2);m=1;n=1

12、;for t=-10:0.1:10for w=-10:0.1:10 A_auto(m,n)=abs(exp(-a/4*t2-1/(4*a)*w2)*(exp(i*w1*t-i*t1*w)+exp(i*w2*t-i*t2*w); A_cross(m,n)=abs(exp(i*wm*t+i*w*tm+i*wd*tm)*(exp(-1/(4*a)*(w+wd)2-a/4*(t-td)2)+exp(-1/(4*a)*(w-wd)2-a/4*(t+td)2);A(m,n)=A_auto(m,n)+A_cross(m,n); n=n+1;end m=m+1; n=1;endfigure; mesh(-10

13、:0.1:10,-10:0.1:10,A);xlabel('time');ylabel('frequency');title('模糊函數');figure; mesh(-10:0.1:10,-10:0.1:10,A_auto);xlabel('time');ylabel('frequency');title('模糊函數信號項');figure;mesh(-10:0.1:10,-10:0.1:10,A_cross);xlabel('time');ylabel('frequenc

14、y');title('模糊函數交叉項');2.4 仿真結果及分析1）Z(t)函數的WV分布波形圖如下2）WV分布信號項波形圖如下3）WV分布交叉項波形圖如下根據結果可以看到，WV分布的兩個信號項是分開的，分別以為中心；WV分布的交叉項則耦合在一起，以為中心。4）模糊函數波形圖5) 模糊函數信號項6) 模糊函數交叉項由結果可以看出，模糊函數的兩個信號項耦合在一起，以原點（0,0）為中心；而交叉項是分開的，分別以為中心。作業(yè)三信號zt=0.25e-t2/2ejt2/2ejt是由xt=0.25e-t2/2與yt=ejt2/2ejt相乘得到。分別求出三者的WV分布，并畫出三維分

15、布圖。該計算過程和作業(yè)二類似，在此不再贅述。3.1 Matlab編程實現syms x ta=10; b=2; o=2;zt=(a/pi)0.25*exp(-a*t2/2+1i*b*t2/2+1i*o*t);xt=(a/pi)0.25*exp(-a*t2/2);yt=exp(1i*b*t2/2+1i*o*t);x1=(a/pi)0.25*exp(-a*(t+x/2)2/2);x2=(a/pi)0.25*exp(-a*(t-x/2)2/2);x12=(a/pi)0.5*exp(-a*(t-x/2)2/2-a*(t-x/2)2/2);y1=exp(1i*b*(t+x/2)2/2+1i*o*(t+x/

16、2);y2=exp(-1i*b*(t-x/2)2/2-1i*o*(t-x/2);z1=(a/pi)0.25*exp(-a*(t+x/2)2/2+1i*b*(t+x/2)2/2+1i*o*(t+x/2);z2=(a/pi)0.25*exp(-a*(t-x/2)2/2+-1i*b*(t-x/2)2/2-1i*o*(t-x/2);z12=(a/pi)0.5*exp(-a*(t+x/2)2/2+1i*b*(t+x/2)2/2+1i*o*(t+x/2)-a*(t-x/2)2/2+-1i*b*(t-x/2)2/2-1i*o*(t-x/2);Wx=simple(fourier(x12)*(a/pi)(1/2);Wz=simple(fourier(z12)*(a/pi)(1/2);Wx=inline(Wx)Wz=inline(Wz) t=-6:0.01:6; f=-6:0.01:6;T,F=meshgrid(t,f);Wxwv= abs(Wx(T,2*pi*F); Wzwv= abs(Wz(T,2*pi*F);figure; mesh(T,F,Wxwv); axis(-6 6 -6 6 -1 3);xlabel(&#