数字滤波器设计

1、数字滤波器设计通信与电子信息当中，在对信号作分析与处理时，常会用到有用信号叠加无用噪声的问题。这些噪声信号有的是与信号同时产生的，有的是在传输过程中混入的， 在接收的信号中，必须消除或减弱噪声干扰，这是信号处理中十分重要的问题。根 据有用信号与噪声的不同特性，消除或减弱噪声，提取有用信号的过程就称为滤波。滤波器的种类很多，实现方法也多种多样，本章利用matlab 来进行数字滤波器的设计。数字滤波器是一离散时间系统，它对输入序列x(n) 进行加工处理后，输出序列y(n) ，并使 y(n) 的频谱与 x(n) 的频谱相比发生某种变化。由 dsp理论得知，无限长冲激响应(iir) 需要递归模型来实现

2、，有限长冲激响应(fir) 滤波器可以采用递归的方式也可采用非递归的方式实现。本章把fir 与iir 滤波器分别用 matlab 进行分析与设计。数字滤波器的结构参看数字信号处理一书。数字滤波器的设计一般经过三个步骤: 1( 给出所需滤波器的技术指标。2( 设计一个 h(z)，使其逼近所需要的技术指标。3( 实现所设计的 h(z)。4.1 iir数字滤波器设计设计 iir 数字滤波器的任务就是寻求一个因果、物理可实现的系统函数h(z) ，j 使它的频响 h(e)满足所希望得到的低通频域指标，即通带衰减a、阻带衰减 a、 pr 通带截频 、阻带截频 。而其它形式的滤波器由低通的变化得到。pr 采

3、用间接法设计 iir 数字滤波器就是按给定的指标，先设计一个模拟滤波器，进而通过模拟域与数字域的变换，求得物理可实现的数字滤波器。从模拟滤波器变换到数字滤波器常用的有: 脉冲响应不变法和双线性变换法。iir 滤波器的设计过程如下, 数字频域指标模拟频域指标设计模拟滤波器h(s) 设计数字滤波器 h(z) 1. 模拟滤波器简介模拟滤波器的设计方法已经发展得十分成熟，常用的高性能模拟低通滤波器有巴特沃斯型、切比雪夫型和椭圆型，而高通、带通、带阻滤波器则可以通过对低通进行频率变换来求得。必须指出，这三种滤波器都是非线性的相频特性。具体设计方法参见数字信号处理一书。(1) 巴特沃斯低通滤波器的特点是:

4、 通、阻带均为单调下降。这种单调下降的特性使得系统的误差分布不均匀。在设计中，如果在通带满足指标，阻带指标就过于好了。或者说，在阻带满足指标，通带指标也过于好了。这样不利于以最小的阶数来满足设计指标。(2) 切比雪夫模拟低通滤波器分成 型和 型: 型是通带等波纹，阻带单调 下降; 型滤波器是通带单调下降，阻带等波纹。切比雪夫的等波纹特性使得 可以用较小的阶数设计出符合指标的滤波器。(3) 椭圆滤波器特点是通、阻带均为等波纹。椭圆滤波器的误差均匀分布在通、 阻带，比较上面几种滤波器，在同样衰耗指标下，椭圆滤波器所要求的阶数最小。 在同样阶数情况下，椭圆滤波器的通带到阻带的变化最陡峭、性能最好。2

5、( 模拟滤波器映射到数字滤波器从模拟滤波器变换到数字滤波器常用的有: 冲激响应不变法和双线性变换法。(1) 冲激响应不变法冲激响应不变法的基本思路为: 设模拟滤波器的单位冲激响应为h(t) ，数 a 字滤波器的单位脉冲响应h(n) ，是对 h(t) 的均匀取样，即 a 拉氏逆变换抽样 z 变换, h(s)h(t) h(n)= h(nt) h(z) 冲激响应不变法由于保持了模拟滤波器的冲激响应在取样时刻的取样值，所以具有较好的时间特性。但由于从模拟域(s 域)变换到数字域 (z 域) 的映射关系不是一一对应的关系，因而频谱混叠难以避免, 所以在实际应用中有它的局限性。3( matlab滤波器设计

6、函数简介使用 matlab 信号处理工具箱中提供的函数可以很容易地实现iir 滤波器设计。 常用的 iir 数字滤波器设计函数有 : (1)iir滤波器阶数 (order) 选择buttord - 巴特沃斯 (butterworth)滤波器阶数选择cheb1ord - 切比雪夫 (chebyshev)i 型滤波器阶数选择cheb2ord - 切比雪夫 (chebyshev)ii型滤波器阶数选择ellipord - 椭圆(elliptic)滤波器阶数选择buttord用来确定巴特沃斯滤波器的阶数: n,wc = buttord (wp,wr,ap,ar,options) ap,ar 为通带最大衰

7、减和阻带最小衰减，以db为单位 ; wp,wr为归一化通带截频和阻带截频。options= low, high , bandpass, stop , s, 分别对应低通、高通、带通、带阻、模拟滤波器，默认情况下，为低通或带通。对于带通和带阻滤波器 wp=wp1,wp2 ， wp1和 wp2为通带的上、下截频， wr=wr1,wr2，wr1和 wr2分别为阻带的上、下截频。若为模拟滤波器定阶，令options = s, 则 wp和 wr以( 归一或非归一 ) 给出，因此可大于1。输出参数 n为滤波器的阶数， wc为归一化 3db截频。若 ap=3db ，则 wc= wp 。其他滤波器的定阶函数类

8、似buttord的格式 : 切比雪夫 i 型滤波器 n,wc = cheb1ord(wp, wr, ap, ar) 切比雪夫 ii型滤波器 n,wc = cheb2ord(wp, wr, ap, ar) 椭圆滤波器 n,wc = ellipord(wp, wr, ap, ar) (2) 完整的滤波器设计函数，根据滤波器阶数和通带截频设计出实际数字滤波器的 h(z) ( 模拟域到数字域的变换只采用双线性变换法): butter - 巴特沃斯 (butterworth)滤波器设计cheby1 - 切比雪夫 (chebyshev)i 型滤波器设计cheby2 - -切比雪夫 (chebyshev)i

9、i型滤波器设计ellip - 椭圆(elliptic)滤波器设计maxflat - 通用的巴特沃斯 (butterworth)低通滤波器设计yulewalk - yule-walker滤波器设计 (直接数字滤波器设计法 ) butter的格式如下 : b,a = butter(n,wc,options) z,p,k = butter(n,wc,options) a,b,c,d = butter(n,wc,options) n为滤波器的阶数， wc为频率归一化后的3db截频:0 wc 衰耗曲线单调增加。 ,0 解:(1) 确定数字频率指标omega0=0.5*pi; % 通带中心频率omegac

10、1=0.45*pi ; omegac2=0.55*pi; apmax=3; %通带截频omegar2=0.6*pi; armin=10; t=5*pi*10(-6); %取样间隔(2) 利用脉冲响应不变法设计数字滤波器1) 按给定数字频率求得相应的模拟角频率wc1=omegac1/t; wc2=omegac2/t; w0=omega0/t; omegar1=omega0-(omegar2-omega0); wr1=omegar1/t; wr2=omegar2/t; b=wc2-wc1; 2) 将指标转换为归一化模拟低通频率指标norm_wr1=(wr12)-(w02)/(b*wr1); nor

11、m_wr2=(wr22)-(w02)/(b*wr2); norm_wc1=(wc12)-(w02)/(b*wc1); norm_wc2=(wc22)-(w02)/(b*wc2); 3) 设计归一化模拟低通滤波器确定归一化截频if (abs(norm_wr1)-abs(norm_wr2)0 norm_wr=abs(norm_wr2); else norm_wr=abs(norm_wr1); end norm_wc=1; 确定归一化模拟低通滤波器阶数n=buttord(norm_wc,norm_wr,apmax ,armin,s); 设计归一化模拟低通滤波器b_lp,a_lp=butter(n,n

12、orm_wc,s); 4) 把归一化模拟低通滤波器转换为模拟带通滤波器b_bp,a_bp=lp2bp(b_lp,a_lp,w0, b); 5) 利用脉冲响应不变法把模拟滤波器转换为数字滤波器fs=1/t; b,a=impinvar(b_bp,a_bp,fs); w=linspace(0,2*pi,500); h=freqz(b,a,w); plot(w,20*log10(abs(h); axis(0,2*pi,-50,20); grid xlabel(frequency (rad); ylabel(gain (db); sys=tf(b,a,t) (2) 双线性变换法双线性变换法克服了脉冲响应

13、不变法存在频谱混叠的缺点, 但在设计过程中一定要注意，由于宽频带压缩到较窄的频率范围，导致了频率畸变。其基本思想是: 首先按给定的指标设计一个模拟滤波器，通过适当的数学变换方法将s 域映射到 z 域，把无限宽的频带变换成有限宽频带。由于在实现数字化以前已经对频带进行了压缩，所以在数字化以后的频响可以做到无频谱混叠。该法可以用来设计各种不同频带要求的数字滤波器，所以获得广泛应用。例 3: 按 例 2 的指标使用双线性变换法设计一个数字带通滤波器。解: (1)给定数字频率指标omega0=0.5*pi; % 通带中心频率omegac1=0.45*pi; omegac2=0.55*pi; apmax

14、=3; %通带截频omegar2=0.6*pi; armin=10; t=5*pi*10(-6); %取样间隔(2) 用双线性变换法设计滤波器1) 按给定数字频率求得相应的模拟角频率: 频率预畸变omegar1=omega0-(omegar2 - omega0); wc1=(2/t)*tan(omegac1/2); wc2=(2/t)*tan(omegac2/2); wr1=(2/t)*tan(omegar1/2); wr2=(2/t)*tan(omegar2/2); w0=(2/t)*tan(omega0/2); b=wc2-wc1; 2) 将指标转换为归一化模拟低通频率指标norm_wr1

15、=(wr12)-(w02)/(b*wr1); norm_wr2=(wr22)-(w02)/(b*wr2); norm_wc1=(wc12)-(w02)/(b*wc1); norm_wc2=(wc22)-(w02)/(b*wc2); 3) 设计归一化模拟低通滤波器确定归一化截频 : if (abs(norm_wr1)-abs(norm_wr2)0 norm_wr=abs(norm_wr2); else norm_wr=abs(norm_wr1); end norm_wc=1; 确定归一化模拟低通滤波器阶数: n=buttord(norm_wc,norm_wr,apmax,armin,s); 设计

16、归一化模拟低通滤波器: b_lp,a_lp=butter(n,norm_wc,s); 4) 把归一化模拟低通滤波器转换为模拟带通滤波器b_bp,a_bp=lp2bp(b_lp,a_lp,w0,b); 5) 转换为数字滤波器fs=1/t; b,a=bilinear(b_bp,a_bp,fs); w=linspace(0,2*pi,500); h=freqz(b,a,w); plot(w,20*log10(abs(h); axis(0,2*pi,-120,20); grid xlabel(frequency (rad); ylabel(gain (db); sys=tf(b,a,t) 例 4 :

17、若信号由 5hz、15hz和 30hz三个正弦频率成分构成。设计一个椭圆滤波器 滤除 5hz和 30hz频率成分。解: 产生该信号fs=100; t=(0:99)/fs; s1=sin(2*pi*t*5); s2=sin(2*pi*t*15); s3=sin(2*pi*t*30); s=s1+s2+s3; 设计该带通滤波器的通带为10hz,20hz b,a=ellip(4,0.1,40,10 20*2/fs); h,w=freqz(b,a,512); 使用 fft得到滤波前后信号的频谱s=fft(s,512); sf=filter(b,a,s); sf=fft(sf,512); 作图subpl

18、ot(3,2,1); plot(t,s); title(time domain); grid ylabel(original signal); subplot(3,2,5); plot(t,sf); xlabel(time (s); ylabe(filted signal); w1=(0:255)/256*(fs/2); subplot(3,2,2); plot(w1,abs(s(1:256); title (frequency domain); grid ylabel(original spectrum); subplot(3,2,4); plot(w*fs/(2*pi),abs(h); y

19、label(filter); grid; subplot(3,2,6); plot(w1,abs(sf(1:256); xlabel(frequency (hz); ylabel(filted spectrum); grid 设计内容 : 1( 利用双线性变换法设计bw 型及 cb型数字低通滤波器，绘出其幅频特性并进 行比较，其指标是 : (1) 通带:apdb, hz ,0.70,1000 ,(2) 阻带:ardb, =1200hz (bw 型) ,10r ,ardb, =1040hz (cb型) ,10r (3) 取样间隔 t=100s , 2( 设计一个数字带通滤波器，使之满足以下指标:

20、 (1) 通带 0.25,0.35, ,(2) =, =0.4 0.2,lu ,(3) =0.8db, 60db ,pr ,22( 设有一信号 x(n)=1+cosn+cosn ，设计各种 iir 数字滤波器以达到下面的43 目的: ,2, 低通滤波器，滤除cosn 的成分，保留的成分为1+cosn; 34 ,2 高通滤波器，滤除1+cosn 的成分，保留的成分为cosn; 43 ,2, 带通滤波器，滤除1+ cosn 的成分，保留的成分为cosn; 34 ,2 带阻滤波器，滤除cosn 的成分，保留的成分为1+cosn; 43(1) 用 matlab命令 butterord定出滤波器的阶次

21、; 用 butter命令设计设个滤波器 ; 画出滤波器的幅度和相位频响; 取得滤波器的系统函数h(z) (2) 试根据 iir 滤波器设计原理和步骤编写出matlab 程序。 4( 四阶归一化低通巴特沃斯模拟滤波器系统函数为 : 1 ha(s), 432ssss， 2.61313.41422.61311 用双线性变换法从ha(s) 设计四阶带通巴特沃斯数字滤波器h(z) ，并画出 bp ,0.65,0.35,其幅度谱。 ( 设 t,1s ，) 2c1c 5( 采用脉冲响应不变法和双线性变换法设计巴特沃斯数字低通滤波器，满足下列指标: , 通带边缘频率 :0.4 ，通带衰减 :0.5db ; ,

22、 阻带边缘频率 :06 ，阻带衰减 :50db 6( 用切比雪夫型设计一个带通iir 数字滤波器，满足下列指标: , j,0.95|()|1.05,he0|0.25, ， ; j,0|()|0.01,he0.35|0.65, ， ; j,0.95|()|1.05,he0.75|, ，求出滤波器的系统函数h(z) ，并画出幅度响应图。 7( 采用脉冲响应不变法设计一个数字低通滤波器，满足下列指标: f=8khz , 通带边缘频率为 1.5khz，波动为 3db，阻带边缘频率为 2khz ，衰减 s 为 40db,通带等波纹但阻带是单调的，求出滤波器的系统函数h(z) ，并画出幅度 响应图。8(

23、采用脉冲响应不变法设计一个数字高通滤波器，满足下列指标: f=1khz , 通带边缘频率为 2khz ，波动为 3db，阻带边缘频率为1.5khz，衰减 s 为 40db, 单调的通带和阻带，求出滤波器的系统函数h(z) ，并画出幅度响应图。4.2 窗函数的特性分析在滤波器设计中和功率谱估计中，窗函数的选择对设计和分析的结果都起着很重要的作用。截短无穷长的序列会造成吉布斯现象，恰当选取窗函数，可以抑制吉布斯现象。下表给出几种常用窗函数表示式和matlab 实现方法。常用窗函数 matlab 实现方法函数名称 matlab 实现 矩形窗 (rectangular) w=boxcar(n) 或 w

24、=ones(n,1) 汉宁窗 (hanning) w=hanning(n) 或 n=0:n-1 w=1/2*(1-cos(2*pi*n/(n-1) 哈明窗 (hamming) w=hamming(n) 或 n=0:n-1 w=0.54-0.46*cos(2*pi*n/(n-1) 布莱克曼窗 (blackman) w=blackman(n) 或n=0:n-1; w=0.42-0.5*cos(2*pi*n/(n-1)+0.08*cos(4 *pi*n/(n-1) bartlett窗( 三角形窗的一种 ) w=bartlett(n) 或 n=0:n-1; w=1-abs(2*(n-(n-1)/2)/

25、(n-1) w=kaiser(n,beta) 凯塞窗(kaiser) 表中前五种窗函数的形状是固定的，因而一但选择了某种窗函数，用它进行谱分析得到的频谱纹波或设计出的滤波器的阻带衰减就确定了。凯泽窗是一种可调窗， 可以通过改变窗函数的形状来控制频谱纹波或阻带衰减指标，因而获得广泛的应用。matlab 中也提供了许多常用的窗函数，调用格式为: w=boxcar(n) w= hanning(n) w = hamming(n) w = blackman(n) w =kaiser(n,beta)，n是窗函数的长度， beta 是控制 kaiser窗形状的参数。例: 画出布莱克曼窗 ; 解:n=22;

26、n=0:n-1; w=0.42-0.5*cos(2*pi*n/(n-1)+0.08*cos(4*pi*n/(n-1);% w = blackman(n); stem(n,w); 设计内容 : 1( 分析并绘出各种窗函数的时域特性; 2( 使用 fft函数做出各窗函数的频域特性, 并从主瓣宽度和旁瓣相对幅度两个角 度进行比较分析 ; 9,11nn3( 一个序列为 x(n)=cos+cos ，使用 fft分析其频谱 : 2020 (1) 使用不同宽度的矩形窗截短该序列为m点长度 , 取 m分别为 : a)m=20; b)m=40; c)m=160 观察不同长度的窗对谱分析结果的影响; (2 ) 使

27、用哈明窗重做 (1 ) ;对两种窗的结果进行比较和分析; 4( 总结窗的不同长度和不同窗对谱分析结果的影响。思考题 : (1) 在信号谱分析中，如何选择窗函数, (2) 在数字系统设计时，如何选择窗函数, (3) 如何选择不同特性的窗函数, 4.3 fir 数字滤波器设计fir 数字滤波器的优点是具有线性相位。在数据通信、图象处理、语音信号处理等实际应用领域，往往都要求线性相位特性，因而使用fir 数字滤波器。 fir滤波器的传递函数表示为: n,1- nh(z)=h(n) z ,n,0 其单位脉冲响应 h(n) 是有限长的，因此这种滤波器称为有限长单位脉冲响应(fir) 滤波器。其设计过程的

28、核心是求出有限的脉冲响应来逼近给定的频响。由于fir 滤波器脉冲响应在满足h(n)=?h(n-n-1)的对称的条件下具有严格的线性相位，使所设计 的数字滤波器无相位失真，因此在设计中可以着重考虑幅频特性如何符合指标的要 求。常用的设计方法有窗口法与频率取样法。1( 窗口法设计 fir 数字滤波器jw 窗口法的设计思想是基于把给定的频响h(e)通过离散时间傅里叶逆变换 d (idtft) 求得脉冲响应 h(n) ，而后利用窗函数对它进行截短和平滑，达到设计一个 d 物理可实现且具有线性相位fir 数字滤波器的目的。窗口法设计过程为 : 加窗( ，,)idtftdtftjwjw h(e)h(n)

29、h(n)h(e) ,dd jwjwh(e) 是理想滤波器的频域系统函数(频响) ，h(e) 是设计出的实际滤波器的频响。 d 加窗的作用是通过把理想滤波器的无限长脉冲响应h(n) 乘以窗函数来产生一个被 d 截短的脉冲响应h(n) ，即 h(n)= h (n) w (n)并对频响进行平滑。 d 例 1: 分别使用矩形窗函数和哈明窗函数设计一个线性相位低通滤波器: ,- j e | | ,cj h(e) = d, 0 | | , 其中 c=1rad,=12s. c,ccsinc(n,)解: 按 idtft求得其理想脉冲响应 :h(n) = d, , 该理想滤波器的脉冲响应h(n) 以为中心对称。

30、根据线性相位条件h(n)= h(n-n-1) ， d 故实际滤波器的脉冲响应长度为n=25 。分别使用矩形窗函数和哈明窗函数对h(n) 进行截短 : ,ccsinc(n,) , 0nn-1 , 1) 使用矩形窗函数得h(n)= 0 其它2) 使用哈明窗函数得 : ,nccsinc(n,),0.54-0.46cos(2) 0nn-1 ,n,1, h(n)= 0 其它 注意，滤波器阶数m=n-1=24 。%rectangular window加矩形窗m=24;n=0:m;omegac=1; %m: order of the filter hd=omegac/pi.*sinc(omegac/pi*(

31、n-m/2); stem(0:m,hd); h1=hd; h1=fft(h1,512); h1_db=20*log10(abs(h1); %hamming window, 哈明窗w=0.54-0.46.*cos(2*pi*n./m); stem(0:m,w); 加哈明窗h2=hd.*w; h2=fft(h2,512); h2_db=20*log10(abs(h2); 做 fir 滤波器的幅频特性曲线omega=0:pi/255:pi; c=plot(omega,h1_db(1:256),omega,h2_db(1:256),r); 上图绘出了使用矩形窗、哈明窗设计的实际滤波器的幅频特性。使用信

32、号处理工具箱中的fir1函数可以很容易地实现fir 滤波器窗口法设计。 它可设计低通、高通、带通、带阻滤波器及多频带滤波器。格式如下: b = fir1(n,wc) b = fir1(n,wc,ftype ) b = fir1(n,wc,window) b = fir1(n,wc,ftype ,window) b = fir1(n,wc,ftype ,window, noscale ) n 为滤波器的阶数， wc为 3db截频:0 wc 1, 1对应取样频率的一半。对带通 和带阻滤波器， wc=w1,w2，为通带上下截频。若wc=w1,w2 ，w3 wn, 表示滤波器具有 n 个频带:0ww1

33、, w2ww3,wnw1 。ftype 指定滤波器类型 : ftype= high , 指定一个截频为wc的高通滤波器 ; ftype= stop 指定一个带阻滤波器，其阻带截止频率为wc=w1,w2; ftype= dc - 0 在多频带滤波器中，使第一个频带0ww1 为阻带; ftype=dc -1 在多频带滤波器中，使第一个频带0ww1 为通带 ; window 为指定的窗向量，若不指定，默认为哈明窗。使用信号处理工具箱中的fir1函数对本例题进行设计 : %rectangular window b1=fir1(24,1/pi,boxcar(25); h1,w=freqz(b1,1,51

34、2); h1_db=20*log10(abs(h1); %hamming window b2=fir1(24,1/pi,hamming(25); h2,w=freqz(b2,1,512); h2_db=20*log10(abs(h2); c=plot(w,h1_db,w,h2_db,r); 2. 频率取样法 fir 滤波器jw 频率取样法的设计思想是基于对给出的理想频响h(e) 进行取样，通过离散d 傅里叶反变换 (idft) 从频谱样点直接求得有限脉冲响应。频率取样法的过程如下 : 2, 频域取样 (,k),idftdtftjwjwn, h(e) h(k) h(n) h(e) ,d jwjw

35、 其中 h(e) 是理想滤波器的数字域系统函数( 频响)，h(e) 是设计出来的实际 d 滤波器系统函数。频域取样法的关键是正确确定频域取样点，在?0 - 2内的样点有如下约束条件 : h(k)=h(n-k) (k)= - (n -k) 例 2: 采用频率取样法设计fir 滤波器 , 使之满足- j e , |? cj h(e) = 0 , | |? 其中， c = 1rad, = 12s. dc解:n 为奇数且 h (n) 为偶对称，故采用n为奇数的 i 型取样。设计程序如下 : n=25; tao=(n-1)/2; omegac=1; m=fix(omegac/(2*pi/n)+1); o

36、mega=0:n-1*2*pi/n; abs_h =ones(1,m), zeros(1,n-2*m+1),ones(1,m-1); h=abs_h.*exp(-j*tao.*omega); h=ifft(h); figure(1) stem(omega,abs(h); title(fir filter) ylabel(mag of frequency response); xlabel(frequency (rad) hold on; i=0:511*2*pi/512; m1=fix(omegac/(2*pi/512)+1); abs_hd=ones(1,m1),zeros(1,512-2*

37、m1+1),ones(1,m1-1); c=plot(i,abs_hd,r); legend(c,ideal mag.,0); hold off; figure(2) stem(0:n-1,real(h); title(impulse response); xlabel(tume n); ylabel(hn); hh=fft(h, 512) hh_db=20*log10(abs(hh); figure(3) omega=0:255*pi/256; plot(omega,hh_db(1:256); grid set(gca, ytick, (-120:20:-20),-3,0); xlabel(

38、frequency (rad); ylabel(normalized magnitude (db); title(frequency response); 例,: 采用 i 型频率取样法设计 fir 滤波器 , 满足 例 2 所给的技术指标。解: n=25;tao=(n-1)/2;omegac=1; m=fix(omegac/(2*pi/n)+1); omega=0:n-1*2*pi/n; abs_h=ones(1,m),zeros(1,n-2*m+1),ones(1,m-1); h=abs_h.*exp(-j*tao.*omega); stem(omega,abs(h); 图示出理想 fir 低通滤波器的频响和 25 个频率取样点。h=ifft(h); stem(0: