Levinson-Durbin算法实验报告 2

1、Levinson-Durbin算法一、实验目的学会LEVINSON-DURBIN算法求解Yule-Walker方程，并通过MATLAB实现程序及仿真。二、实验原理1、线性预测分析的基本原理由于语音样点之间存在相关性，所以可以用过去的样点值来预测现在或未来的样点值。如下图所示线性预测图1 线性预测图示由上图可得，从而可以通过使实际语音x（n）和线性预测结果之间的误差e（n）在某个准则下达到最小值来决定唯一的一组预测系数。而这组系数就能反映语音信号的特性，可以作为语音信号特征参数来用于语音编码、语音合成和语音识别等应用中去。2、 Yule-Walker（Y-W）方程由估计值和实际信号值的误差可有

2、根据e(n)最小均方误差准则，来决定唯一的一组预测系数，即由此可得到Y-W方程：取遍k值之后可有以下：由相关函数的偶函数性质可有在已知自相关函数的前提下，根据e(n)均方误差最小的原则来求解ak，本实验中采用Levinson-Durbin算法。3、 Levinson-Durbin算法简介Levinson-Durbin算法首先由一阶AR模型开始，按照前面的Y-W方程可有，一阶AR模型（p=1）的Y-W方程是然后增加一阶，即令p=2，可得到：由上式可解出：然后令以此类推，可以得到一般的递推公式：式中的称为反射系数，而3、 实验过程实验开始先按照前文所述的Levinson算法的递推公式编写子函数fu

3、nction a_p,E_V=Levinson_algorithm(x,p)其中a_p是输出的预测系数，E_V是输出的预测误差均放值。该子函数实现用Levinson-Durbin算法求解Yule-Walker方程，用迭代方法求解与谱估计有关的参量。利用子函数估计AR（2）模型中的参数值，估计其功率谱，并与已知的信号功率谱进行比较。参数设置：已知信号为白噪声通过线性系统H(z)=1/(1+a1*z(-1)+a2*z(-2)产生信号向量，p=90，先假设已知模型的输入系数为a0=1 0.78 0.92。实验结果由图可知，两图的频谱基本一致，即可知该子函数可以很好的估计出模型参数。附录：函数程序%子

4、函数程序%Levinson-Durbin 函数%参量：A(R)模型的输入参数为信号x和滤波器阶数p，输出参数为模型预测系数a_p和预 测功率误差E_p；function a_p,E_V=Levinson_algorithm(x,p) %输入信号的互相关函数计算 N=length(x); for i=1:N Rx(i)=x(1:N-i+1)*(x(i:N)/N; %计算输入信号的自相关函数值 end%levinson算法迭代过程 a(1)=1; a(2)=-Rx(2)/Rx(1); %自适应滤波器的前两个系数 for k=1:p-1 E_V(k+1)=Rx(1)+a(2:k+1)*Rx(2:k+

5、1); %预测误差更新 D(k+1)=a(1:k+1)*(fliplr(Rx(2:k+2); %为扩大方程矩阵中的Dk的更新 gama(k+2)=D(k+1)/E_V(k+1); %反射系数gama的更新 E_V(k+2)=E_V(k+1)*(1-(gama(k+2)2); %预测误差的更新%预测系数的更新 a_yuce(1)=1; %由k阶系数预测第k+1阶系数 for q=1:k a_yuce(q+1)= a(q+1)-gama(k+2)* a(k-q+2); end a_yuce(k+2)=-gama(k+2); a=a_yuce; %准备下一次迭代 end a_p=a; %算法迭代完成

6、输出模型系数 E_p=E_V; %输出系统预测误差%函数应用实例%运用前面所编写的levinson算法程序计算已知信号的预测误差,一次检验所设计的算法准确性。%给出一个已知的模型，让一直信号经过该模型之后利用函数估计模型系数%已知模型设为x(n)=a1*u(n-2)+a2*u(n-1)+a3*u(n)clear;close all;clc;%产生已知的信号u=randn(1,3000); %产生均值为0，方差（功率）为1，数据长度为3000的高斯白噪声a0=1 0.72 0.88; %已知模型的输入系数x=filter(1,a0,u); %产生的信号模型是x(n)=a0(1)*u(n-2)+a

7、0(2)*u(n-1)+a0(3)*u(n); %其中的a0(i)为a0中的数值%由已知模型画出信号的功率谱函数图w=linspace(-pi,pi,3000); %将角频率w从-pi到pi平均分成2000等份for m=1:3000 c=w(m); s(m)=1/(abs(1+a0(2:3)*exp(-j*c*(1:2)2; %信号功率谱s，其中系数为给定的模型参数endfigure(1);plot(w,s); %画出一直信号的频谱图title(根据已知的模型参数求得信号的频谱); %下面利用所编写的算法求得预测模型系数p=50; %用于估计的模型参数a_p,E_V=Levinson_algorithm(x,p); %调用子函数for m=1:3000 c=w(m); sx(m)=1/(abs(1+a_p