现代信号处理大作业.docx
《现代信号处理大作业.docx》由会员分享，可在线阅读，更多相关《现代信号处理大作业.docx（16页珍藏版）》请在冰豆网上搜索。
现代信号处理大作业
现代信号处理大作业
姓名：
潘晓丹
学号：
0140349045班级：
A1403492
作业1
LD算法实现AR过程估计
1.1AR模型
p阶AR模型的差分方程为：
p
x（n）•yaix（n—i）=w（n）,其中w（n）是均值为o的白噪声。
i二
Yule-Walker方程递推
AR过程的线性预测方法为：
先求得观测数据的自相关函数，然后利用
求得模型参数，再根据公式求得功率谱的估计。
Yule-Walker方程可写成矩阵形式:
J（0）
rxx
（1）
rxx
（2）…
rxx（-P）■
-
-11
f
rxx
（1）
「xx（O）
rxx
（1）…
rxx（-P+1）
ap
（1）
0
rxx
（2）
rxx
（1）
rxx（0）…
「xx（-P+2）
ap
（2）
=
0
］「xx（p）
rxx（P-1）
「xx（P—2）…
rxx（0）一
ap（P）_
ii
0一
1.2LD算法介绍
Levinson-Durbin算法可求解上述问题，其一般步骤为：
1）计算观测值各自相关系数rxx（j）,j二0,1,…，P；5二rxx（0）；i=1；
2）利用以下递推公式运算：
i」
rxx（i）+迟a」（j）订』-j）
ai（i）二ki
ai（j）-aiJ（j）-kiaiJ（i-j）,j-1,2>...>i-1
2
:
i=:
i二（1-ki）
3）i=i+1，若i>p，则算法结束；否则，返回⑵。
1.3
以AR模型:
matlab编程实现
为例，Matlab程序代码如下:
clear;clc;
var=1;
noise=var*randn（1,10000）;
p=2;
coefficient=[1-0.50.5];
x=filter（1,coefficient,noise）;
divide=linspace（-pi,pi,200）;
forii=1:
200
w=divide（ii）;
S1（ii）=var/（abs（1+coefficient（2:
3）*exp（-j*w*（1:
2））'）F2;
end
[a_pvar_p]=Levinson_Durbin（x,p）;
forii=1:
200
w=divide（ii）;
Sxx（ii）=var_p/（abs（1+a_p（2:
p+1）*exp（-j*w*（1:
p））'））A2;
end
figure;
subplot（2,2,1）;plot（divide,S1,'b'）;
gridon
xlabel（'w'）;ylabel（'功率'）;title（'AR功率谱'）;
subplot（2,2,2）;plot（divide,Sxx,'r-'）;
gridon
xlabel（'w'）;ylabel（'功率'）;title（'L-D算法估计'）;
subplot（2,2,3）;plot（divide,S1,'b'）;
holdon
plot（divide,Sxx,'r--'）;
holdoff
gridon
xlabel（'w'）;ylabel（'功率'）;title（'AR功率谱和算法比较'）;
子函数：
Levinson_Durbin.m
function[a_pvar_p]=Levinson_Durbin（x,p）
N=length（x）;
forii=1:
N
Rxx（ii）=x（1:
N-ii+1）*（x（ii:
N））'/N;
end
a
（1）=1;
a
（2）=-Rxx
（2）/Rxx
（1）;
fork=1:
p-1%Levinson-Durbinalgorithm
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））A2）*var（k+1）;
a_temp
（1）=1;
forkk=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;
end
a_p=a;
var_p=var（p+1）;
%predictioncoeffecients
%predictionerrorpower
1.4仿真结果
1）
阳功率i亂和LB登法比啟
p=2时，仿真结果图如下
预测系数：
误差功率：
var_p=1.0194
2）p=20时，仿真结果图如下
预测系数:
误差功率：
var_p=0.9998
3）
ARU]^谙和LD負:
去比较
W
p=50时，仿真结果图如下
预测系数:
误差功率：
var_p=0.9955
1.5结果分析
由不同阶数（P值）得到的仿真结果可得：
当P的阶数较低时，L-D算法估计AR模型对功率谱估计的分辨率较低，有平滑的效果，从P=2的仿真结果可以看出估计得到的功率谱与原始功率谱基本吻合，且曲线平滑没有毛刺；随着阶数增大，采用L-D算法进行估计后，得
到的功率谱会产生振荡，从仿真可以看到，当阶数P较高为50时，估计得到的功率谱与原
始功率谱基本吻合，但估计得到的功率谱曲线不平滑，有急剧的振荡。
从LD算法得到的预测系数可得，不论阶数P（>2）取何值，通过该算法得到的预测参数与原始目标函数中的一致，其余各个参数均接近0。
因此，L-D算法得到可计算得到较为精确
的预测值或估计值。
作业
非平稳信号由两个高斯信号叠加而成，为
1
2:
-2
z（t）=（—）[exp（（t-ti）j-it）exp（（t-t2）j2t）]
TL22
其中tl“2」'1J，2，分别求出z（t）的WV分布及其模糊函数，画出二者的波形图，指
出并分析其信号项和交叉项。
2.1WV分布
由WV分布的定义知：
Wt「）z（t）z*（t）ejd.
22
i
若记Zi（t）=（】）4exp（（t-ti）2j.J）,
JT2
1
Z2（t）=
（二）4exp（（t-t2）2j2t）
二2
则W（t,）=WU（t,■）WCross（t,）,
2）Z2*（t-2^d
）z/*（t）e」d
其中，
W/uto（t,■）二"Zi（t）Zi（t）ejd."Z2（t
22
Wross（t,J二"Zi（t）Z2*（t）ejd.■Z2（t
22
由此，我们计算得到：
信号项：
2121
W/uto（t,■）=2exp[-:
（t-t1）（；-；「）]2exp[-:
（t-12）（；•；，2）]
aa
交叉项：
21
Wross（t,）=4exp[-：
（t-tm）（；：
，；：
m）]COS[（；：
，；：
，mXd（t-tm）'d]
CL
t1t2-'1
其中，tm,'m,td-t1_t2/_'2
22
2.2模糊函数
由模糊函数的定义知:
*,可二：
：
z（t）z*（t0如
®22
i
若记Zi（t）=（—）4exp（（t_ti）2jit）
兀2
1
Z2（t）=
（一）4exp（（t—t2）2j2t）
JT2
则Ae8）=Aauto（£8）+Across（y，8），
其中，
Aauto（,讨二「Zi（t—）Zi（t）ej勺dt，Z2（t—）Z2（t）ej^dt
22…22
Across（,刃二Zi（t尹2（t-：
）e心dt「Z2（t；）Zi（t-；）e心dt
22…22
由此，我们计算得到：
信号项：
212
i.-jtp）exp（j2-jtp）]
心o（「）.exp（盲「石'）[exp（
交叉项:
12'212'2
Across（”d）=exp[j（F-「t■dtm）]{exp[
石u■-d）盲（「td）]exp[w（—d）--Ctd）]}
2.3Matlab编程实现
取11=10,t2=4,x=8,，2
4,〉=20，进行matlab编程如下
clearall;
formatlong;
alpha仁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;
fort=0:
0.1:
8
foromega=-6:
0.1:
12
W_auto（m,n）=2*（exp（-a*（t-t1）A2-1/a*（omega-w1F2）+exp（-a*（t-t2）A2-1/a*（omega-w2）A2））;
W_cross（m,n）=4*exp（-a*（t-tm）A2-1/a*（omega-omegam）A2）*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;
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分布交叉项'）;
formatlong;%模糊函数
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;
fort=-10:
0.1:
10
forw=-10:
0.1:
10
A_auto（m,n）=abs（exp（-a/4*tA2-1/（4*a）*wA2）*（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）A2-a/4*（t-td）A2）+exp（-1/（4*a）*（w-wd）A2-a/4*（t+td）A2）））;
A（m,n）=A_auto（m,n）+A_cross（m,n）;
n=n+1;
end
m=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（'模糊函数交叉项'）;
2.4仿真结果及分析
1）Z（t）函数的WV分布波形图如下
2）WV分布信号项波形图如下
卅山玮唁号曲
3）WV分布交叉项波形图如下
■....
i^ciisncr
根据结果可以看到，WV分布的两个信号项是分开的，分别以（ti,「1）和（t2「2）为中心;
WV分布的交叉项则耦合在一起，以（tmJ'm）为中心。
it
4）模糊函数波形图
5）模糊函数信号项
囈舉的劉百号顶
fr»qu«ncy10
6）模糊函数交叉项
frequency'in-10
由结果可以看出，模糊函数的两个信号项耦合在一起，以原点（0,0）为中心；而交叉项是
作业三
相乘得到。
信号-是由一与
分别求出三者的WV分布，并画出三维分布图。
该计算过程和作业二类似，在此不再赘述。
3.1Matlab编程实现
symsxt
a=10;b=2;o=2;
zt=（a/pi）A0.25*exp（-a*tA2/2+1i*b*tA2/2+1i*o*t）;
xt=（a/pi）A0.25*exp（-a*tA2/2）;
yt=exp（1i*b*tA2/2+1i*o*t）;
x1=（a/pi）A0.25*exp（-a*（t+x/2）A2/2）;
x2=（a/pi）A0.25*exp（-a*（t-x/2）A2/2）;
x12=（a/pi）A0.5*exp（-a*（t-x/2）A2/2-a*（t-x/2）A2/2）;
y1=exp（1i*b*（t+x/2）A2/2+1i*o*（t+x/2））;
y2=exp（-1i*b*（t-x/2）A2/2-1i*o*（t-x/2））;
z1=（a/pi）A0.25*exp（-a*（t+x/2）A2/2+1i*b*（t+x/2）A2/2+1i*o*（t+x/2））;
z2=（a/pi）A0.25*exp（-a*（t-x/2）A2/2+-1i*b*（t-x/2）A2/2-1i*o*（t-x/2））;
z12=（a/pi）A0.5*exp（-a*（t+x/2）A2/2+1i*b*（t+x/2）A2/2+1i*o*（t+x/2）-a*（t-x/2）A2/2+-1i*b*（t-x/2）A2/2-1i*o*（t-x/2））;
Wx=simple（fourier（x12）*（（a/pi）A（1/2）））;
Wz=simple（fourier（z12）*（（a/pi）A（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（[-66-66-13]）;
xlabel（'时间'）;
ylabel（'频率'）;
title（'X（t）WV分布'）;
figure;
mesh（T,F,Wzwv）;
axis（[-66-66-13]）;
xlabel（'时间'）;
ylabel（'频率'）;
title（'Z（t）WV分布'）;
3.2仿真结果及分析
1）X（t）函数的WV分布
2）注意到y信号的WV分布是奇异函数，这里便并没有将其画出来。
手算得到其WV分布为：
（）
3）Z（t）函数的WV分布
升级会员 