第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阶自回归(automaticregression,AR)线性模型，它可用一个全极点IIR滤波器脉冲响应来逼近，全极点滤波器的传递函数具有如下形式：nznazazH)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=[10.10.10.10.1]，用线性预报法建立其AR模型。%Samp8_1randn('state',0);%设置随机函数的状态a=[10.10.10.10.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（LinearPredictorCoefficients）函数求得结果该程序的运行结果为：aa=1.00000.18490.12790.11140.1839a1=1.00000.18490.12790.11140.1839程序中两种方法求得的结果完全相同。若信号不包含随机噪声，用上面的程序求得信号AR模型和所采用全极点滤波器模型完全相同。即使加入随机噪声，得到结果也与原模型参数a接近。这说明，如果一个信号x是IIR滤波器的脉冲响应，可以利用函数lpc建立的向量a作为滤波器参数化模型。8.1.2Prony法（ARMA模型）如果已知某滤波器（或系统）的脉冲响应序列，可以采用Prony方法求得该滤波器（IIR）传递函数的系数。MATLAB信号处理工具箱提供了Prony法建模的函数prony,其调用格式为：[b,a]=prony(h,nb,na)式中，h为时间域脉冲响应序列；nb,na分别为滤波器分子和分母多项式阶数；b,a分别为滤波器分子和分母多项式系数向量。这就是说，要想建立起ARMA模型，必须知道传递函数分子和分母多项式的最高阶次。滤波器传递函数表示为：nanbznaazaaznbbzbbzAzBzH)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.00480.01930.02890.01930.0048a=1.0000-2.36952.3140-1.05470.1874bb=0.00480.01930.02890.01930.0048aa=1.0000-2.36952.3140-1.05470.1874程序运行结果表明，用函数prony建模和原始滤波器模型完全相同，显然函数prony所建立的新滤波器的脉冲响应和信号h是完全吻合的。因此，若已知一个系统的脉冲响应和该系统分子分母多项式阶数，可以运用该方法求得该系统的传递函数。8.1.3Steiglity-McBride法(ARMA模型)已知某滤波器（或系统）的脉冲响应序列建立滤波器传递函数的另一种方法是Steiglity-McBride法。这种方法利用Steiglity-McBride迭代建模，使系统b(z)/a(z)的脉冲响应x’和输出信号x的均方差最小，即02,)()(minibaixix（8-4）MATLAB信号处理工具箱函数stmcb实现这种方法建模，可用于滤波器设计和系统辨识。其调用格式为：[b,a]=stmcb(x,nb,na[,niter[,ai]])式中,x为系统的脉冲响应；nb,na分别为系统传递函数b(z)/a(z)分子和分母多项式的阶次；niter为迭代计算次数；缺省值为5；ai为分母系数的初步估计，可缺省；b,a分别为系统（IIR型滤波器）分子和分母多项式系数向量。系统(IIR滤波器)传递函数具有与（8-3）式相同的形式。函数的另一种调用形式适用于系统的输入和输出信号作为函数输入：[b,a]=stmcb(x,u,nb,na[,niter[,ai]])式中，x和u分别为系统（滤波器）输出和输入信号，其余参数同上。函数的这种用法可以根据滤波器的输入和输出（而不是前面的脉冲响应）求得滤波器传递函数。【例8-3】建立一个6阶Butterworth滤波器，截止频率为0.2，求其脉冲响应，并用脉冲响应数据和函数stmcb建立新滤波器模型，比较两个滤波器的频率特性。%Samp8_3[b,a]=butter(6,0.2)%产生6阶截止频率为0.2的Butterworth滤波器h=filter(b,a,[1zeros(1,100)]);%采用输入为脉冲函数的方法求得系统的脉冲响应[bb,aa]=stmcb(h,6,6)%采用stmcb方法求得系统传递函数多项式系数运行结果为：b=0.00030.00200.00510.00680.00510.00200.0003a=1.0000-3.57945.6587-4.96542.5295-0.70530.0838bb=0.00030.00200.00510.00680.00510.00200.0003aa=1.0000-3.57945.6587-4.96542.5295-0.70530.0838结果表明，函数stmcb自滤波器脉冲响应建立的模型和原滤波器模型完全相同。为了比较lpc，prony和stmcb这三个函数,下例给出了三种方法建模的误差【例8-4】采用例8-1中的滤波器，用阶数3模拟比较函数lpc、prony和stmcb三种建模方法的误差。%Samp8_4randn('seed',0);%设置随机数种子x=impz(1,[10.10.10.1],10)+randn(10,1)/20;%a1=lpc(x,3);%运用lpc函数建立传递函数系数[b2,a2]=prony(x,3,3);%运用prony函数建立传递函数系数[b3,a3]=stmcb(x,3,3);%运用stmcb函数建立传递函数系数%计算误差disp('函数LPC输出:')err1=x-impz(1,a1,10);%x与lpc函数确定模型脉冲响应的每个数的误差sumerr1=sum(err1.^2)%x与lpc函数确定模型脉冲响应的总误差disp('函数PRONY输出:')err2=x-impz(b2,a2,10);%x与prony函数确定模型脉冲响应的每个数误差sumerr2=sum(err2.^2)%x与prony函数确定模型脉冲响应的总误差disp('函数STMCB输出:')err3=x-impz(b3,a3,10);%x与stmcb函数确定模型脉冲响应的每个数误差sumerr3=sum(err3.^2)%x与stmcbc函数确定模型脉冲响应的总误差程序输出结果为：函数LPC输出:sumerr1=0.0219;函数PRONY输出:sumerr2=0.0147;函数STMCB输出:sumerr3=0.0039该例表明：在相同的阶数条件下，函数stmcb精度最好，prony次之，lpc较差。8.2频率域建模前面一节我们根据系统（滤波器）的时间域信号（脉冲响应或输入、输出信号）建立系统传递函数。但如果已知某系统（滤波器）的频率响应，能否建立起该系统（滤波器）的模型呢？这就涉及到频率域建模。频率域建模就是由滤波器（系统）的频率响应建立其模型。根据频率域类型分为面向模拟滤波器的s域内建模和面向数字滤波器的z域内建模。8.2.1模拟滤波器的s域建模模拟滤波器传递函数的一般形式为：)1()2()1()1()2()1()()()(11naasasanbbsbsbsAsBsHnananbnb（8-5）MATLAB信号处理工具箱函数invfreqs用于由给定的频率响应数据求形如上式的模拟滤波器传递函数。调用格式为：[b,a]=invfreqs(h,w,nb,na[,wt[niter]])式中，h为复频率响应向量；w为对应的频率向量；nb,na分别为模拟滤波器分子和分母多项式阶次；wt为与w同维数的权系数向量，调整对频率的拟合误差；niter为迭代次数。若已知复频率响应的幅值向量mag和相位向量phase，在MATLAB中要采用h=mag.*exp(j*phase)将其复合为复数形式。【例8-5】已知滤波器原始传递函数中b=[1232],a=[12321],试求其频率响应，并利用频率响应数据求原始滤波器的传递函数。%Samp8_5disp('滤波器初始模型:');b=[1232]a=[12321][h,w]=freqs(b,a,64);%计算64个点的模拟滤波器频率响应disp('INVFREQS函数得到的模型:')[bb,aa]=invfreqs(h,w,3,4)%运用invfreqs函数求滤波器系数程序的运行结果为：滤波器初始模型:b=1232a=12321INVFREQS函数得到的模型:bb=1.00002.00003.00002.0000aa=1.00002.00003.00002.00001.0000运行结果表明，由函数invfreqs从频率响应数据辨识的滤波器模型和原始模型完全相同。8.2.2数字滤波器的z域内建模数字滤波器传递函数的一般形式为: