现代信号处理地大作业LD算法以及WV变换.docx
《现代信号处理地大作业LD算法以及WV变换.docx》由会员分享，可在线阅读，更多相关《现代信号处理地大作业LD算法以及WV变换.docx（15页珍藏版）》请在冰豆网上搜索。
现代信号处理地大作业LD算法以及WV变换
现代信号处理大作业
姓名：
学号：
专业：
电子科学与技术
一.L-D算法的仿真实现
1、问题描述
用Matlab实现Levinsion-Durbin算法。
2、算法分析
2.1、LD算法原理
由于语音样点之间存在相关性，所以可以用过去的样点值来预测现在或未来的样点值。
如下图所示
图1线性预测图示
由上图可得
，从而可以通过使实际语音x（n）和线性预测结果
之间的误差e（n）在某个准则下达到最小值来决定唯一的一组预测系数
。
而这组系数就能反映语音信号的特性，可以作为语音信号特征参数来用于语音编码、语音合成和语音识别等应用中去。
由估计值和实际信号值的误差,可有
根据e（n）最小均方误差准则，来决定唯一的一组预测系数
，即:
由此可得到Y-W方程：
，取遍k值之后有以下：
由相关函数的偶函数性质有：
在已知自相关函数的前提下，根据e（n）均方误差最小的原则来求解a，本实验中采用Levinson-Durbin算法。
2.2、LD算法的实现
Levinson-Durbin算法首先由一阶AR模型开始，按照前面的Y-W方程可有，一阶AR模型（p=1）的Y-W方程是
然后增加一阶，即令p=2，可得到：
由上式可解出：
然后令
以此类推，可以得到一般的递推公式：
式中的
称为反射系数，而
和
是预测误差的均方误差值，因此
必须大于等于0，这样
应满足
，进而得到
，即预测误差随递推次数增加而减少。
3、程序实现思路
3.1、观测值的选取。
将均值为0.方差为1，数据长度3000的白噪声信号通过线性系统：
H（z）=1/（1+a1*z^（-1）+a2*z^（-2）），设置模型的输入系数为a0=[10.720.88]；由输出信号作为观测信号x（n）。
3.2、通过LD算法函数估计从一阶到50阶的系数。
3.3、用Matlab自带的aryule函数预测50阶系数，与自己编写的程序结果对比。
4、程序如下：
4.1、主函数：
clc;clear;closeall;
%%给出一个已知的模型，让已知信号经过该模型之后利用函数估计模型系数
%已知模型设为x（n）=a1*noise（n-2）+a2*noise（n-1）+a3*noise（n）
len=3000;%数据长度
noise=randn（1,len）;%产生均值为0，方差（功率）为1，数据长度为3000的高斯白噪声
a0=[10.720.88];%已知模型的输入系数
x=filter（1,a0,noise）;%产生的信号模型是x（n）=a0
（1）*noise（n-2）+a0
（2）*noise（n-1）+a0（3）*noise（n）;
%%利用L-D算法预测模型参数
m=50;
[a,P]=Levinson_Durbin（x,m）;
aa=[ones（m,1）,a];%预测出来的系数
%预测误差功率图
figure
plot（1:
m,P（1:
m）,'rh'）
ylim（[0,6]）
title（'预测误差功率随迭代次数的变化'）;
xlabel（'迭代次数'）
ylabel（'预测误差功率'）
%利用aryule验证
y=aryule（x,m）;
figure
subplot（1,2,1）
plot（1:
51,aa（m,:
）,'r'）
title（'编写程序运行得出的系数'）
subplot（1,2,2）
plot（1:
51,y,'k'）
title（'aryule运行得出的系数'）
4.2、LD算法子函数：
%%Levinson-Durbinalgorithm
function[a,P]=Levinson_Durbin（x,m）
%x为输入信号
%m为阶数
a=zeros（m,m）;
N=length（x）;
P=zeros（1,m+1）;
fori=1:
N
P
（1）=P
（1）+x（i）^2;
end
P
（1）=P
（1）/N;%计算预测误差功率的初始值P_0
f=zeros（m+1,N）;%m阶前向预测误差
g=zeros（m+1,N）;%m阶后向预测误差
f（1,:
）=x;
g（1,:
）=x;
K=zeros（1,m）;
%开始迭代：
fori=1:
m
numerator=0;
denominator=0;
forj=i+1:
N
numerator=numerator-f（i,j）*g（i,j-1）;
denominator=denominator+（f（i,j）^2+g（i,j-1）^2）/2;
end
K（i）=numerator/denominator;%反射系数
forj=1:
i-1
a（i,j）=a（i-1,j）+K（i）*a（i-1,i-j）;
end
a（i,i）=K（i）;
ifi>=2
P（i）=（1-K（i）^2）*P（i-1）;%更新预测误差功率
end
forj=i+1:
N
f（i+1,j）=f（i,j）+K（i）*g（i,j-1）;
g（i+1,j）=K（i）*f（i,j）+g（i,j-1）;
end
end
4.3、仿真结果：
显然，我们想要预测估计的是a0=[10.720.88]，预测结果为（预测了1-50阶，后面阶数过多不易展示，因此只列出前1-6阶预测结果）：
表1：
LD算法预测系数表（1-6阶）
1
0.382709
0
0
0
0
0
1
0.719284
0.879451
0
0
0
0
1
0.688918
0.854616
-0.03453
0
0
0
1
0.688855
0.856167
-0.03328
0.001816
0
0
1
0.688871
0.855869
-0.02561
0.007983
0.008953
0
1
0.689077
0.856053
-0.0262
0.027638
0.024773
0.022965
以及预测误差功率随阶数变化图：
图1、LD算法预测误差功率随阶数的变化
显然从2阶开始预测误差功率就显著减小几乎不变。
用自带函数aryule计算系数，并对比50阶的各系数（横坐标为第几个系数，纵坐标为系数的值）：
图2、预测系数的对比
可以验证LD算法程序是有效的。
二．WV变换
1、问题描述
建立一个非平稳信号，由三个线性调频信号叠加而成：
其中
。
求
的WV分布，指出并分析其信号项和交叉项。
2、WV分布的分析
信号z（t）的Winger-Ville分布定义为：
（1）
由二次叠加原理，Wigner-Ville分布的信号项和交叉项分别为：
（2）
其中，
（3）
由WV分布的定义得：
（4）
在积分公式
（5）
中，令
和
，得：
（6）
同理，有：
（7）
故WV分布的信号项为：
（8）
又由交叉项的定义知：
（9）
利用积分公式，代入相应的参数，得：
（10）
故WV分布的交叉项为：
（11）
式中
（12）
则信号
的WV分布为：
（13）
3、程序设计
closeall;
formatlong;
symsxt;
aifa=10;%设置初始参数
aifa=input（'Pleaseinputaifa:
'）;%ÉèÖÃ
t1=input（'Pleaseinputt1:
'）;
t2=input（'Pleaseinputt2:
'）;
t3=input（'Pleaseinputt3（t1>t2>t3）:
'）;
w1=input（'Pleaseinputw1:
'）;
w2=input（'Pleaseinputw2:
'）;
w3=input（'Pleaseinputw3（w1>w2>w3）:
'）;
%求取z（t+x/2）
z1=exp（-aifa/2*（t+x/2-t1）^2+1j*w1*（t+x/2））+exp（-aifa/2*（t+x/2-t2）^2+1j*w2*（t+x/2））+exp（-aifa/2*（t+x/2-t3）^2+1j*w3*（t+x/2））;
%求取z（t-x/2）的共轭z2=exp（-aifa/2*（t-x/2-t1）^2-1j*w1*（t-x/2））+exp（-aifa/2*（t-x/2-t2）^2-1j*w2*（t-x/2））+exp（-aifa/2*（t-x/2-t3）^2-1j*w3*（t-x/2））;
%为求WV分布信号项作准备
z3=exp（-aifa/2*（t+x/2-t1）^2+1j*w1*（t+x/2））*exp（-aifa/2*（t-x/2-t1）^2-1j*w1*（t-x/2））;
z4=exp（-aifa/2*（t+x/2-t2）^2+1j*w2*（t+x/2））*exp（-aifa/2*（t-x/2-t2）^2-1j*w2*（t-x/2））;
z5=exp（-aifa/2*（t+x/2-t3）^2+1j*w3*（t+x/2））*exp（-aifa/2*（t-x/2-t3）^2-1j*w3*（t-x/2））;
z=z1*z2;
Z_auto=z3+z4+z5;
Wz=simplify（fourier（z）*（（aifa/pi）^（1/4）））
Wz_auto=simplify（fourier（Z_auto）*（（aifa/pi）^（1/4）））
Wz_cross=Wz-Wz_auto
Wz=inline（Wz）
Wz_auto=inline（Wz_auto）
t=0:
0.01:
3;
f=-3.0:
0.01:
3;
[T,F]=meshgrid（t,f）;
W=double（Wz（T,2*pi*F））;
W_auto=double（Wz_auto（T,2*pi*F））;
figure（'name','WV分布'）;
mesh（T,F,W）;
figure（'name','WV分布信号项'）;
mesh（T,F,W_auto）;
W_cross=W-W_auto;
figure（'name','WV分布交叉项'）;
mesh（T,F,W_cross）;
4、程序运行结果
1）值
时，WV分布为
图1、WV分布
图3、WV分布的信号项（左）与交叉项（右）
2.值
时，WV分布为
图4、WV分布
图5、WV分布的信号项（左）与交叉项（右）
升级会员 