现代数字信号处理LevinsonDurbin法和Burg法文档格式.docx
- 文档编号:19399039
- 上传时间:2023-01-05
- 格式:DOCX
- 页数:15
- 大小:232.54KB
现代数字信号处理LevinsonDurbin法和Burg法文档格式.docx
《现代数字信号处理LevinsonDurbin法和Burg法文档格式.docx》由会员分享,可在线阅读,更多相关《现代数字信号处理LevinsonDurbin法和Burg法文档格式.docx(15页珍藏版)》请在冰豆网上搜索。
由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法
其中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阶的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,'
);
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);
a(i,i+1)=-g(i);
s(p+1)=(1-g(p)^2)*s(p);
%%功率谱
W=linspace(-pi,pi,2000);
S=zeros(1,length(W));
2000
spot=W(i);
Ssum=1;
p
Ssum=Ssum+a(p,k+1)*exp(-j*spot*k);
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);
512点80阶Levinson估计'
Burg法:
%%信号产生
%%初始化参数
%设置阶数
erP=zeros(p+1,N);
%ek+(n)
erN=zeros(p+1,N);
%ek-(n)
erP(1,:
)=x;
erN(1,:
%
%σ2
N
s
(1)=s
(1)+x(i)^2;
s
(1)=s
(1)/N;
temp1=0;
temp2=0;
temp1=temp1+2*x(i)*x(i-1);
temp2=temp2+x(i)^2+x(i-1)^2;
g
(1)=temp1/temp2;
%%迭代
erP(i,1)=x
(1);
erN(i,1)=x
(1);
fork=2:
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);
s(i)=(1-g(i-1)^2)*s(i-1);
temp1=0;
temp2=0;
form=(i+1):
temp1=temp1+erP(i,m)^2+erN(i,m-1)^2;
temp2=temp2+2*erP(i,m)*erN(i,m-1);
g(i)=temp2/temp1;
form=2:
a(i,m)=a(i-1,m)-g(i)*a(i-1,i-m+2);
信号'
80阶Burg法估计'
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 现代 数字信号 处理 LevinsonDurbin Burg