现代数字信号处理LevinsonDurbin法和Burg法.docx
《现代数字信号处理LevinsonDurbin法和Burg法.docx》由会员分享，可在线阅读，更多相关《现代数字信号处理LevinsonDurbin法和Burg法.docx（15页珍藏版）》请在冰豆网上搜索。
现代数字信号处理LevinsonDurbin法和Burg法
现代数字信号处理实验报告
Levinson-Durbin算法与Burg算法
Levinson-Durbin算法
一、实验要求:
使用burg法计算x（n）的功率谱:
其中wn是高斯白噪声，阶数选用80阶估计，数据采样分别取128和512。
二、实验原理:
实际应用的随机过程大多数可以用有理传输函数模型很好逼近，输入激励u（n）是均值为零、方差为σ2的白噪声序列，线性系统传输函数为：
其中bk是前馈（或动平均）支路的系数，ak是反馈（或自回归）支路的系数，当除了b0=1之外其余的MA系数都为0，则p阶自回归模型，即AR模型:
此时传输函数为：
传输模型的输出功率谱为：
Wold分解定理认为任何广义平稳随机过程都可分解为一个完全随机的部分和一个确定的部分，如果功率谱完全连续，那么任何ARMA过程或者MA过程都可以用一个无限阶AR过程表示。
而Yule-Walker方程则把模型p阶系数和方差σ2的关系：
也就是只要估计出P+1个自相关函数值，就可以由该式解出P+1个模型参数。
Levinson-Durbin算法的运算数量级为p2.这是一种按阶次进行递推的算法，即首先以AR（0）和AR
（1）模型参数作为初始条件，计算AR
（2），如此类推。
由k阶到k+1阶的计算公式为:
三、编程实现
仿真工具使用MATLAB2015a，根据上述的算法原理，代码可以分为信号采样模块、自相关系数计算模块、初始值模块、算法迭代模块、功率谱计算模块和画图模块，接下来进行依次分析。
信号采样模块如下，其中N是采样数据个数，利用wgn函数直接产生高斯白噪声，wgn（1,N,0）表示随机产生一个强度为0dB的1行N列的高斯白噪声矩阵，矩阵n记录了N个时间值，矩阵x则记录了N个输入信号采样值。
自相关函数计算模块如下，因为阶数为80，所以需要计算R（0）~R（80）的值，此处直接采用xcorr函数进行自相关系数的求解，[b,c]=xcorr（x,p,'biased'）表示计算x序列时间差为[-p,p]的自相关系数值，并依次返回时间差值到矩阵b中，而返回相应相关系数到矩阵c中，取[0,p]所对应的相关系数值，即得相关系数矩阵r。
初始值模块如下，此模块用于模型参数的初始化，a表示AR参数，s表示σ2，d表示D，g表示
。
且σ02=R（0），D0=R
（1），a1,0=1，所以
=D0/σ02=R
（1）/R（0），进而a1,1=-
。
此处需要注意的是，MATLAB矩阵标号按列从1开始，但是σ2、D、a的下标均从0开始。
算法迭代模块如下，首先确定k-1阶的a、σ2、D和
，然后根据原理部分的公式依次计算k阶的σ2、D和
，再根据k阶的
计算k阶的a参数，最后再根据p阶σ2计算p+1阶σ2。
功率谱计算模块如下，在[-pi,pi]范围内的2000个等间隔频率点上均匀采样，直接用linspace函数生成，接着利用书中的式4-27计算功率谱。
以512点估计为例，画图模块如下，此处画出输入信号以及利用Levinson-Durbin算法求得的功率谱，将两幅图显示在一个figure中。
四、仿真结果及分析
由图可知，大致在0.9与0.12附近有峰值，这也与原题目中的0.3pi与0.32pi相对应，由于噪声的存在，使得除了峰值点的其余频率的功率谱不为0，但接近于0，这是因为代码中设置的高斯白噪声的强度为0dB比较小。
另外可以看到N=128和N=512对应的峰值不一样，这也是由于编程方式使得噪声在取了不同的n时噪声信号又重置了。
Burg法
一、实验要求:
使用burg法计算x（n）的功率谱:
其中wn是高斯白噪声，分别使用50阶和80阶估计。
二、实验原理:
在实际应用中，常需根据信号的有限个取样值来估计AR模型的参数，应用较多的是Yule-Walker法或自相关法，协方差法，Burg法。
与自相关法的计算效率很高，而保证预测误差滤波器是最小相位的，但数据两端要附加零取样值，实质上等于数据加窗，这会使参数估计的精度下降，当数据段很短时，加窗效应更严重，是一种直接估计AR参数的方法。
Burg法则一方面希望利用已知数据两端以外的未知数据（但它相对于这些未知数据不做主观臆测），另一方面有设法保证预测误差滤波器是最小相位的，Burg法与自相关法和协方差法不同，它不直接估计AR参数，二是先估计反射系数，然后利用Levinson递推算法由反射系数求得AR参数。
Burg法首先要估计反射系数，所使用的准则是前向和后向预测误差功率估计的平均值最小准则。
在这里，预测误差功率估计仍然用时间平均来代替集合平均。
因此，Burg法估计反射系数的准则为：
前向和后向预测误差滤波器的工作都是在数据段上进行的数据段两端不需要补充零。
一般情况下，求出γk后，即可使用Levinson递推算法中的公式4.25由k-1阶AR参数计算出k阶AR参数。
Burg法的公式可归纳如下：
Burg法估计AR（p）模型参数的具体计算方法如下：
1）确定初始条件
2）确定k-1阶AR参数（迭代计算时，k的值从1开始选取）Ak-1（z），σ2k-1,k-1≤n≤N-1;
3）计算γk;
4）计算Ak（z）;
5）计算
和
，k≤n≤N-1;
6）计算k阶均方误差
7）回到第二部，进行下一次迭代。
一般来说，如果处理的数据来自AR过程，那么采用Burg法可以获得精确地AR谱估计。
三、编程实现:
仿真工具同样使用MATLAB2015a，根据原理，代码也可以分为信号采样模块、初始值模块、算法迭代模块、功率谱计算模块和画图模块，其中信号采样模块、功率谱计算模块以及画图模块与Levinson-Durbin算法代码大致相同，代码中各字母含义也与之前类似，此处不再赘述，下面将着重分析初始值模块、算法迭代模块。
要计算的初始参数包括e0+（n）、e0-（n）、σ02、
，其中e0+（n）和e0-（n）均等于采样的信号矩阵，用矩阵erP和erN表示，σ02则等于采样矩阵的平方均值，根据e0+（n）和e0-（n）可以计算出
，根据
和a1,0=1进一步可计算出a1,1，具体代码如下。
此处同样需要注意的是，MATLAB矩阵标号按列从1开始，但是σ2、a、e0+（n）和e0-（n）的下标均从0开始。
算法迭代模块如下，首先确定k-1阶的a、σ2、ek-1+（n）、ek-1-（n），然后根据k-1阶的ek-1+（n）、ek-1-（n）计算k阶的ek+（n）、ek-（n）和
，再根据k阶的
计算k阶的a参数和σ2，最后再根据p阶σ2计算p+1阶σ2。
四、仿真结果及分析
由图可知，大致在0.9与0.12附近有峰值，但由于相离较近难以辨别，这也与原题目中的0.3pi与0.32pi相一致，一高一低现象是随机造成的，因为单次的序列进行计算会有随机性，但总体是会趋向于一致。
由于是按照书上4.8.3节标准公式所写，放大可以发现，AR谱估计中存在两个靠的很近的谱峰，即有谱线分裂现象，Burg算法得到的谱估计，其谱峰位置移动有可能达到原位置的16%，如果要进行修正，那么可以选择相应的算法进行改进，例如复数据和反射系数修正等等。
而AR模型阶数的确定，可以用最终预测误差（FPE）准则或者AIC信息准则来确定，阶数选的太低，功率谱显得平滑，阶数太高，会出现虚假峰，因此需要选择合适的阶数，此试验中阶数已经确定了。
附录（完整代码）
Levinson-Durbin算法：
%%清除工作空间和变量空间
clear;
closeall;
clc;
%%信号采样
%采样数据个数
N=512;
n=0:
（N-1）;
%高斯白噪声
w=wgn（1,N,0）;
%输入信号x（n）
x=cos（0.3*pi*n）+cos（0.32*pi*n）+w;
%%计算自相关系数矩阵，R（0）-R（80）
p=80;
[b,c]=xcorr（x,p,'biased'）;
r=[];
fori=1:
length（b）
ifc（i）>=0
r=[rb（i）];
end
end
%%算法迭代
a=zeros（p,p+1）;%a参数
s=zeros（1,p+1）;%σ2
d=zeros（1,p）;%Dk
g=zeros（1,p）;%
%初始化参数
a（1,1）=1;
s
（1）=r
（1）;
d
（1）=r
（2）;
g
（1）=r
（2）/r
（1）;
a（1,2）=-g
（1）;
%迭代
fori=2:
p
s（i）=（1-g（i-1）^2）*s（i-1）;
d（i）=0;
fork=1:
i
d（i）=d（i）+a（i-1,k）*r（i+2-k）;
end
g（i）=d（i）/s（i）;
a（i,1）=1;
form=1:
（i-1）
a（i,m+1）=a（i-1,m+1）-g（i）*a（i-1,i-m+1）;
end
a（i,i+1）=-g（i）;
end
s（p+1）=（1-g（p）^2）*s（p）;
%%功率谱
W=linspace（-pi,pi,2000）;
S=zeros（1,length（W））;
fori=1:
2000
spot=W（i）;
Ssum=1;
fork=1:
p
Ssum=Ssum+a（p,k+1）*exp（-j*spot*k）;
end
S（i）=s（p+1）/（abs（Ssum）^2）;
End
%%画图模块
subplot（1,2,1）;
plot（n,x）;
title（'512点信号'）;
subplot（1,2,2）;
plot（W,S）;
title（'512点80阶Levinson估计'）;
Burg法：
%%清除工作空间和变量空间
clear;
closeall;
clc;
%%信号产生
N=512;
n=0:
（N-1）;
w=wgn（1,N,0）;
x=cos（0.3*pi*n）+cos（0.32*pi*n）+w;
%%初始化参数
p=80;%设置阶数
erP=zeros（p+1,N）;%ek+（n）
erN=zeros（p+1,N）;%ek-（n）
erP（1,:
）=x;
erN（1,:
）=x;
g=zeros（1,p）;%
s=zeros（1,p+1）;%σ2
fori=1:
N
s
（1）=s
（1）+x（i）^2;
end
s
（1）=s
（1）/N;
a=zeros（p,p+1）;
a（1,1）=1;
temp1=0;
temp2=0;
fori=2:
N
temp1=temp1+2*x（i）*x（i-1）;
temp2=temp2+x（i）^2+x（i-1）^2;
end
g
（1）=temp1/temp2;
a（1,2）=-g
（1）;
%%迭代
fori=2:
p
erP（i,1）=x
（1）;
erN（i,1）=x
（1）;
fork=2:
N
erP（i,k）=erP（i-1,k）-g（i-1）*erN（i-1,k-1）;
erN（i,k）=erN（i-1,k-1）-g（i-1）*erP（i-1,k）;
end
s（i）=（1-g（i-1）^2）*s（i-1）;
temp1=0;temp2=0;
form=（i+1）:
N
temp1=temp1+erP（i,m）^2+erN（i,m-1）^2;
temp2=temp2+2*erP（i,m）*erN（i,m-1）;
end
g（i）=temp2/temp1;
a（i,1）=1;
form=2:
i
a（i,m）=a（i-1,m）-g（i）*a（i-1,i-m+2）;
end
a（i,i+1）=-g（i）;
end
s（p+1）=（1-g（p）^2）*s（p）;
%%功率谱
W=linspace（-pi,pi,2000）;
S=zeros（1,length（W））;
fori=1:
2000
spot=W（i）;
Ssum=1;
fork=1:
p
Ssum=Ssum+a（p,k+1）*exp（-j*spot*k）;
end
S（i）=s（p+1）/（abs（Ssum）^2）;
End
%%画图模块
subplot（1,2,1）;
plot（n,x）;
title（'信号'）;
subplot（1,2,2）;
plot（W,S）;
title（'80阶Burg法估计'）;
升级会员 