现代信号处理的大作业：LD算法以及WV变换

1现代信号处理大作业姓名：学号：专业：电子科学与技术2目录现代信号处理大作业....................................................................................................1一.L-D算法的仿真实现............................................................................................31、问题描述.........................................................................................................32、算法分析.........................................................................................................32.1、LD算法原理.......................................................................................32.2、LD算法的实现...................................................................................43、程序实现思路.................................................................................................54、程序如下：.....................................................................................................54.1、主函数：..............................................................................................54.2、LD算法子函数：.................................................................................74.3、仿真结果：..........................................................................................8二．WV变换...............................................................................................................101、问题描述.......................................................................................................102、WV分布的分析.............................................................................................103、程序设计.......................................................................................................124、程序运行结果...............................................................................................133一.L-D算法的仿真实现1、问题描述用Matlab实现Levinsion-Durbin算法。2、算法分析2.1、LD算法原理由于语音样点之间存在相关性，所以可以用过去的样点值来预测现在或未来的样点值。如下图所示(),(1),...,(1)xnpxnpxn图1线性预测图示由上图可得plpllnxanx1)()(，从而可以通过使实际语音x（n）和线性预测结果)(nx之间的误差e（n）在某个准则下达到最小值来决定唯一的一组预测系数pla。而这组系数就能反映语音信号的特性，可以作为语音信号特征参数来用于语音编码、语音合成和语音识别等应用中去。由估计值和实际信号值的误差,可有010()()()()()(),1ppplplpllenxnxnxnaxnlaxnla根据e(n)最小均方误差准则，来决定唯一的一组预测系数pla，即:min)()(eE122plpllnxanxEn,由此可得到Y-W方程：线性预测4pkklkRalpl,...,2,100minp0，取遍k值之后有以下：ppppaaaRpRpRpRRRpRRR...0...1............1...01...100...010min由相关函数的偶函数性质有：ppppaaaRpRpRpRRRpRRR...0...1............1-...01...100...010min在已知自相关函数的前提下，根据e(n)均方误差最小的原则来求解a，本实验中采用Levinson-Durbin算法。2.2、LD算法的实现Levinson-Durbin算法首先由一阶AR模型开始，按照前面的Y-W方程可有，一阶AR模型（p=1）的Y-W方程是010101011021121112111xxxxxxxxxxxxxxrarraarrrr该方程解出：然后增加一阶，即令p=2，可得到：00a1012101210222221arrrrrrrrrxxxxxxxxxxxxxxxxxx由上式可解出：2122222112211212111221/22aaaaararaxxxx5然后令,...3,2p以此类推，可以得到一般的递推公式：11121(1)(1)()2221220101,2,...1pxxpxxkppppppkpkpppkpppxxrparpkkkaaakakrExnkp式中的pk称为反射系数，而2p和21p是预测误差的均方误差值，因此21pk必须大于等于0，这样pk应满足1pk，进而得到221pp，即预测误差随递推次数增加而减少。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;%数据长度6noise=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];%预测出来的系数%预测误差功率图figureplot(1:m,P(1:m),'rh')ylim([0,6])title('预测误差功率随迭代次数的变化');xlabel('迭代次数')ylabel('预测误差功率')%利用aryule验证y=aryule(x,m);figuresubplot(1,2,1)plot(1:51,aa(m,:),'r')title('编写程序运行得出的系数')subplot(1,2,2)plot(1:51,y,'k')title('aryule运行得出的系数')74.2、LD算法子函数：%%Levinson-Durbinalgorithmfunction[a,P]=Levinson_Durbin(x,m)%x为输入信号%m为阶数a=zeros(m,m);N=length(x);P=zeros(1,m+1);fori=1:NP(1)=P(1)+x(i)^2;endP(1)=P(1)/N;%计算预测误差功率的初始值P_0f=zeros(m+1,N);%m阶前向预测误差g=zeros(m+1,N);%m阶后向预测误差f(1,:)=x;g(1,:)=x;K=zeros(1,m);%开始迭代：fori=1:mnumerator=0;denominator=0;forj=i+1:Nnumerator=numerator-f(i,j)*g(i,j-1);denominator=denominator+(f(i,j)^2+g(i,j-1)^2)/2;endK(i)=numerator/denominator;%反射系数forj=1:i-18a(i,j)=a(i-1,j)+K(i)*a(i-1,i-j);enda(i,i)=K(i);ifi=2P(i)=(1-K(i)^2)*P(i-1);%更新预测误差功率endforj=i+1:Nf(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);endend4.3、仿真结果：显然，我们想要预测估计的是a0=[10.720.88]，预测结果为（预测了1-50阶，后面阶数过多不易展示，因此只列出前1-6阶预测结果）：表1：LD算法预测系数表（1-6阶）10.3827090000010.7192840.879451000010.6889180.854616-0.0345300010.6888550.856167-0.033280.0018160010.6888710.855869-0.025610.0079830.008953010.6890770.856053-0.02620.0276380.0247730.022965以及预测误差功率随阶数变化图：9图1、LD算法预测误差功率随阶数的变化显然从2阶开始预测误差功率就显著减小几乎不变。用自带函数aryule计算系数，并对比50阶的各系数（横坐标为第几个系数，纵坐标为系数的值）：图2、预测系数的对比可以验证LD算法程序是有效的。10二．WV变换1、问题描述建立一个非平稳信号，由三个线性调频信号叠加而成：12224112233()()[exp