现代信号处理的大作业LD算法以及WV变换.docx
- 文档编号:7051077
- 上传时间:2023-01-16
- 格式:DOCX
- 页数:14
- 大小:248.36KB
现代信号处理的大作业LD算法以及WV变换.docx
《现代信号处理的大作业LD算法以及WV变换.docx》由会员分享,可在线阅读,更多相关《现代信号处理的大作业LD算法以及WV变换.docx(14页珍藏版)》请在冰豆网上搜索。
现代信号处理的大作业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分布的信号项(左)与交叉项(右)
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 现代 信号 处理 作业 LD 算法 以及 WV 变换
![提示](https://static.bdocx.com/images/bang_tan.gif)