北京交通大学(数字信号处理研究生课程)ch5-4FIR优化设计

1、近代数字信号处理(Advanced Digital Signal Processing) 电子信息工程学院信号与图像处理研究室FIR数字滤波器线性相位FIR数字滤波器的性质窗函数法设计FIR数字滤波器频率取样法设计线性相位FIR数字滤波器线性相位FIR数字滤波器的优化设计FIR数字滤波器的设计误差准则4种类型的线性相位滤波器统一表示离散加权平方误差准则设计FIR积分加权平方误差准则设计FIR等波纹FIR滤波器设计FIR数字滤波器的优化设计常用误差准则所需系统的幅度函数为 D(W)设计所得系统的幅度函数为 A (W)加权误差函数定义 e(W)=W(W)A (W) - D(W)其中W(W)0，是由

2、设计者定义的加权函数。常用误差准则 加权离散Lp 误差 加权积分平方误差 最大最小准则(minimax criterion)或Chebyshev准则4种类型的线性相位滤波器统一表示 为使四种类型的线性相位滤波器能有统一的优化算法，将它们的幅度函数统一表示为 线性方程组的最小均方误差解 x=x1,x2,xnT当mn时，线性方程组的最小均方误差解定义为 可证线性方程组的最小均方误差解的等价方程为 可用MATLAB命令 Ac 获得方程组的解最小均方误差解离散加权平方误差准则设计FIR加权离散平方误差定义为定义误差向量e利用矩阵乘法可得离散加权平方误差准则设计FIR离散加权平方误差准则设计FIR离散加

3、权平方误差准则设计FIRC矩阵的MATLAB实现例：离散加权平方误差准则设计M=63 (II型)，Wp=0.5p，Ws=0.6p的线性相位低通滤波器。M=63;J=(M-1)/2;Wp=0.5*pi;Ws=0.6*pi;Np=100;Ns=100;weightp=1;weights=1; %通带和阻带的加权系数WpN=linspace(0,Wp,Np);WsN=linspace(Ws,pi-0.05,Ns);Wm=WpN WsN;d=ones(1,Np) ,zeros(1,Ns);W=diag(weightp*ones(1,Np),weights*ones(1,Ns);C=cos(Wm*0:J

4、);Q=diag(cos(Wm*0.5);g=(W*Q*C)(W*d);g=g;h=(g(J:-1:2)+g(J+1:-1:3)/4;h=g(J+1)/4 h g(1)*0.5+g(2)/4;h=h fliplr(h);W=linspace(0,pi-0.1,512);H=freqz(h,1,W);plot(W/pi,20*log10(abs(H);00.50.61-80-54-300Normalized frequencyGain, dB00.50.61-110-80-50-200Normalized frequencyGain, dB阻带的加权系数weights=1阻带的加权系数weigh

5、ts=100积分加权平方误差准则设计FIR积分加权平方误差准则设计FIR积分加权平方误差准则设计FIRMATLAB实现 h = firls(M,f,a,w) M: 表示滤波器的阶数。 f: 如FIR滤波器的B个频带分别为 pf1Wpf2， pf3Wpf4， ， pf2B-1Wpf2B 则f是一个有2B个元素的向量，其值为 f=f1, f2，, f2B a: 也是一个有2B个元素的向量，表示滤波器在上述各频带边界的幅度。滤波器各频带的幅度是2个边界点的幅度的线性内插。 w：是B个元素的向量，表示各频带的加权值。缺省时表示各频带的加权值相同。 例：用积分加权平方误差准则设计M=63 (II型)，

6、Wp=0.5p，Ws=0.6p的线性相位低通滤波器。M=63;Fp=0.5; Fs=0.6;h = firls(M, 0 Fp Fs 1, 1 1 0 0);w0=Fp*pi Fs*pi;mag0=freqz(h,1,w0);fprintf(Ap=%f dB As=%f dBn,-20*log10(abs(mag0(1),-20*log10(abs(mag0(2);w=linspace(0,pi,1000);mag=freqz(h,1,w);plot(w/pi,-20*log10(abs(mag); grid on;xlabel(Normalized frequency);00.10.20.3

7、0.40.50.60.70.80.91-80-55-200Normalized frequencyGain, dB 积分加权平方误差准则设计的FIR。离散加权平方误差准则设计的FIRAp= 0.034 dB As= 46.97 dB等波纹FIR滤波器设计 设计方法：利用Chebyshev误差准则设计FIR滤波器 设计出的FIR滤波器的特点： 幅度响应在通带和阻带都会呈现等波波纹的波动 PM(Parks-McClellan ,1972)算法 。 Remez(1957)交换算法。也有人将其称为Remez算法。等波纹FIR滤波器设计 设计准则：确定系数gk,0k J，使得误差达到最小。其中 交替定理

8、 最佳逼近函数G(W)存在并惟一的充要条件为在集合I的范围内至少存在J+2个极值点 W0 W1 WJ+1，使得 e (Wk ) = -e(Wk+1 )， 0 k J, |e (Wk )| = e , 0 k J+1如果J+2个极值点已找到，由交替定理共有J+2个方程，从中可解出J+1个系数gk, 0 k J和误差e 。等波纹FIR滤波器设计 等波纹线性相位FIR滤波器的设计步骤：(1) 用Kaiser提出的经验公式估计滤波器的阶数并由M及类型(I型、II型)确定J。(2) 选定初始极值点Wk; 0kJ+1。等波纹FIR滤波器设计 等波纹线性相位FIR滤波器的设计步骤：等波纹FIR滤波器设计(3

9、)计算e 。可证e 为其中 等波纹线性相位FIR滤波器的设计步骤：等波纹FIR滤波器设计(4)计算函数G(W)在I 中的抽样值。 G(W)在Wk; 0kJ+1点的值为用Lagrange内插获得G(W)在其他点上的值(5) 寻找新的极值点 ，并求出最大误差值 等波纹线性相位FIR滤波器的设计步骤：等波纹FIR滤波器设计(7) 解方程获得gk。(8) 由gk求出hk。等波纹FIR滤波器的MATLAB实现 M,fo,ao,w = firpmord (f,a,dev) (remezord) h = firpm (M,fo,ao,w) (remez)f：如FIR滤波器的B个频带分别为 0Wpf1， pf

10、2Wpf3，pf2(B-1)-2Wpf2(B-1)-1，pf2B-2Wp则f是一个有2B-2个元素的向量(B2)，其值为f=f1, f2,f3, , f2B-2a：是一个B个元素的向量,分别表示FIR滤波器在B个频带中的幅度值。一般对通带取值为1，阻带取值为0。dev：是一个B个元素的向量,分别表示FIR滤波器在B个频带中的波动值。例：设计满足下列指标的等波纹线性相位FIR低通滤波器。 Wp=0.5p，Ws=0.6p, dp=ds=0.0017Fp=0.5;Fs=0.6;ds=0.0017;dp=ds;f=Fp Fs; a=1 0; dev=dp ds;M,fo,ao,w = firpmord

11、(f,a,dev);h = firpm(M+1,fo,ao,w);w0=Fp*pi Fs*pi;mag0=freqz(h,1,w0);fprintf(Ap=%f dB As=%f dBn,-20*log10(abs(mag0(1),-20*log10(abs(mag0(2);w=linspace(0,pi-0.01,1000);mag=freqz(h,1,w);plot(w/pi,20*log10(abs(mag); grid on;xlabel(Normalized frequency);ylabel(Gain, dB);用PM算法实现的MATLAB程序：M=59 设计结果（设计指标为Wp=

12、0.5p，Ws=0.6p, dp=ds=0.0017） Ap= 0.016 dB As= 54.68 dBAp= -20*lg (1-0.0017) =0.0148 dBAs=-20*lg(0.0017) =55.3910例：设计满足下列指标的等波纹线性相位FIR带通滤波器。 Ws1=0.2p，Wp1=0.3p, Wp2=0.6p, Ws1=0.7p, dp=0.1，ds=0.01 Fs1=0.2;Fp1=0.3;Fp2=0.6;Fs2=0.7;f=Fs1 Fp1 Fp2 Fs2;a=0 1 0;Rp=0.1;Rs=0.01;dev=Rs Rp Rs; M,fo,ao,w = firpmord

13、(f,a,dev);h = firpm(M+5,fo,ao,w);w0=Fs1*pi, Fp1*pi, Fp2*pi, Fs2*pi;mag0=freqz(h,1,w0);fprintf(As1=%f Ap1=%f Ap2=%f As2=%f dBn,-20*log10(abs(mag0(1),-20*log10(abs(mag0(2),-20*log10(abs(mag0(3),-20*log10(abs(mag0(4);w=linspace(0.01,pi-0.01,1000);mag=freqz(h,1,w);plot(w/pi,20*log10(abs(mag); grid onxla

14、bel(Normalized frequency);ylabel(Gain, dB);用PM算法实现的MATLAB程序：00.20.30.60.71-40-200Normalized frequencyGain, dBM=25时滤波器的增益响应。 设计结果 M=25 （Ws1=0.2p，Wp1=0.3p, Wp2=0.6p, Ws2=0.7p, dp=0.1，ds=0.01） Ap=-20lg(1-dp)=0.915 dB阻带衰减不满足指标As1=36.3 Ap1=1.4Ap2=1.4 As2=36.3 As=-20lg(ds)=40dB 设计结果 M=30 （Ws1=0.2p，Wp1=0.3

15、p, Wp2=0.6p, Ws1=0.7p, dp=0.1，ds=0.01） 00.20.30.60.71-40-200Normalized frequencyGain, dBM=30，由remez得出的滤波器的增益响应阻带衰减满足指标As1=41.08 Ap1=0.8 Ap2=0.8 As2=41.08Ap=-20lg(1-dp)=0.915 dBAs=-20lg(ds)=40dB问题研讨IIR DF: Ap=1.0 dB As=50.0dB问题研讨FIR DF: Ap=0.98 dB As=49.84 dB% FIR Filter Design and its applicationx,F

16、s,nBit=wavread(E:2013DSPyourn.wav);N=length(x);subplot(2,2,1)plot(x);title(Input signal xk);% Display input signal frequency spectrumw=linspace(0,pi,1000);X=freqz(x,1,w);subplot(2,2,2)plot(w/pi,abs(X); grid on;title(Input signal frequency spectrum X(ejW);问题研讨%Design FIR Digital Low-pass FilterFp=0.7;Fs=0.75;ds=0.003;dp=0.1;f=Fp Fs; a=1 0; dev=dp ds;M,fo,ao,w = firpmord(f,a,dev);h = firpm(M,fo,ao,w);w=linspace(0,pi,1000);H=freqz(h,