郑州大学随机信号处理大作业

包括levense递推关系式和burg算法

基于AR模型的功率谱估计

摘要

频域分析又称谱分析，主要研究信号在中的各种特征。功率谱密度函数描述随机过程的功率随频率的平均分布。功率谱估计问题就是根据一组有限观测值来确定该过程谱的内容。对于平稳随机过程来说，功率谱理论上的数值不可能实现，只能用有限的观测值来估计或者逼近真实值，从而叫真实的反应问题的本质。当然估计结果的好坏，与采用的处理方法有很大关系。

本文概要介绍了古典谱估计方法以及其基本原理，详细介绍了基于AR模型的功率谱估计的原理、算法。古典谱估计主要基于离散傅里叶变换（DFT），包括周期图法和相关函数法及在此基础上改进的方法；基于AR模型的功率谱估计，是参数法做功率谱估计的一种，简而言之，就是通过计算模型的参数而不是做FFT来做功率谱估，本文主要介绍的基于AR模型的功率谱估计的方法，主要有Levinson-Durbin快速算法和Burg算法。本文通过MATLAB编程分别利用古典法和基于AR模型的功率谱估计法对信号加噪声进行功率谱估计，并通过不同的点数以及不同的频率间隔的比较来分析古典法和基于AR模型的功率谱估计法做出相应的分析。

关键词：AR模型 功率谱估计 Levinson-Durbin快速算法 Burg算法

第1章 引言

谱估计在信号处理方法中有很重要的作用，它不仅是分析、了解信号所含有用信息的工具，也是信号内在本质的一种表现形式。古典谱估计方法存在着分辨率低和方差性能不好的问题。现代谱分析是针对经典谱估计的不足，近几年提出的谱估计新方法，其基本思路是:在进行谱估计过程中，对有限数据以外的数据不做任何确定性假设，利用己知数据，采用外推或预测方法，由观测数据推求m N以后的数据。在这种情况下，首先必须寻找一个产生随机信号的模型，通过观测的数据把模型的参数估计出来，进而求得所需的功率谱。基于模型的功率谱估计方法，不作自相关函数Rxx(m) 0, m > n的假设，从而避免了功率泄漏，提高了分辨率。

功率谱估计应用范围很广，日益受到各学科和应用领域的极大重视。以傅 立叶变换为基础的传统(或经典)谱估计方法，虽然具有计算效率高的优点，但 却有着频率分辨率低和旁瓣泄漏严重的固有缺点。这就迫使人们大力研究现代 谱估计方法。现代谱估计方法是以参数模型为基础的方法。

包括levense递推关系式和burg算法

第2章 谱估计理论基础

2.1 几种主要的谱估计方法

信号处理的核心，说到底就是如何保证在信号受到干扰产生失真的情况下， 正确恢复原有信号，提取有用信息。而功率谱(简称谱)估计就是信号处理的一个重要分支，它应用范围很广，日益受到各学科和应用领域的极大重视。它是在频率域研究随机信号的统计规律，其中心目的为了估计出随机干扰信号并将其去除，提取有用信号，对语音、声纳、雷达等信号处理，有着重要的意义，广泛应用于通信、控制、地球物理，它的研究对象主要是零均值平稳高斯过程。以傅立叶变换为基础谱估计一般称为的传统(或经典)谱估计方法，传统谱估计法又可以分为直接法和间接法，后来由于FFT的出现，直接法和间接法往往被结合起来使用。在信号分析方法中，傅立叶变换是较常用的数学工具。时间信号经过傅立叶变换后，可得到它的频谱，平稳随机过程的相关函数和它的功率谱密度是一对傅立叶变换对。

近几年，在信号功率谱密度估计方面出现了许多新的算法，其中应用最广泛的算法是1967年由Burg提出的最大嫡谱估计，这些新方法连同演变出来的各种算法不下几十种，统称为现代谱估计。它是相对经典谱估计方法而言的。其比较有名的是:Levinson-Durbin算法自回归模型、Burg算法的最大嫡分析、正反向线性预测的最小二乘算法、自回归模型、滑动平均模型、自回归一滑动平均模型Pisarenko谐波分解法、Prony提取极点法、Prony谱线分解法以及Capon最大似然估计法。根据算法的基本属性，把这些算法归纳在图2. l

2.2 傅里叶变换

傅立叶变换(DFT)认为:有限长的数据段可看作无限长的取样序列进行加窗 截断后的结果。不论是数据加窗还是自相关函数加窗，在频率域都会发生“泄露”现象，即功率谱主瓣的能量泄露到旁瓣中去，这样，弱信号的主瓣很容易被强信号的旁瓣淹没或畸变，造成谱的模糊与失真。为了降低旁瓣，很多学者在选择窗口函数的形式上和窗口处理函数方法上想办法，但是所有的旁瓣抑制

包括levense递推关系式和burg算法

技术都是以损失谱分辨率为代价的。谱分析应用中，谱分辨率和低旁瓣是同样重要的指标，在作DFT时，人们常在有效数据后面补一些零，使得补零后的数据N为2的整数次幂以便于快速傅立叶变换，而那些认为在傅立叶变换之前对数据段补零可以改善周期图的谱分辨率是一种模糊的错误概念，因为补零后，虽然谱线能够增密，但因为补零没有增加任何新的信息，因此不可能提高分辨率。谱分辨率的极限只能由取样数据段的长度决定。

2.3 AR模型

2.3.1 引言

为了克服经典谱估计的缺点，近年来在实现高分辨率谱估计技术方面取得了很大的进展，提出了许多功率谱估计的参数方法，也就是现代谱估计的基本方法。其基本思想是在进行谱估计过程对所观测的有限数据以外的数据不作任何确定性假设。在对这些数据如何产生具有一定先验知识的基础上，采取外推或预测的方法，从观测的数据推求m大于等于N以后的数据。显然，在这种情况下，首先必须寻求一个产生随机信号的模型。通过观测数据把模型的参数估计出来，进而求得所需的功率谱。所以参数法是基于模型(信号参数模型)的功率谱估计方法。这种方法既不需要加窗、又不对相关函数做Rxx(m) 0,m N的假设，从而避免功率泄漏，提高了分辨率。

本节所讨论的信号模型，认为所观测的数据x (n)是由某输入序列u (n)激励线性离散物理可实现系统的相应输出，如图2. 2所示，按离散系统的系统函数其通式可写成为

相应的差分方程为

在功率谱估计中，若观测的数据x (n)是平稳随机过程，则该系统的输入u (n) 也可以认为是平稳的，因而观测数据的功率谱为

一般信号的模型的输入数据是不能观测到的，在已知模型参数或频响|H(e

jw)|

包括levense递推关系式和burg算法

的条件下，为了求得模型输出功率谱的估值，最简单的方法就是假设输入u (n)为零均值的白噪序列。离散时间白噪声自相关函数与功率谱分别为

故最后求得的观测数据的功率谱为

由此可见，基于模型的现代谱估计过程一般可以分为三个步骤进行:

(1)为被估计的随机过程确定或选择一个合理的模型。这有赖于对随机过程进行 理论分析和实验研究。

(2)根据已知观测数据x(n),0 n N 1，估计模型的参数{ak}和{{bk )。这涉及到 对各种算法的研究。通常，模型参数的数据量比观测数据的数据量少很多，因此，为数据压缩创造了条件。

(3)用估计得到的模型参数代入式(2. 3. 4)计算功率谱的估值。

由式(2. 3. 1)可见，产生随机序列的模型其系统函数是有理函数，具有q个零点 和p个极点称为自回归滑动平均模型，用ARMA (p, q)表示。若q=0,b0 1 ，系统函数H (z) =1/A (z)只有极点称为自回归((AR)模型，用AR (p)表示。若A (z) =1，则系统函数H(z)=B(z)只有零点称为滑动平均(MA)模型，用MA (q)表示。在这三种模型中，AR模型得到普遍应用，其原因是AR模型的参数计算是线性方程比较简单，它与建立在外推自相关函数时保持原概论空间的嫡最大的最大嫡法是等价的。同时很适于表示很窄的频谱，在谱估计时，由于递推特性所以所需的数据很短，而MA模型表示窄谱时一般需要数量很多的参数，ARMA模型虽然所需的参数数量最少，但参数估计的算法是非线性方程组，其运算远比AR模型复杂。况且根据Wold证明，任意ARMA或MA信号模型可以用无限阶或阶数足够大的AR模型来表示的.所以本论文对现代谱估计方法着重介绍利用AR模型进行功率谱估计的理论和算法。

2.3.2 AR模型的Yule-walker方程

以AR模型为基础的谱估计可由式(2. 3. 4)计算，只是这里的q=0,b0 0 m。这就需要知道模型的阶数P和P个AR系数，以及模型激励源的方差 2。为此，必须把这些参数和己知(或估计得到)的自相关函数联系起来，这就是著Yule-walker方程。

包括levense递推关系式和burg算法

2.3.3 Levinson-Durbin快速递推法

用线性方程组的常用算法(例如高斯消元法)求解(2. 3. 9)式，需要的运算量数量级为p3。但若利用系数矩阵的对称性和Toeplitz性质，则可构成一些高效算法，Levinson-Durbin算法是其中最著名、应用最广泛的一种。这种算法的运算量量级为p2。这是一种按阶次进行递推的算法，即首先以AR (0)和AR (1)模型参数作

包括levense递推关系式和burg算法

为初始条件，计算AR (2)模型参数;然后根据这些参数计算AR (3)模型的参数，等等，一直到计算出AR (p)模型参数为止。这样，当整个迭代计算结束后，不仅求得了所需要的p阶AR模型的参数，而且还得到了所有各低阶模型的参数。定义

am(k)为p阶AR模型在阶次为m时的第k个系数，k =1,2...m，m =1, 2... 。 m

为m阶时的前向预测的最小误差功率(此处省去了“min"，且 m m2)。由(2. 3. 9)式，当m =1时，有

包括levense递推关系式和burg算法

第3章 AR模型参数提取方法

3.1.1 Yule-Walker算法

用最小平方时间平均准则代替集合平均准则，有

式中，e p(n)可由长度为P+1的预测误差滤波器冲激响应序列

(1,ap1,ap2, ,app)与长度为N的数据序列((x(0), x(1) ,...}x(N-1))进行卷积得到。因而e p(n)序列的长度为N+P，这就决定了式(3. 1. 1)中求和的项数。显然，在计算卷积时，在数据段xN(n)的两端，实际上添加了若干零取样值。实际上，e p(n)是由xN(n)经过冲激响应为api(i 0,1, ,p;ap0 1)的滤波器得到的。只要xN(n)的第一个数据x (0)进入滤波器，滤波器便输出第一个误差信号取样值

ep (N p 1)。这意味着，已知数据x(n)(0 x(n) N 1)是通过对无穷长数据序

包括levense递推关系式和burg算法

列x(n)( n )加窗得到的。将e

p

(n) apix(n i)代入式(3.1.1)得

i 0

p

构成的N阶取样自相关矩阵。因此，用时间平均最小化准则同样可以导出 Yule-Walker方程组，不过方程组中的R要用反取代。取样自相关矩阵左是正定的，因而能够保证所得到的预测误差滤波器是最小相位的，因而也能保证反射系数的模值都小于1，这是使滤波器稳定的充分条件。

用Yule-Walker算法估计AR模型参数的具体步骤为 (1)根据式(3. 1. 3)计算自相关函数。

(2)用Levinson-Durbin算法求解Yule-Walker方程组

(3)用AR参数的估计值，可以计算功率谱

3.1.2 Burg算法

为了克服L-D算法中因估计相关函数给功率谱带来的影响，Burg (1968)提出 一种新的算法，其基本思想是直接从观测的数据利用线性预测器得前向和后向预测的总均方误差之和为最小的准则来估计预测系数，进而通过L-D算法的递推公式求出AR模型优化的参数。

前、后向预测误差的定义式分别为

将式((2. 3. 12b) am(k)分别代入以上两式，得

包括levense递推关系式和burg算法

式中

m am(m)

，表示m阶预测器结构模型用格型实现时的反射系数

m

。

等于直接实现第m次预测系数。

按前后向预测误差均方(平方)误差的总和，有

将式((3. 1. 6)代入式((3. 1. 7)并对反射系数求偏导数，且令得总平均

方误差为最小时得反射系数等于

。求

上式分母中的两项分别是前向与后向最小方误差，可以证明总的最小平方误差为

按L-D递推公式可求得

综上所述，Burg算法步骤为

（1）根据式(2. 3. ?)和前、后向预测误差的定义知道

这时均方差等于预测器平均输出功率，故有

（2）按式(2. 3. 12b)及(3. 1. 8)和(3. 1. 9)递推关系，估计P阶预测系数((AR模型参数)及最小预测误差的均方值，即

包括levense递推关系式和burg算法

（3）按式(3. 1. 6)递推高一阶前、后向预测误差，即

（4）求得的AR模型参数估值，得到功率谱估计

Burg算法的优点在于①频率分辨率高;②所得的AR模型稳定;③计算效率高。然而，Burg方法也有公认的不足。首先，在高信噪比时，它的谱线呈现出分裂;其次对于高阶模型来讲，该方法也可能引入假谱峰;再次，对于噪声中的正弦信号，Burg方法对正弦信号的初始相位呈现出敏感性。特别是对短的数据记录更是如此。

第4章 AR模型功率谱估计流程图

4.1 Levinson-Durbin快速递推法做功率谱估计流程图

包括levense递推关系式和burg算法

4.2 Burg算法做功率谱估计流程图

第5章 通过MATLAB编程做谱估计

待估计的信号为yn=cos(2*pi*f1*n)+cos(2*pi*f2*n)+wn;其中f1,f2为信号两个不同的频率，wn为高斯白噪声，数据采样点数为N

5.1.1 f1=0.1,f2=0.4，N=512时的古典谱估计结果，程序见附录（1.1）

信号加噪声

5

-50

200

400

600

0.2

0.4

0.6

0.8

1

Bartlett法 (l=2）

用相关法作功率谱估计

20100-100

0.2

0.4

0.6

0.8

1

0.2

Bartlett法 (l=4）

0.40.60.81

用周期图法做谱估计

500

Welch(汉宁窗，L=7,64点交叠

包括levense递推关系式和burg算法

5.1.2 f1=0.1,f2=0.4,N=512的AR模型谱估计结果，Levinson-Durbin快速递推法（阶数为4）程序见附录（1.2）和Burg算法（阶数为4）程序见附录（1.3）

43210-1-2-3-4

levinson-Durbin 算法

Burg算法

5.2.1 f1=0.1,f2=0.11，N=512时的古典谱估计结果，程序见附录（1.1修改f2值即可） 信号加噪声Bartlett法 (l=2）

020040060000.20.40.60.81

用相关法作功率谱估计Bartlett法 (l=4）

00.20.40.60.8100.20.40.60.81

用周期图法做谱估计Welch(汉宁窗，L=7,64点交叠

包括levense递推关系式和burg算法

5.2.2 Levinson-Durbin快速递推法(阶数为88)和Burg算法（阶数为70），N=512，程序见附录（1.2）和（1.3），修改其中f2及阶数值即可。

Burg算法

5.3.1 Levinson-Durbin快速递推法(阶数为4)和Burg算法（阶数为4），

f1=0.1,f2=0.4,N=16，程序见附录（1.2）和（1.3），修改其中f2及阶数值即可,N值即可。

levinson-Durbin 算法

Burg算法

第6章 比较分析各种功率谱估计

6.1 古典法与AR模型法比较

6.1.1 在相同的条件下，f1=0.1,f2=0.4，N=512，观察5.1.1中图和5.1.2中图， 我们可以发现在频率间隔比较大是古典法和AR模型法都能估计出信号的频率，但估计的质量是有差别的，总结如下：

（1） 直接用相关函数法和周期图法谱估计质量很差，尤其是周期图法。

（2） 改进后的古典法Bartlett功率谱质量较好，Bartlett分段越多，谱估计好。

包括levense递推关系式和burg算法

（3） 改进后的古典谱估计中Welch法谱估计质量最好。

（4） 基于AR模型的功率谱估计质量明显优于古典谱估计法。

6.1.2在相同的条件下，f1=0.1,f2=0.11，N=512，观察5.2.1和5.2.2中图，我们可以发现但频率间隔变小时，用古典谱估计法已经估计不出两个频率点；用AR模型谱估计法，当阶次一定大时，仍能够分辨出两频率点。由此我们可以得出这样的结论：古典谱估计分辨率低，当频率间隔小到一定程度时，就分辨不出两接近的频率点来；AR模型频率分辨率较高，只要选择合适的阶次，就能分辨出两接近的频率点。 6.1.3 分析

古典谱估计存在着加窗截断的现象，古典谱估计不是真是谱的一致估计，即使改进的古典谱估计法也依然存在着加窗问题，这导致分辨率低和方差性能不好。基于AR模型的谱估计避免了古典法中的加窗截断，所以做出的谱估计很接近真实值。（调试程序中发现AR模型阶次越低，做出的谱估计越平滑，阶次越高，越不平滑，毛刺越多）

6.2 Levinson-Durbin快速递推法和Burg算法数据样点N不同时比较 6.2.1 N=16与N=512比较

在N=16条件下调试程序过程中用Levinson-Durbin快速递推法得出的图形不稳定，每次调试结果都不大相同，用Burg算法数得出的图形很稳定，变化很小。N=512是用两种方法得出的图形变化都不大。 6.2.2 分析

（1）Levinson-Durbin快速递推法需要先估计自相关函序列。 （2）Levinson-Durbin快速递推法数据之外的值假设为0. （3）Burg算法不需要估计自相关序列。

（4）Burg算法不对已知数据之外的数据作假设。

（5）Burg算法总是保证预测误差滤波器是最小相位的。

因此，Burg算法求得的AR模型最稳定；由于Burg算法不需要自相关函数，所以性能优于自相关法。尤其在短数据时，Burg算法明显优越，具有较高的普分辨率。

第7章 总结

谱估计分为古典谱估计和现代谱估计，现代谱估计方法主要是以随机过程的参数模型为基础的，因此，也可以将其称为参数模型方法或模型方法。在短据条件下，参数法谱估计比非参数法具有更高的分辨率。这些算法是通过不同的方法进行功率谱估计。它们不是直接用数据估计功率谱，而是通过对数据建模，然后再估计出系统模型参数。最常用的系统模型是全极点模型，这种模型称为自回归模型(AR模型)，因而，这些谱估计方法有时又称为AR法谱估计。AR法谱估计是用功率谱峰值来表示数据的频率。由于估计AR模型参数通常只需解一组线性方程，因此AR模型得到了深入的研究和广泛应用。

不同的谱估计方法估计出的AR模型的参数又有一些细微的差别，估计 出的功率谱不完全相同，下面用一表格对几种AR模型谱估计算法做一简要总结:

包括levense递推关系式和burg算法

表7. 1

参考文献

[1] 张玲华，郑宝玉. 随机信号处理. 北京：清华大学出版社，2003.9 [2] 钟麟，王峰. MATLAB仿真技术与应用教程. 北京：国防工业出版社，2004.1 [3] 刘志刚 . 基于现代谱佑计理论的信噪分离方法及其应用研究. 成都理工大学硕士学位论文,2006.5

[4] 高西全，丁玉美. 数字信号处理（第三版）.西安：西安电子科技大学出版社，2008.8

附录：

（1.1）古典法 clear all

%产生两个正弦信号与高斯白噪声的叠加 N=512; %定义总点数N f1=0.1;f2=0.4; n=0:(N-1);

wn=randn(1,N); %产生高斯白噪声 %y(n)=cos(2*pi*f1*n)+cos(2*pi*f2*n);

yn=cos(2*pi*f1*n)+cos(2*pi*f2*n)+wn; %信号加噪声 subplot(321);

plot(n,yn); %画图 title('信号加噪声');

%用相关法作功率谱估计

xk=fft(yn,2*N-1); %做2N-1点fft

包括levense递推关系式和burg算法

xk1=abs(xk).^2/N;

rm=ifft(xk1,2*N-1); %用循环卷积发求出自相关序列 M=512; for i=1:M

R(i)=rm(i); end

xk2=fft(R,M); %对相关序列做ifft subplot(323) k=1:M; f=(k-1)./M;

plot(f,10*log10(abs(xk2)));

title('用相关法作功率谱估计'); %用周期图法做谱估计

xk3=fft(yn,N); %做N点fft

xk4=abs(xk3).^2/N; %模平方，取平均 k=1:N; f=(k-1)./N; subplot(325)

plot(f,10*log10(abs(xk4))); title('用周期图法做谱估计'); %Bartlett法 (l=2）

for i=1:N/2 %把信号分成2段

yn1(1,i)=cos(2*pi*f1*i)+cos(2*pi*f2*i);

yn2(1,i)=cos(2*pi*f1*(i+N/2))+cos(2*pi*f2*(i+N/2)); if i==N/2 break; end end

yn1=yn1+randn(1,N/2); %第1段加噪声 yn2=yn2+randn(1,N/2); %第2段加噪声 xk5=fft(yn1,N/2); %做fft xk6=fft(yn2,N/2); %做fft

xk7=(abs(xk5).^2+abs(xk6).^2)./(N); %模平方取均值 k=1:N/2;

f=(k-1)./(N/2); subplot(322)

plot(f,10*log10(abs(xk7))); title('Bartlett法 (l=2）'); %Bartlett法 (l=4）

for i=1:N/4 %把信号分成4段

yn3(1,i)=cos(2*pi*f1*i)+cos(2*pi*f2*i)

包括levense递推关系式和burg算法

yn4(1,i)=cos(2*pi*f1*(i+128))+cos(2*pi*f2*(i+128)); yn5(1,i)=cos(2*pi*f1*(i+256))+cos(2*pi*f2*(i+256)); yn6(1,i)=cos(2*pi*f1*(i+384))+cos(2*pi*f2*(i+384)); if i==N/4; break; end end

yn3=yn3+randn(1,N/4); %第1段加噪声 yn4=yn4+randn(1,N/4); %第2段加噪声 yn5=yn5+randn(1,N/4); %第3段加噪声 yn6=yn6+randn(1,N/4); %第4段加噪声 xk3=fft(yn3,N/4); xk4=fft(yn4,N/4); xk5=fft(yn5,N/4);

xk6=fft(yn6,N/4); %对各段做fft

xk7=(abs(xk3).^2+abs(xk4).^2+abs(xk5).^2+abs(xk6).^2)./(N); %模平方取均值 k=1:N/4;

f=(k-1)./(N/4); subplot(324)

plot(f,10*log10(abs(xk7))); title('Bartlett法 (l=4）');

%Welch(汉宁窗，L=7,64点交叠） U=0; for i=1:7

for j=1:128

wh(i,j)=0.5*(1-cos(2*pi*j/(127))); %汉宁窗

U= U+wh(i,j).^2/7; %窗函数平均功率

yn7(i,j)=(cos(2*pi*f1*(j+64*(i-1)))+cos(2*pi*f2*(j+64*(i-1)))).*wh(i,j); %信号分为重叠64点的7段 end end for i=1:7

yn7(i,:) = yn7(i,:)+randn(1,128); %对每段加噪声 xk_7(i,:)=fft( yn7(i,:),N/4); %对每段做fft end

xk8=zeros(1,128); for i=1:7

xk8=(xk8+ abs(xk_7(i,:)).^2./U); %求功率谱估计 end

xk8=xk8/7; k=1:N/4;

f=(k-1)./(N/4); subplot(326)

包括levense递推关系式和burg算法

plot(f,10*log10(abs(xk8)));

title('Welch(汉宁窗，L=7,64点交叠');

（1.2） levinson-Durbin 算法 clear all %for i=1:25

%产生两个正弦信号与高斯白噪声的叠加 N=512; %定义总点数N p=4; %阶数

f1=0.1;f2=0.4; %定义两个频率 n=0:(N-1);

wn=randn(1,N); %产生高斯白噪声 %y(n)=cos(2*pi*f1*n)+cos(2*pi*f2*n);

yn=cos(2*pi*f1*n)+cos(2*pi*f2*n)+wn; %信号加噪声 %subplot(311); %plot(n,yn);

%title('信号加噪声');

%用相关法作功率谱估计 xk=fft(yn,2*N-1); xk1=abs(xk).^2/N;

rm=ifft(xk1,2*N-1); %用循环卷积法求出自相关序列 %levinson-Durbin 算法

aa(1,1)=-rm(2)/rm(1) %求初值

p(1)=rm(1)*(1-abs(aa(1,1)^2)) %求初值 sum=0;m=2; %初始化 for m=2:(4) %循环的阶数

for i=1:m-1

sum=sum+aa(m-1,i)*rm(m+1-i); end

aa(m,m)=-(rm(m+1)+sum)/p(m-1); %第m阶的反射系数

p(m)=p(m-1)*(1-aa(m,m)^2); %第m阶的预测误差功率 for i=1:m-1

aa(m,i)=aa(m-1,i)+aa(m,m)*aa(m-1,m-i); %Levinson关系式 end

sum=0; %在下次循环前给sum赋初值 end %levinson-Durbin 算法结束 B=(p(m)).^0.5;

aa(m,:) %利用levinson-Durbin 算法求出的参数 A=[1,aa(m,:)];

H=freqz(B,A) %求频响 k=1:512;

f=(k-1)./(1024); %subplot(312);

包括levense递推关系式和burg算法

figure

plot(f,10*log10(abs(H).^2)); %画图 title('levinson-Durbin 算法'); %end

(1.3) Burg 算法 %for i=1:10 clear all

%产生两个正弦信号与高斯白噪声的叠加 N=512; %定义总点数N p=4; % order

f1=0.1;f2=0.4; %定义两个频率 n=0:(N-1);

wn=randn(1,N); %产生高斯白噪声 %y(n)=cos(2*pi*f1*n)+cos(2*pi*f2*n);

yn=cos(2*pi*f1*n)+cos(2*pi*f2*n)+wn; %信号加噪声 %subplot(321); %plot(n,yn);

%title('信号加噪声');

%用相关法作功率谱估计 %xk=fft(yn,2*N-1); %xk1=abs(xk).^2/N;

%rm=ifft(xk1,2*N-1); %循环卷积法求出相关序列 %Burg suan fa

suma=0;sumb=0; %初始化 for n=1:N %赋初值 ef(1,n)=yn(n); eb(1,n)=yn(n); end

for m=2:(p+1) %循环的阶数 for n=m:N

suma=suma+ef(m-1,n)*eb(m-1,n-1);

sumb=sumb+ef(m-1,n)^2+eb(m-1,n-1)^2; end

k(m-1)=-2*suma/sumb; %根据前后向预测均方误差值和最小来求反射系数

suma=0;sumb=0; %每次循环前赋零值 for n=2:N

ef(m,n)=ef(m-1,n)+k(m-1)*eb(m-1,n-1); %前向误差 eb(m,n)=eb(m-1,n-1)+k(m-1)*ef(m-1,n); %后向误差 end end

aa(1,1)=k(1); %赋初值 for m=2:p

包括levense递推关系式和burg算法

for i=1:m-1

%aa(m-1,m-1)=k(m-1);

aa(m,i)=aa(m-1,i)+k(m)*aa(m-1,m-i); %利用levinson-Durbin 算法求出的参数

aa(m,m)=k(m); end

end aa(m,:); %求出的参数 A=[1,aa(m,:)]; B=1;

H=freqz(B,A); %求频响 k=1:512;

f=(k-1)./(1024); figure

plot(f,10*log10(abs(H).^2)); title('Burg算法'); %end

% 画图

本文来源：https://www.bwwdw.com/article/fg34.html