第8章参数化建模

《第8章参数化建模》由会员分享，可在线阅读，更多相关《第8章参数化建模（8页珍藏版）》请在装配图网上搜索。

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

2、内的信息，如脉冲响应或系统输入和输出时间序列，求数字滤波器有理传递函数分子和分母多项式系数。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 regre

3、ssion, AR)线性模型，它可用一个全极点IIR滤波器脉冲响应来逼近，全极点滤波器的传递函数具有如下形式：1H(z)1n(8-2)1 +a(2)z +a( n+1)zMATLAB言号处理工具箱函数Ipc利用自回归(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) 调用函数levi n

4、son 实现Levin so n-Durbi n 递推算法，求出系数 a。函数levi nson 的调用格 式为：a=levi nson(r,n)其中，r为x的自相关序列；n为AR模型阶次。【例8-1】设信号x是一个带白噪声的 4阶IIR滤波器的脉冲响应，IIR滤波器的传递函数 无零点，其分母多项式系数为a=1 0.1 0.1 0.1 0.1，用线性预报法建立其AR模型。%Samp8_1randn( state,0);%设置随机函数的状态a=1 0.1 0.1 0.1 0.1;%滤波器分母多项式系数x=impz(1,a,10)+ra ndn (10,1)/20;%求得带噪声的滤波器脉冲响应%采

5、用第一种方法r(1:length(x)-1)=;%该序列的前面部分(相当于x的长度-1)受到边界的影响，扣除 aa=levi nson (r,4)%用Levis on-Durbin递推算法采用自相关序列求出系数 aa%采用第二种方法a1=lpc(x,4)% 直接调用 lpc ( Lin ear Predictor Coefficie nts )函数求得结果该程序的运行结果为：aa =1.00000.18490.12790.11140.1839a1 =1.00000.18490.12790.11140.1839程序中两种方法求得的结果完全相同。若信号不包含随机噪声，用上面的程序求得信号AR模型和

6、所采用全极点滤波器模型完全相同。即使加入随机噪声，得到结果也与原模型参 数a接近。这说明，如果一个信号x是IIR滤波器的脉冲响应，可以利用函数lpc建立的向量a作为滤波器参数化模型。8. 1. 2 Prony 法(ARM/模型)如果已知某滤波器(或系统)的脉冲响应序列，可以采用Prony方法求得该滤波器(IIR) 传递函数的系数。MATLAB言号处理工具箱提供了Prony法建模的函数prony,其调用格式为：b,a=pr ony(h,nb,na)式中，h为时间域脉冲响应序列；nb,na分别为滤波器分子和分母多项式阶数；b,a分别为滤波器分子和分母多项式系数向量。这就是说，要想建立起ARMA莫型

7、，必须知道传递函数分子和分母多项式的最高阶次。滤波器传递函数表示为：(8-3)、 B(z) b(1)+b(2)z+b(nb + 1)zbHz1naA(z) a(1) + a (2)z 十+a( na + 1)zProny法是用一个IIR滤波器脉冲响应逼近信号h。【例8-2】建立一个4阶Butterworth 数字滤波器，归一化截止频率为 0.2，求其脉冲响应， 再用脉冲响应系数和prony函数求滤波器模型，并与原始模型进行比较。%Samp8_2b,a=butter(4,0.2) % 设计4阶截止频率为 0.2(归一化)的Butterworth 滤波器h=impz(b,a,26); % 求26个

8、点组成的滤波器的脉冲响应bb,aa=prony(h,4,4) %运用脉冲响应h建立4阶系统的分子和分母多项式系数运行结果为：b = 0.0048 0.0193 0.0289 0.0193 0.0048a = 1.0000 -2.3695 2.3140 -1.0547 0.1874bb = 0.0048 0.0193 0.0289 0.0193 0.0048aa = 1.0000 -2.3695 2.3140 -1.0547 0.1874程序运行结果表明，用函数prony建模和原始滤波器模型完全相同，显然函数prony所建立的新滤波器的脉冲响应和信号h是完全吻合的。因此，若已知一个系统的脉冲响应

9、和该系统分子分母多项式阶数，可以运用该方法求得该系统的传递函数。8. 1. 3 Steiglity-McBride法(ARMA模型)已知某滤波器(或系统)的脉冲响应序列建立滤波器传递函数的另一种方法是Steiglity-McBride 法。这种方法利用 Steiglity-McBride迭代建模，使系统 b(z)/a(z) 的脉冲响应x和输出信号x的均方差最小，即(8-4)MATLAB言号处理工具箱函数 stmcb实现这种方法建模， 可用于滤波器设计和系统辨识。 其调用格式为：b,a=stmcb(x ,nb,n a, niter,ai)式中,x为系统的脉冲响应；nb,na分别为系统传递函数b(

10、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 滤波器

11、，截止频率为 0.2 ，求其脉冲响应，并用脉冲 响应数据和函数 stmcb 建立新滤波器模型，比较两个滤波器的频率特性。%Samp8_3b,a=butter（6,0.2） % 产生 6 阶截止频率为 0.2 的 Butterworth 滤波器 h=filter（b,a,1 zeros（1,100）; % 采用输入为脉冲函数的方法求得系统的脉冲响应 bb,aa=stmcb（h,6,6） % 采用 stmcb 方法求得系统传递函数多项式系数 运行结果为：b = 0.0003 0.0020 0.0051 0.0068 0.0051 0.0020 0.0003 a = 1.0000 -3.5794 5

12、.6587 -4.9654 2.5295 -0.7053 0.0838 bb = 0.0003 0.0020 0.0051 0.0068 0.0051 0.0020 0.0003 aa = 1.0000 -3.5794 5.6587 -4.9654 2.5295 -0.7053 0.0838结果表明，函数 stmcb 自滤波器脉冲响应建立的模型和原滤波器模型完全相同。 为了比较 lpc ， prony 和 stmcb 这三个函数 , 下例给出了三种方法建模的误差 【例 8-4 】采用例 8-1 中的滤波器，用阶数 3 模拟比较函数 lpc 、 prony 和 stmcb 三种建模 方法的误差。

13、%Samp8_4randn(seed,0); % 设置随机数种子x=impz(1,1 0.1 0.1 0.1,10)+randn(10,1)/20; %lpc 函数建立传递函数系数 运用 prony 函数建立传递函数系数 运用 stmcb 函数建立传递函数系数与 lpc 函数确定模型脉冲响应的每个数的误差与 lpc 函数确定模型脉冲响应的总误差a1=lpc(x,3); % 运用b2,a2=prony(x,3,3); % b3,a3=stmcb(x,3,3); % %计算误差disp(函数LPC输出:) err1=x-impz(1,a1,10); %xsumerr1=sum(err1.A2) %

14、x disp(函数 PRONY输出:)err2=x-impz(b2,a2,10); %x sumerr2=sum(err2.A2) %x与 prony 函数确定模型脉冲响应的每个数误差 与 prony 函数确定模型脉冲响应的总误差disp(函数 STMCB俞出:) err3=x-impz(b3,a3,10); %x与 stmcb 函数确定模型脉冲响应的每个数误差sumerr3=sum（err3.A2） %x 与 stmcbc 函数确定模型脉冲响应的总误差 程序输出结果为：函数 LPC输出：sumerr1=0.0219; 函数 PRON输出：sumerr2 =0.0147;函数 STMC输出：s

15、umerr3=0.0039该例表明：在相同的阶数条件下，函数 stmcb 精度最好， prony 次之， lpc 较差。8.2 频率域建模前面一节我们根据系统（滤波器）的时间域信号（脉冲响应或输入、输出信号）建立系 统传递函数。但如果已知某系统（滤波器）的频率响应，能否建立起该系统（滤波器）的模 型呢？这就涉及到频率域建模。频率域建模就是由滤波器（系统）的频率响应建立其模型。 根据频率域类型分为面向模拟滤波器的 s 域内建模和面向数字滤波器的 z 域内建模。8.2.1模拟滤波器的s域建模模拟滤波器传递函数的一般形式为:H(s)B(s)AS)b(1)snb b(2)snbb(nb 1)a(1)s

16、na a(2)sna J - a(na 1)(8-5)MATLAB!号处理工具箱函数invfreqs用于由给定的频率响应数据求形如上式的模拟滤 波器传递函数。调用格式为：b,a=in vfreqs(h,w ,nb,n a,wt n iter)式中，h为复频率响应向量；w为对应的频率向量；nb,na分别为模拟滤波器分子和分母多 项式阶次；wt为与w同维数的权系数向量，调整对频率的拟合误差；niter为迭代次数。若已知复频率响应的幅值向量 mag和相位向量phase，在 MATLAB中要采用 h=mag.*exp(j*phase) 将其复合为复数形式。【例8-5】已知滤波器原始传递函数中b=1 2

17、 3 2,a=1 2 3 2 1,试求其频率响应，并利用频率响应数据求原始滤波器的传递函数。%Samp8_5disp(滤波器初始模型:);b=1 2 3 2a=1 2 3 2 1h,w=freqs(b,a,64); % 计算64个点的模拟滤波器频率响应disp(INVFREQS函数得到的模型：)bb,aa=invfreqs(h,w,3,4) % 运用 invfreqs 函数求滤波器系数程序的运行结果为：滤波器初始模型：b=1 2 3 2a=1 2 3 2 1INVFREQS 函数得到的模型:bb =1.0000 2.0000 3.0000 2.0000aa=1.0000 2.0000 3.00

18、00 2.0000 1.0000运行结果表明，由函数in vfreqs 从频率响应数据辨识的滤波器模型和原始模型完全相 同。8.2.2数字滤波器的z域内建模数字滤波器传递函数的一般形式为H四A(z)MATLAB!号处理工具箱函数 波器的传递函数。调用格式为：b(1)+b(2)z+b(nb + 1)zb( 86)-a(1)a(2)za(na 1)zain vfreqz用于由给定的频率响应数据求形如上式的数字滤b,a=i nvfreqz(h,w, nb, na,wt, ni ter)式中，h为复频率响应向量；w为对应的频率向量；nb,na分别为数字滤波器分子和分母多 项式阶次；wt为与w同维数的权

19、系数向量，调整对频率的拟合误差；niter为迭代次数。函数invfreqs 和invfreqz利用等误差方法由给定的频率响应数据来辨识最佳模型，当niter缺省时，函数利用下面的算法：n2min 送 wt(k) h(k)A(w(k) B(w(k)( 8-7)b,a k式中，A w(k)和B w(k)分别表示多项式a,b在频率w(k)处的Fourier变换，wt为权系 数向量，n为频率点数(h和w的长度)。这种算法不能保证辨识模型的稳定性。当迭代次数 niter确定后，函数用输出误差迭代算法进行运算：minb,aE wt(k) h(k)-k 二1A(w(k)B(w(k)(8-8)这样可保证辨识模

20、型的稳定性，且提高辨识模型的精度。【例8-6】已知原始滤波器为4阶Butterworth低通数字滤波器，归一化截止频率为0.4。用其频率响应数据产生一个三阶滤波器。试比较采用等误差和输出误差迭代算法的结果。%Samp8_6b,a=butter(4,0.4); %产生 4 阶 Butterworth滤波器h,w=freqz(b,a,64); %计算频率响应bb,aa=invfreqz(h,w,3,3); %采用 invfreqz求解滤波器的系数wt=o nes(size(w); % 产生权向量niter=30; %迭代次数bbb,aaa=i nvfreqz(h,w,3,3,wt, niter);

21、 %采用权向量和迭代次数求滤波器系数h1,w1=freqz(bb,aa,w); %计算模型1的频率响应h2,w2=freqz(bbb,aaa,w); %计算模型2的频率响应plot(w/pi,abs(h),k,w1/pi,abs(h1),k:,w2/pi,abs(h2),b.); %绘制频率响应xlabel(归一化频率);ylabel( 振幅)legend(理想模型,未采用权向量和迭代次数,采用权向量和迭代次数);grid on sumerr1=sum(abs(h-h1).A2) %输出模型 1 误差sumerr2=sum(abs(h-h2).A2) % 输出模型 2 误差图8-1例8-6的两

22、种算法估计模型的幅频特性的比较程序命令窗口的输出为： sumerr1=0.0200 ； sumerr2 = 0.0096 。程序运行的图形显示为 图 8-1 。程序运行结果表明，输出误差迭代法（预估算法）具有较好的拟合精度。第八章 习题1. 用MATLAB!数产生一个5阶低通Chebyshev I型数字滤波器（归一化截止频率为0.4 ）的脉冲响应序列，再用 Prony 法建模并与原滤波器模型进行比较。2. 用MATLAB!数产生一个5阶低通Butterworth 数字滤波器（归一化截止频率为 0.4） 并计算其脉冲响应序列，再用 Steiglitz-McBride 法建模并与原滤波器模型进行比 较。3. 用MATLABi数产生一个低通 Chebyshev I型高通数字滤波器，设计指标为：通带截 止频率Fc=2.5kHz,通带波纹为1dB,阻带衰减大于 20dB。试求（1）高通数字滤波 器的系统函数； （2）滤波器的频率响应； （3）由频率响应数据辨识滤波器模型并与 原滤波器模型进行比较。