现代数字信号处理-Levinson-Durbin法和Burg法

现代数字信号处理实验报告Levinson-Durbin算法与Burg算法目录Levinson-Durbin算法................................................................................................1一、实验要求:........................................................................................................1二、实验原理:........................................................................................................1三、编程实现.........................................................................................................2四、仿真结果及分析..............................................................................................3Burg法......................................................................................................................5一、实验要求:........................................................................................................5二、实验原理:........................................................................................................5三、编程实现:........................................................................................................6四、仿真结果及分析..............................................................................................7附录（完整代码）.....................................................................................................9Levinson-Durbin算法：.........................................................................................9Burg法：.............................................................................................................101Levinson-Durbin算法一、实验要求:使用burg法计算x(n)的功率谱:𝑥𝑛=cos(0.3𝜋𝑛)+cos(0.32𝜋𝑛)+𝑤(𝑛)其中wn是高斯白噪声，阶数选用80阶估计，数据采样分别取128和512。二、实验原理:实际应用的随机过程大多数可以用有理传输函数模型很好逼近，输入激励u(n)是均值为零、方差为σ2的白噪声序列，线性系统传输函数为：𝐻(𝑧)=𝐵(𝑧)𝐴(𝑧)=∑𝑏𝑘𝑧−𝑘𝑞𝑘=0∑𝑎𝑘𝑧−𝑘𝑝𝑘=0其中bk是前馈（或动平均）支路的系数，ak是反馈（或自回归）支路的系数，当除了b0=1之外其余的MA系数都为0，则p阶自回归模型，即AR模型:𝑥(𝑛)=−∑𝑎𝑘𝑥(𝑛−𝑘)+𝑢(𝑛)𝑝𝑘=1此时传输函数为：𝐻(𝑧)=1𝐴(𝑧)=11+∑𝑎𝑘𝑧−𝑘𝑝𝑘=1传输模型的输出功率谱为：𝑆𝑥𝑥(𝑒𝑗𝜔)=𝜎2|A(𝑒𝑗𝜔)|2=𝜎2|1+∑𝑎𝑘𝑒−𝑗𝜔𝑘𝑝𝑘=1|2Wold分解定理认为任何广义平稳随机过程都可分解为一个完全随机的部分和一个确定的部分，如果功率谱完全连续，那么任何ARMA过程或者MA过程都可以用一个无限阶AR过程表示。而Yule-Walker方程则把模型p阶系数和方差σ2的关系：R𝑥𝑥(𝑚)={−∑𝑎𝑘𝑅𝑥𝑥(𝑚−𝑘)+𝜎2,𝑚=0𝑝𝑘=1−∑𝑎𝑘𝑅𝑥𝑥(𝑚−𝑘),𝑚0𝑝𝑘=1也就是只要估计出P+1个自相关函数值，就可以由该式解出P+1个模型参数。Levinson-Durbin算法的运算数量级为p2.这是一种按阶次进行递推的算法，即首先以AR(0)和AR(1)模型参数作为初始条件，计算AR(2)，如此类推。由k阶到k+1阶的计算公式为:𝜎𝑘2=𝑅(0)+∑𝑎𝑘,𝑖𝑅(𝑖)𝑘𝑖=0𝐷𝑘=∑𝑎𝑘,𝑖𝑅(𝑘+1−𝑖),𝑎𝑘,0=1𝑘𝑖=0𝛾𝑘+1=𝐷𝑘𝜎𝑘22𝜎𝑘+12=(1−𝛾𝑘+12)𝜎𝑘2𝑎𝑘+1,𝑖=𝑎𝑘,𝑖−𝛾𝑘+1𝑎𝑘,𝑘+1−𝑖,𝑖=1,2…,𝑘;𝑎𝑘+1,𝑘+1=−𝛾𝑘+1三、编程实现仿真工具使用MATLAB2015a，根据上述的算法原理，代码可以分为信号采样模块、自相关系数计算模块、初始值模块、算法迭代模块、功率谱计算模块和画图模块，接下来进行依次分析。信号采样模块如下，其中N是采样数据个数，利用wgn函数直接产生高斯白噪声，wgn(1,N,0)表示随机产生一个强度为0dB的1行N列的高斯白噪声矩阵，矩阵n记录了N个时间值，矩阵x则记录了N个输入信号采样值。0:1;1,,0;0.3**0.32**;nNwwgnNxcospincospinw自相关函数计算模块如下，因为阶数为80，所以需要计算R(0)~R(80)的值，此处直接采用xcorr函数进行自相关系数的求解，[b,c]=xcorr(x,p,'biased')表示计算x序列时间差为[-p,p]的自相关系数值，并依次返回时间差值到矩阵b中，而返回相应相关系数到矩阵c中，取[0,p]所对应的相关系数值，即得相关系数矩阵r。[80;,,,'';;1:    0        ;    ]pbcxcorrxpbiasedrforilengthbifcirrbiendend初始值模块如下，此模块用于模型参数的初始化，a表示AR参数，s表示σ2，d表示D，g表示γ。且σ02=R(0)，D0=R(1)，a1,0=1，所以γ1=D0/σ02=R(1)/R(0)，进而a1,1=-γ1。此处需要注意的是，MATLAB矩阵标号按列从1开始，但是σ2、D、a的下标均从0开始。,1;1,1;1,;1,;1,11;11;12;12/1;1,21;azerosppszerospdzerospgzerospasrdrgrrag算法迭代模块如下，首先确定k-1阶的a、σ2、D和γ，然后根据原理部分的公式依次计算k阶的σ2、D和γ，再根据k阶的γ计算k阶的a参数，最后再根据p阶σ2计算p+1阶σ2。322:()(1(1))*(1);0;1:      1,*2;   /;   ,11;      1:1      ,11,1*1,1;   ,foripsigisidiforkididiaikrikendgidisiaiformiaimaimgiaiimendaii21;(1)(1())*();giendspgpsp功率谱计算模块如下，在[-pi,pi]范围内的2000个等间隔频率点上均匀采样，直接用linspace函数生成，接着利用书中的式4-27计算功率谱。2,,2000;1,;1:2000();1;1:(,1)*(**);()(1)/(());WlinspacepipiSzeroslengthWforispotWiSsumforkpSsumSsumapkexpjspotkendSispabsSsumend以512点估计为例，画图模块如下，此处画出输入信号以及利用Levinson-Durbin算法求得的功率谱，将两幅图显示在一个figure中。1,2,1;,;'512';1,2,2;,;'51280'();()subplotplotnxtitlesubplotplotWStitleLevinson点信号点阶估计四、仿真结果及分析4由图可知，大致在0.9与0.12附近有峰值，这也与原题目中的0.3pi与0.32pi相对应，由于噪声的存在，使得除了峰值点的其余频率的功率谱不为0，但接近于0，这是因为代码中设置的高斯白噪声的强度为0dB比较小。另外可以看到N=128和N=512对应的峰值不一样，这也是由于编程方式使得噪声在取了不同的n时噪声信号又重置了。5Burg法一、实验要求:使用burg法计算x(n)的功率谱:𝑥𝑛=cos(0.3𝜋𝑛)+cos(0.32𝜋𝑛)+𝑤(𝑛)其中wn是高斯白噪声，分别使用50阶和80阶估计。二、实验原理:在实际应用中，常需根据信号的有限个取样值来估计AR模型的参数，应用较多的是Yule-Walker法或自相关法，协方差法，Burg法。与自相关法的计算效率很高，而保证预测误差滤波器是最小相位的，但数据两端要附加零取样值，实质上等于数据加窗，这会使参数估计的精度下降，当数据段很短时，加窗效应更严重，是一种直接估计AR参数的方法。Burg法则一方面希望利用已知数据两端以外的未知数据（但它相对于这些未知数据不做主观臆测），另一方面有设法保证预测误差滤波器是最小相位的，Burg法与自相关法和协方差法不同，它不直接估计AR参数，二是先估计反射系数，然后利用Levinson递推算法由反射系数求得AR参数。Burg法首先要估计反射系数，所使用的准则是前向和后向预测误差功率估计的平均值最小准则。在这里，预测误差功率估计仍然用时间平均来代替集合平均。因此，Burg法估计反射系数的准则为：ε=∑[𝑒𝑝+(𝑛)2+𝑒𝑝−(𝑛)2]=𝑚𝑖𝑛𝑁−1𝑛=𝑝前向和后向预测误差滤波器的工作都是在数据段上进行的数据段两端不需要补充零。一般情况下，求出γk后，即可使用Levinson递推算法中的公式4.25由k-1阶AR参数计算出k阶AR参数。Burg法的公式可归纳如下：𝛾𝑘=2∑𝑒𝑘+1+𝑁−1𝑛=𝑘(𝑛)𝑒𝑘−1−(𝑛−1)∑[(𝑒𝑘−1+(𝑛))2+(𝑒𝑘−1−(𝑛−1))2]𝑁−1𝑛=𝑘[𝑒𝑘+(𝑛)𝑒𝑘−(𝑛)]=[1−𝛾𝑘−𝛾𝑘1][𝑒𝑘−1+(𝑛)𝑒𝑘−1−(𝑛−1)]𝑎𝑘,𝑖=𝑎𝑘−1,𝑖−𝛾𝑘𝑎𝑘−1,𝑘−𝑖,𝑖=1,2,3,…,𝑘,𝑎𝑘−1,0=1.Burg法估计AR（p）模型参数的具体计算方法如下：1）确定初始条件𝑒0+(𝑛)=𝑒0−(𝑛)=𝑥(𝑛),σ02=1�