现代信号处理大作业

姓名：潘晓丹。

学号：0140349045

班级：a1403492

ld算法实现ar过程估计。

p阶ar模型的差分方程为：

其中是均值为0的白噪声。

ar过程的线性**方法为：先求得观测数据的自相关函数，然后利用yule-walker方程递推求得模型参数，再根据公式求得功率谱的估计。

yule-walker方程可写成矩阵形式：

levinson-durbin算法可求解上述问题，其一般步骤为：

1) 计算观测值各自相关系数；；i=1；

2) 利用以下递推公式运算：

3) i=i+1，若i>p，则算法结束；否则，返回(2)。

以ar模型：为例，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(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 on

xlabel('w');ylabel('功率');title('ar 功率谱');

subplot(2,2,2);plot(divide,sxx,'r-')

grid on

xlabel('w');ylabel('功率');title('l-d算法估计');

subplot(2,2,3);plot(divide,s1,'b');

hold on

plot(divide,sxx,'r--'

hold off

grid on

xlabel('w');ylabel('功率');title('ar功率谱和算法比较');

子函数：levinson_

function [a_p var_p] =levinson_durbin(x,p)

n = length(x);

for ii=1:n

rxx(ii) =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;

for kk=1:k

a_temp(kk+1) =a(kk+1)+reflect_coefficient(k+1+1)*a(k+1-kk+1);

enda_temp(k+1+1) =reflect_coefficient(k+1+1);

a = a_temp;

enda_p = a; %prediction coeffecients

var_p = var(p+1); prediction error power

1）p=2时，**结果图如下。

**系数：误差功率：var_p=1.0194

2）p=20时，**结果图如下。

**系数：误差功率：var_p=0.9998

3）p=50时，**结果图如下。

**系数：误差功率：var_p=0.9955

由不同阶数（p值）得到的**结果可得：当p的阶数较低时，l-d算法估计ar模型对功率谱估计的分辨率较低，有平滑的效果，从p=2的**结果可以看出估计得到的功率谱与原始功率谱基本吻合，且曲线平滑没有毛刺；随着阶数增大，采用l-d算法进行估计后，得到的功率谱会产生振荡，从**可以看到，当阶数p较高为50时，估计得到的功率谱与原始功率谱基本吻合，但估计得到的功率谱曲线不平滑，有急剧的振荡。

从ld算法得到的**系数可得，不论阶数p(>2)取何值，通过该算法得到的**参数与原始目标函数中的一致，其余各个参数均接近0。因此，l-d算法得到可计算得到较为精确的**值或估计值。

非平稳信号由两个高斯信号叠加而成，为。

其中，分别求出的wv分布及其模糊函数，画出二者的波形图，指出并分析其信号项和交叉项。

由wv分布的定义知：

若记 ,则，其中，由此，我们计算得到：

信号项：交叉项：

其中，由模糊函数的定义知：

若记 ,则，其中，由此，我们计算得到：

信号项：交叉项：

其中，取，进行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:8

for omega=-6:0.1:12

w_auto(m,n)=2*(exp(-a*(t-t1)^2-1/a*(omega-w1)^2)+exp(-a*(t-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;

endm=m+1;

n=1;end

figure;

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],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;

for t=-10:0.1:10

for w=-10:0.1:

10 a_auto(m,n)=abs(exp(-a/4*t^2-1/(4*a)*w^2)*(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;

endm=m+1;

n=1;end

figure;

mesh([-10: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('frequency');

title('模糊函数交叉项');

1）z(t)函数的wv分布波形图如下。

2）wv分布信号项波形图如下。

3）wv分布交叉项波形图如下。

根据结果可以看到，wv分布的两个信号项是分开的，分别以为中心；

wv分布的交叉项则耦合在一起，以为中心。

4）模糊函数波形图。

5) 模糊函数信号项。

6) 模糊函数交叉项。

由结果可以看出，模糊函数的两个信号项耦合在一起，以原点（0,0）为中心；而交叉项是分开的，分别以为中心。

信号是由与相乘得到。分别求出三者的wv分布，并画出三维分布图。

该计算过程和作业二类似，在此不再赘述。

syms xt

a=10;b=2;o=2;

zt=(a/pi)^0.25*exp(-a*t^2/2+1i*b*t^2/2+1i*o*t);

xt=(a/pi)^0.25*exp(-a*t^2/2);

yt=exp(1i*b*t^2/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/2));

y2=exp(-1i*b*(t-x/2)^2/2-1i*o*(t-x/2));

现代信号处理大作业

现代信号处理大作业

现代信号处理作业

现代信号处理作业

其他用户还读了