现代信号处理大作业王成志1

《现代信号处理》大作业
姓名：王成志
学号：1140349078
一. L D 迭代算法的matlab 实现
1.1 Levinson-Durbin 算法介绍
功率谱估计大致可以分为经典谱估计和现代功率谱估计，经典谱估计方法存在着以下三点缺陷：
（1）数据加窗或自相关加窗，都隐含着假定在窗外未观测到的数据或自相关系数为零，该假设不切实际。

（2）要性能好往往需要较长的数据，但实际数据长度有限
（3）窗函数容易造成谱的模糊。

采用AR 模型的现代谱估计方法可以克服这些不足。

其中LD 递推算法可以在计算机上方便实现。

LD 递推算法具体计算步骤如下：
（1） Yule-Walker 方程的矩阵形式(1)所示：
⎥⎥
⎥⎥
⎥
⎦⎤
⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢
⎢⎢⎢⎣⎡--+-----001)0()2()1()()1()1()0()1()()2()
1()0(2,1,
σk k k xx xx xx xx x xx xx xx xx xx xx xx a a r k r k r k r k r r r r k r r r r 系数矩阵xx H
xx R R =,为Hermitian 矩阵，对角线上元素相同，即为Topliez 矩阵。

（2） P-1阶Yule-Walker 方程为：
2
1
111(0)
(1)(1),1(1)(0)(2)0,1(1)(2)
(0)0x x x p p x x x p x x x R R R p a R R R p a p R p R p R σ-----⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥
⎢⎥-⎢
⎥⎢⎥⎢⎥=⎢
⎥⎢⎥⎢⎥⎢
⎥⎢⎥⎢
⎥---⎢⎥⎣⎦⎣
⎦⎣⎦ 其中，2
2
11{()}p p E e l σ--=为误差功率。

写成联立方程：
2
1
11,0,0()0,1,,1p p
p k x
k m a R m k m p σ---=⎧=-=⎨=-⎩∑ 取共轭得：
2
1
*
*
11,0
,0()0,1,,1p p
p k
x
k m a
R m k m p σ---=⎧=-=⎨=-⎩∑
变量替换，并利用*
()()x x R l R l =得：
2
1
*
11,10
,1()0,0,,2
p p
p p k
x k m p a
R m k m p σ-----=⎧=--=⎨
=-⎩∑ 表示成矩阵：
*1*1
210(0)
(1)(1),1
0(1)(0)(2),2(1)(2)
(0)1x x x p x x x p p x x x R R R p a p R R R p a p R p R p R σ-----⎡⎤⎡⎤⎡⎤-⎢⎥
⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥=⎢⎥⎢
⎥⎢⎥⎢⎥⎢
⎥⎢
⎥--⎢⎥⎣⎦⎣⎦⎣
⎦ 求解得：
*.1,1,,0,
,p k p k p p p k a a K a k p ---=+=
22*
1p p p p K σσ-=+∆ 22
10p p p K σ-=∆+
,p p p K a =
2
22*22
111[][1]p p p p p p p K K K σσσσ---=+-=-
（3） 当k=1时，即一阶递推为：
⎥⎦
⎤
⎢⎣⎡=⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡-01)0()1()1()0(211,1σa R R R R x x x x
求解可得：
)
1()0()0()
1( ,11,1211,10,1x x x x R a R R R a a +=-
==σ
（4） 对于2≥p 时，递推为:
10,≡p a , *
,1,1,k p p p k p k p a
K a a ---+=, ]1[2
21
2p p p K -=-σσ 21
,-∆-
==p p
p p p a K σ
∑-=--+
=∆1
1
,1)()(p k x k
p x p k p R a
p R
矩阵R x 已知，可得到各阶AR 模型系数为：
)0())1(1( ,)0()1()1(2
111xx xx xx r a r r a -=-
=ρ
1
1
111
)
()()()(--=--∑-+-
=∆-
=k k l xx k xx k k
k l k r l a k r k a ρρ
1,,2,1)
()()()(*
11-=-+=--k i i k a k a i a i a k k k k
12
))(1(--=k k k k a ρρ
1.2实验结果
(1) 输入p=3，rr = [70,60,50,40] 时，求得AR 模型估计参数为:
a =
1.0000 -0.8571 0 0 1.0000 -0.5275 -0.3846 0 1.0000 -0.7572 -0.6996 0.5972 各阶求得的方差为：
sigma = 18.5714 15.8242 10.1801
3阶时，a 3 (1)= -0.7572 a 3 (2)= -0.6996 a 3 (3)= -0.5972
(2) 输入p=5，rr = [30,45,26,33,47,43]时，AR 模型估计参数为：
a =
1.0000 -1.5000 0 0 0 0 1.0000 0.2800 -1.1867 0 0 0 1.0000 0.8227 -1.3147 -0.4573 0 0 1.0000 1.9708 1.9858 -
2.5226 -2.5105 0 1.0000 1.0869 1.0977 -1.8235 -1.8166 0.3521 各阶求得的方差为： sigma =
37.5000 15.3067 12.1054 64.1881 56.2316
5阶时， a 5 (1)= 1.0869 a 5(2)= 1.0977 a 5(3)= -1.8235 a 5(4)= -1.8166 a 5(5)= 0.3521
二． 一维平稳信号由两个高斯信号叠加而成
12241122()()[exp(())exp(())]22
z t t t j t t t j t ααα
ωωπ=--++--+，其中12,
t t >12ωω>，分别求出()z t 的WV 分布及其模糊函数，画出二者的波形
图，指出并分析其信号项和交叉项。

(1) 信号()z t 的WV 分布求取公式为
2(,)()()22j f z W t f z t z t e d πτττ
τ
∞
*--∞=+-⎰
在这里我们使用两种方法来对WV 分布进行求取，分别为手工求法和matlab 求法。

A. 求解过程
由二次叠加原理知，Wigner-Ville 分布的自项和交叉项分别为
12(,)(,)(,)auto z z W t w W t w W t w =+ 1221,,(,)(,)(,)cross z z z z W t w W t w W t w =+
式中 14
2111exp ()2z t t jw t ααπ⎛⎫⎡⎤
=--+ ⎪⎢⎥⎝⎭⎣⎦
1
4
2222exp ()2z t t jw t ααπ⎛⎫⎡⎤
=--+ ⎪⎢⎥⎝⎭⎣⎦
由Wigner-Ville 分布的定义知
111(,)()()22j w z W t w z t z t e d τττ
τ
∞
*--∞=+-⎰
12
2211exp ()()4j w w t t d ααττατπ∞
-∞⎛⎫
⎡⎤=----- ⎪
⎢⎥⎝⎭
⎣⎦
⎰ 在积分公式(
)22
exp 2AC B Ax Bx C dx A ∞
-∞
⎛⎫
--±-=
- ⎪⎝
⎭⎰中 令4
A α
=
，()1B j w w =-和()2
1C t t α=-，则得
()()122111(,)2exp z W t w t t w w αα⎡
⎤=----⎢⎥⎣⎦
类似地，我们有()()222221(,)2exp z W t w t t w w αα⎡
⎤=----⎢⎥⎣⎦
故()()2
221
1(,)2exp auto i i i W t w t t w w αα=⎡
⎤=----⎢⎥⎣⎦∑
又由交叉项定义知 12,12(,)()()22
j w z z W t w z t z t e d τττ
τ∞
*--∞=+-⎰
()()()()222121212121exp ()4222t t jw j w w t t t t j w w t d ααατττ∞-∞⎧⎫⎡⎤⎡⎤=
-+--++--+-+-⎨⎬⎢⎥⎣⎦⎣⎦⎩⎭
于是，积分公式(
)22
exp 2AC B Ax Bx C dx A ∞
-∞
⎛⎫
--±-=
- ⎪⎝
⎭⎰中的参数为
4
A α
=
()()1212114
2
2B t t j w w w α
⎡⎤=
--
-+⎢⎥⎣⎦
()()()22
12122C t t t t j w w t α⎡⎤=
-+-+-⎣
⎦
由此得()()(){}
1222,1(,)2exp exp z z m m m d d W t w t t w w j w w t w t αα⎡
⎤=------+⎡⎤⎣⎦⎢⎥⎣⎦
故Wigner-Ville 分布的交叉项
12,(,)2Re (,)cross z z W t w W t w ⎡⎤=⎣⎦
()()()2214exp cos m m m d d t t w w w w t w t αα⎡
⎤=-----+⎡⎤⎣⎦⎢⎥⎣⎦
式中()1212m t t t =
+和()121
2
m w w w =+
分别为两个谐波信号时延的平均值和频率的平均值，而12d t t t =-和
12d w w w =-，分别是两个谐波信号时延差与频率差。

将上述自项和交叉项相加就可以得到WV 分布
(,)(,)(,)z auto cross W t w W t w W t w =+
()()()()()2
22221112exp 4exp cos i i m m m d d i t t w w t t w w w w t w t αααα=⎡⎤⎡
⎤=----+-----+⎡⎤⎣⎦⎢⎥⎢⎥⎣⎦⎣⎦∑
B. 程序实现
由上面步骤可以发现，手工求取WV 分布非常繁琐。

在这里，我使用matlab 来求取WV 分布并画图。

为了进行比较，我还同时画出了WV 分布信号项的图和交叉项的图。

程序运行结果为
当参数值121210,3,2,3,2t t w w α=====时，WV 分布表达式：
Wz=502186294043021/1407374883553280*exp(-10*t^2)*pi^(1/2)*5^(1/2)*2^(1/2)*(exp(50*t-505/8+5/2*i+i*t+1/2*w-i*w-1/10*w^2)+exp(40*t-202/5-1/10*w^2+2/5*w)+exp(-1/10*w^2-909/10+3/5*w+60*t)+exp(50*t-505/8-5/2*i-i*t+1/2*w+i*w-1/10*w^2))
得到的结果图为：
WV 分布图 WV 分布信号项图
WV 分布交叉项图
当设置如上参数时，WV 分布图中信号项和交叉项混叠在一起。

我们只能把信号项图和交叉项图分别显示出来才能看清楚。

当参数值121235,3,2,3,2t t w w α=====时，WV 分布表达式：
Wz=93950452832135/492581209243648*exp(-35*t^2)*35^(1/2)*pi^(1/2)*(exp(210*t-11034/35-1/35*w^2+6/35*w)+exp(140*t-4904/35-1/35*w^2+4/35*w)+exp(175*t-3065/14+5/2*i+i*t+1/7*w-i*w-1/35*w^2)+exp(175*t-3065/14-5/2*i-i*t+1/7*w+i*w-1/35*w^2))
程序中，Wz 中原来的分量f 用2w f π=⨯⨯代替。

得到的结果图为：
WV 分布图 WV 分布信号项图
WV 分布交叉项图
当设置如上参数时,从图4可以看到, 信号项部分和交叉项部分被分了开来，结合上一组参数的经验，我们很容易对哪部分是信号项和哪部分是交叉项做出判断。

而图5和图6的分布图正好证实了我们的判断。

当我们对其他参数进行改变时，WV 分布图也会发生变化。

(2) 信号()z t 的模糊函数求取公式为
2(,)()()22j tv z A v z t z t e dt
πττ
τ∞
*-∞=+-⎰
和求取WV 分布类似，我同样使用两种方法来求取模糊函数，分别为手工求法和matlab 求法。

A ． 求解过程
信号的()z t 模糊函数的信号项和交叉项分别为
121221,,(,)(,)(,)(,)(,)(,)
auto z z cross z z z z A A A A A A τθτθτθτθτθτθ=+=+
由模糊函数公式知
1111
2
222111(,)()()22
exp exp 242j t z A z t z t e dt
j t t j t t dt θττ
τθααθτωταααπ∞
*--∞∞-∞=+-⎡⎤⎛⎫⎛⎫⎛⎫=+-++- ⎪ ⎪ ⎪⎢⎥⎝⎭⎝⎭⎝⎭⎣⎦
⎰⎰
又由积分公式(
)22
exp 2AC B Ax Bx C dx A ∞
-∞
⎛⎫--±-=
- ⎪⎝⎭
⎰代入得： ()122111(,)exp exp 4
4z A j t ατθτθωτθα⎡⎤
⎛⎫=-+-⎡⎤ ⎪⎢⎥⎣⎦⎝⎭⎣⎦ 类似地，我们有()222221(,)exp exp 44z A j t ατθτθωτθα⎡⎤
⎛⎫=-+-⎡⎤ ⎪⎢⎥⎣⎦⎝⎭⎣⎦ 故可以得到：
()(){}
2211221(,)exp exp exp 4
4auto A j t j t ατθτθωτθωτθα⎡⎤
⎛⎫=-+-+-⎡⎤⎡⎤ ⎪⎢⎥⎣⎦⎣⎦⎝⎭⎣⎦ 对于交叉项，我们有：
(
)()()()()12,1222
112222212121212,22exp exp 2222221exp 224j t
z z j t
A z t z t e dt t t j t t t j t e dt
t t t j t t t t t t θθτττθατταττωωαααθωω∞
*-∞
∞-∞⎛⎫⎛⎫=+- ⎪ ⎪⎝⎭⎝⎭
⎡⎤⎡⎤⎛⎫⎛⎫⎛⎫⎛⎫-+-++⨯-----⎢⎥⎢⎥ ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭⎝⎭⎢⎥⎢⎥⎣⎦⎣
⎦⎡⎤⎡⎤-++++---+++⎢⎥⎣⎣⎦⎰()1212j dt ωωτ∞-∞⎧⎫++⎨⎬⎦⎩⎭
于是，积分公式(
)22
exp 2AC B Ax Bx C dx A ∞
-∞
⎛⎫
--±-=
- ⎪⎝⎭
⎰中的参数为 ()()()()1222
121
4
2
4d d m A B t t j C t t t t j αα
θωαωτ
==++
+⎡⎤=
-+++⎣
⎦
由此得()()()1222,1
(,)exp exp 44z z d d m m d m A t j t t ατωθωτωτθωα⎡⎤=-+--++⎡⎤⎣⎦⎢⎥⎣⎦ 类似地有：
()()()2122,1
(,)exp exp 44z z d d m m d m A t j t t ατωθωτωτθωα⎡⎤=---+++⎡⎤⎣⎦⎢⎥⎣⎦ 所以模糊函数的交叉项
()()()()()1221,,2222(,)(,)(,)
11exp exp 4444exp cross z z z z d d d d m m d m A A A t t j t t τθτθτθααθωτθωτααωτθω=+⎧⎫⎡⎤⎡⎤=-+--+---+⎨⎬⎢⎥⎢⎥⎣⎦⎣⎦⎩⎭⨯++⎡⎤⎣⎦
所以，该信号的模糊函数为：
的平均值，而12d t t t =-和12d ωωω=-，分别是两个谐波信号时延差与频率差。

B ． 程序实现
在这一部分，我使用matlab 来求取模糊函数并画图。

为了进行比较，我还同时画出模糊函数信号项的图和交叉项的图。

程序运行结果为
当参数值12121,2,1,2,1t t w w α=====时，模糊函数表达式：
Az=5081767996463981/18014398509481984*exp(-1/4*t^2)/pi^(1/2)*(exp(i*t+i*v-1/4*v^2)+exp(-1/2*t-1/2-3/2*i+3/2*i*t+1/2*v+3/2*i*v-1/4*v^2)+exp(2*i*t+2*i*v-1/4*v^2)+exp(1/2*t-1/2+3/2*i+3/2*i*t-1/2*v+3/2*i*v-1/4*v^2))
程序中，我们使用t 来代替手工求法中的τ，用v 来代替手工求法中的θ。

得到的结果图为：
模糊函数幅度分布图 模糊函数信号项幅度分布图
模糊函数交叉项幅度分布图 当参数值121210,2,1,2,1t t w w α=====时，模糊函数表达式：
Az=502186294043021/5629499534213120*exp(-5/2*t^2)/pi^(1/2)*10^(1/2)*(exp(5*t-101/40+3/2*i+3/2*i*t-1/20*v+3/2*i*v-1/40*v^2)+exp(-5*t-101/40-3/2*i+3/2*i*t +1/20*v+3/2*i*v-1/40*v^2)+exp(i*t+i*v-1/40*v^2)+exp(2*i*t+2*i*v-1/40*v^2)) 程序中，我们使用t 来代替τ
得到的结果图为：
模糊函数幅度分布图 模糊函数信号项幅度分布图
模糊函数交叉项幅度分布图
三． 令信号
221-t 422()j t j t z t e e e αβαπΩ⎛⎫= ⎪⎝⎭是由X (t )=21-t 42e ααπ⎛⎫ ⎪⎝⎭ 和 Y (t )=2
02j t j t e e βΩ相乘而得到的，Y(t)是线性调频Chirp 信号，分别求出X （t ）,Y(t)和Z(t)的WV 分布并画出三者的三维分布图。

A ． 求解过程
WV 分布的计算公式如下：
2(,)()()22j f z W t f z t z t e d πτττ
τ∞*--∞=+-⎰
带入X （t ）,Y(t)和Z(t)后
得到X （t ）的WV 分布如下式： ()()()
22212212(,)()()22
22exp exp 22exp 22j f x j f W t f x t x t e d t t e d t j f d K t j f πτπτττ
τττααατπαταπτπδαπ∞*--∞∞--∞∞-∞
=+-⎛⎫⎛⎫⎛⎫⎛⎫+- ⎪ ⎪ ⎪ ⎪⎛⎫⎝
⎭⎝⎭ ⎪ ⎪=- ⎪ ⎪ ⎪⎝⎭ ⎪ ⎪⎝⎭⎝⎭
⎛⎫=-+ ⎪⎝⎭=+⎰⎰⎰
得到Y （t ）的WV 分布如下式： ()()22212200120(,)()()22
22exp exp j 2222exp 2j f y j f W t f y t y t e d t t j t t e d j t f d πτπτττ
τττββατττπατβπτπ∞*--∞∞--∞∞-∞
=+-⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫+- ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎛⎫⎛⎫⎛⎫⎝⎭⎝⎭ ⎪ ⎪ ⎪ ⎪=-+Ω++Ω- ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭ ⎪ ⎪ ⎪ ⎪⎝⎭⎝
⎭⎝⎭⎝⎭⎛⎫=-+Ω+ ⎪⎝⎭⎰⎰()
02K t f δβπ=+Ω+⎰ 得到Z （t ）的WV 分布如下式：
()()()()221200(,)()()22
()()()()2222
exp exp 22j f z j f W t f z t z t e d x t y t x t y t e d t j t f d K t f j t πτπτττττττττααττβπτπδβπα∞*--∞∞**--∞∞-∞
=+-=++--⎛⎫=--+Ω+ ⎪⎝⎭=+Ω+-⎰⎰⎰
在本题中，Y （t ）频率0Ω为30HZ ，调频斜率βπ=，有Nyquist 采样定理得采样频率100s f Hz =。

信号长度选为奇数点901，衰减控制因子2απ=，信号持续时段为5s ~5s -。

B ． 程序实现
运行程序得到X （t ）,Y(t)和Z(t)的WV 三维分布如下：
X(t)的WV 三维分布 Y(t)的WV 三维分布
Z(t)的WV 三维分布
附：实验代码
2.非平稳信号
1
22
4
1122
()()[exp(())exp(())]
22
z t t t j t t t j t
ααα
ωω
π
=--++--+的WV分布
及模糊函数求解代码。