第8章 参数化建模

第8章 参数化建模

由前几章的讨论可见，无论是模拟滤波器还是数字滤波器，均是根据输入信号和系统的传递函数求得系统的输出信号。反过来，能否通过系统的输入信号和输出信号，或滤波器的频率响应或脉冲响应求得系统的传递函数呢？这在一定的假设下是可以得到的。这种技术涉及本章要讲的参数化建模。所谓参数化建模就是根据未知系统的某些信息（如脉冲响应、频率响应或输入输出序列）建立该系统的有理传递函数模型。这项技术用于求一个信号、系统或过程的模型参数，并广泛应用于语言分析、数据压缩、高分辨谱估计、信号处理等领域。参数化建模分为时间域建模和频率域建模两类，本章将分别介绍。

8.1 时间域建模

时间域建模就是给定系统时间域内的信息，如脉冲响应或系统输入和输出时间序列，求数字滤波器有理传递函数分子和分母多项式系数。MATLAB 中提供了三种时间域建模函数，下面分别叙述。

8.1.1 线性预报法(AR 模型)

如果一个信号的各个采样值不是独立的，其任一时刻k 的采样值x(k)是它过去的n 个采样值及一个白噪声序列k 时刻的值u(k)线性组合而成（线性预报）,即

x(k)=-a(2)x(k-1)-a(3)x(k-2)…-a(n)x(k-n-1)-a(n+1)x(k-n)+u(k)

（8-1）

如果知道参数个数n ，根据这个信号序列可以得到这n 个系数的值。数字信号处理中称这个信号x 符合n 阶自回归(automatic regression, AR)线性模型，它可用一个全极点IIR 滤波器脉冲响应来逼近，全极点滤波器的传递函数具有如下形式： n z

n a z a z H --++++=)1()2(11)(1 （8-2） MATLAB 信号处理工具箱函数lpc 利用自回归(AR)模型的相关法来求出AR 模型系数a ，也是全极点滤波器模型。其调用格式为：

[a,g]=lpc(x,n)

其中，x 为信号，是一个实时间序列；n 为AR 模型阶次；a 为AR 模型系数，a=[1,a(2),…,a(n+1);g 为AR 模型的增益。这就是说，要建立起AR 模型，必须知道其阶次。 函数lpc 的算法可简要地分为几个步骤：（1）调用函数xcorr 求信号x 的相关估计。（2）调用函数levinson 实现Levinson-Durbin 递推算法，求出系数a 。函数levinson 的调用格式为：

a=levinson(r,n)

其中，r 为x 的自相关序列；n 为AR 模型阶次。

【例8-1】设信号x 是一个带白噪声的4阶IIR 滤波器的脉冲响应，IIR 滤波器的传递函数无零点，其分母多项式系数为a=[1 0.1 0.1 0.1 0.1]，用线性预报法建立其AR 模型。 %Samp8_1

randn('state',0); %设置随机函数的状态

a=[1 0.1 0.1 0.1 0.1]; %滤波器分母多项式系数

x=impz(1,a,10)+randn(10,1)/20; %求得带噪声的滤波器脉冲响应

%采用第一种方法

r=xcorr(x); %求信号x 的自相关r ，参考第9章的该函数用法。

r(1:length(x)-1)=[];%该序列的前面部分(相当于x 的长度-1)受到边界的影响，扣除 aa=levinson(r,4) %用Levison-Durbin 递推算法采用自相关序列求出系数aa

%采用第二种方法

a1=lpc(x,4) %直接调用lpc （Linear Predictor Coefficients ）函数求得结果

该程序的运行结果为：

aa = 1.0000 0.1849 0.1279 0.1114 0.1839

a1 = 1.0000 0.1849 0.1279 0.1114 0.1839

程序中两种方法求得的结果完全相同。若信号不包含随机噪声，用上面的程序求得信号AR 模型和所采用全极点滤波器模型完全相同。即使加入随机噪声，得到结果也与原模型参数a 接近。这说明，如果一个信号x 是IIR 滤波器的脉冲响应，可以利用函数lpc 建立的向量a 作为滤波器参数化模型。

8.1.2 Prony 法（ARMA 模型）

如果已知某滤波器（或系统）的脉冲响应序列，可以采用Prony 方法求得该滤波器（IIR ）传递函数的系数。MATLAB 信号处理工具箱提供了Prony 法建模的函数prony,其调用格式为：

[b,a]=prony(h,nb,na)

式中，h 为时间域脉冲响应序列；nb,na 分别为滤波器分子和分母多项式阶数；b,a 分别为滤波器分子和分母多项式系数向量。这就是说，要想建立起ARMA 模型，必须知道传递函数分子和分母多项式的最高阶次。

滤波器传递函数表示为：

()na nb

z

na a z a a z nb b z b b z A z B z H ----++++++++==)1()2()1()1()2()1()()(11 （8-3） Prony 法是用一个IIR 滤波器脉冲响应逼近信号h 。

【例 8-2】建立一个4阶Butterworth 数字滤波器，归一化截止频率为0.2，求其脉冲响应，再用脉冲响应系数和prony 函数求滤波器模型,并与原始模型进行比较。

%Samp8_2

[b,a]=butter(4,0.2) %设计4阶截止频率为0.2(归一化)的Butterworth 滤波器 h=impz(b,a,26); %求26个点组成的滤波器的脉冲响应

[bb,aa]=prony(h,4,4) %运用脉冲响应h 建立4阶系统的分子和分母多项式系数 运行结果为：

b = 0.0048 0.0193 0.0289 0.0193 0.0048

a = 1.0000 -2.3695 2.3140 -1.0547 0.1874

bb = 0.0048 0.0193 0.0289 0.0193 0.0048

aa = 1.0000 -2.3695 2.3140 -1.0547 0.1874

程序运行结果表明，用函数prony 建模和原始滤波器模型完全相同，显然函数prony 所建立的新滤波器的脉冲响应和信号h 是完全吻合的。因此，若已知一个系统的脉冲响应和该系统分子分母多项式阶数，可以运用该方法求得该系统的传递函数。

8.1.3 Steiglity-McBride 法(ARMA 模型)

已知某滤波器（或系统）的脉冲响应序列建立滤波器传递函数的另一种方法是Steiglity-McBride 法。这种方法利用Steiglity-McBride 迭代建模，使系统b(z)/a(z)的脉冲响应x ’和输出信号x 的均方差最小，即 ∑∞='-02

,)()(min i b a i x i x （8-4） MATLAB 信号处理工具箱函数stmcb 实现这种方法建模，可用于滤波器设计和系统辨识。其调用格式为：

[b,a]=stmcb(x,nb,na[,niter[,ai]])

式中,x 为系统的脉冲响应；nb,na 分别为系统传递函数b(z)/a(z)分子和分母多项式的阶次；