FIR数字滤波器的设计--等波纹最佳逼近法

FIR 数字滤波器的设计--等波纹最佳逼近法

一、等波最佳逼近的原理简介

等波纹最佳逼近法是一种优化设计法，即最大误差最小化准则，它克服了窗函数设计法和频率采样法的缺点，使最大误差（即波纹的峰值）最小化，并在整个逼近频段上均匀分布。用等波纹最佳逼近法设计的FIR 数字滤波器的幅频响应在通带和阻带都是等波纹的，而且可以分别控制通带和阻带波纹幅度，这就是等波纹的含义。最佳逼近是指在滤波器长度给定的条件下，使加权误差波纹幅度最小化。与窗函数设计法和频率采样法比较，由于这种设计法使滤波器的最大逼近误差均匀分布，所以设计的滤波器性能价格比最高。阶数相同时，这种设计法使滤波器的最大逼近误差最小，即通带最大衰减最小，阻带最小衰减最大；指标相同时，这种设计法使滤波器阶数最低。

等波纹最佳逼近法的设计思想 。用)(ωd H 表示希望逼近的幅度特性函数，要求设计线性相位FIR 数字滤波器时,)(ωd H 必须满足线性相位约束条件。用()ωH 表示实际设计的滤波器的幅度特性函数。定义加权误差函数()ωε为

()()()()[]ωωωωεH H W d -=

式中，()ωW 为幅度误差加权函数，用来控制不同频带（一般指通带和阻带）的幅度逼近精度。等波纹最佳逼近法的设计在于找到滤波器的系数向量()n h ，使得在通带和阻带内的最大绝对值幅度误差()ωε为最小，这也就是最大误差最小化问题。

二、等波纹逼近法设计滤波器的步骤和函数介绍

1.根据滤波器的设计指标的要求：边界频率，通带最大衰减，阻带最大衰等估计滤波器阶数n ，确定幅度误差加权函数()ωW

2.采用Parks-McClellan 算法，获得所设计滤波器的单位脉冲响应()n h

实现FIR 数字滤波器的等波纹最佳逼近法的MATLAB 信号处理工具函数为firpm 和firpmord 。 firpm 函数采用数值分析中的多重交换迭代算法求解等波纹最佳逼近问题，求的满足等波纹最佳逼近准则的FIR 数字滤波器的单位脉冲响应()n h 。firpmord 根据逼近指标,计算采用Parks-McClellan 算法等波纹最佳逼近滤波器的最低阶数，误差加权向量w,归一化边界频率向量f 。

3对firpm 和firpmord 的说明

firpm

函数功能：采用Parks-McClellan算法设计FIR滤波器

函数格式：hn=firpm(n,f,m,w)

n是滤波器的阶

hn是数字滤波器的单位脉冲响应，其长度为n+1

f是希望滤波器的边界频率向量，要求f是单调增向量，并且从0开始，以1结束，1对于数字频率π

m表示希望滤波器在频率点()k f上m是与f对应的希望滤波器的幅度向量，m和f的长度相等，()k

的幅频响应，m和f给出了希望滤波器的幅度特性。

w i表示对m中第i个频率段幅度逼近精度的加权值。w缺w是误差加权向量，其长度为f的一半。)(

省时，函数默认w全为1，即每个频率段的逼近误差加权值相同。

firpmord函数

函数功能：根据逼近指标，计算采用Parks-McClellan算法等波纹最佳逼近滤波器的最低阶数n，误差加权向量w和归一化边界频率f。其返回参数作为firpm函数的调用参数。

函数格式：[n,f,m,w]=firpmord(f,m,rip,fs)

f可以是归一化边界数字频率向量，也可以是模拟边界频率向量，但必须以0开始，以1结束或fs/2结束，并且其中省略了0和fs/2两个频率点。

fs是时域采样频率，单位Hz。fs缺省时,函数默认fs=2Hz.但这是f的长度（包括省略0和fs/2两个频率点）是m的两倍，即m中的每个元素表示f给定的一个逼近频段上的希望逼近的幅度值。

rip表示f和m描述的各逼近段允许的波纹振幅（幅频响应最大偏差）。

4.注意事项：

省略fs时，f必须是归一化的数字频率

有时计算的阶数n略小，使设计结果达不到指标要求，这时要取n+1或n+2

三、程序

低通滤波器设计

% Lowpass filter: wp=0.4pi, ws=0.6pi, N=26

% peak passband ripple is 0.01, peak stopband ripple is 0.001

clc;clear;clear all;

wp_l=0.4*pi; ws_l=0.6*pi; N1=26; deltal_1=0.01; deltal_2=0.001; %初始化参数

f_l=[0 wp_l ws_l pi]/pi; %设定归一化的firpm函数参数

deltal=[deltal_1,deltal_2];

a_l=[1 1 0 0]; %幅值

[n,f,m,w_l]=firpmord(f_l(2:3),[1 0],deltal); %计算权值

fil_l=firpm(N1,f_l,a_l,w_l); %生成滤波器

figure(1);

stem(fil_l,'.');

title('impulse response of lowpass filter');xlabel('time sequence') ;ylabel('amplitude');

fft_l=fft(fil_l,1024);

db_l=20*log10(abs(fft_l)); %对数表示

figure(2);

plot([0:length(fft_l)-1]/length(fft_l),db_l); %将横坐归化到0-1,单位为1/1024

title(' Log-magnitude response of lowpass filter');xlabel('Normalized frequency');ylabel('Magnitude response (dB)');

axis([0 0.5 -100 10]); %设坐标轴的范围，只取0-2pi的一半

afl=abs(fft_l); %取傅变的幅值

afl_p=zeros(1,1024/2);afl_s=zeros(1,1024/2); %初始化两个一行512列的零序列

afl_p (1:fix((wp_l/(2*pi))*1024))=afl(1:fix((wp_l/(2*pi))*1024))-1; %将通带幅值赋给序列一

afl_s (fix((ws_l/(2*pi))*1024)+3:end)=afl(fix((ws_l/(2*pi))*1024)+3:end/2); %将阻带幅值赋给序列二figure(3);

plot([0:1024/2-1]/1024,afl_p+afl_s); %将通带和阻带的波形显示到一幅图像上

title('passband and stopband approximation error of lowpass filter');xlabel('Normalized frequency ');ylabel('Passband and stopband ripple');

带通滤波器设计

% Bandpass filter:ws1=0.2pi, wp1=0.28pi, wp2=0.72pi, ws2=0.8pi, N=80

% peak passband ripple is 0.01, each peak stopband ripple is 0.001

ws_b1=0.2*pi; wp_b1=0.28*pi; wp_b2=0.72*pi; ws_b2=0.8*pi; N2=80; deltab_1=0.001; deltab_2=0.01; deltab_3=0.001;

deltab=[deltab_1,deltab_2,deltab_3];

f_b=[0 ws_b1 wp_b1 wp_b2 ws_b2 pi]/pi;

a_b=[0 0 1 1 0 0];

[n,f,m,w_b]=firpmord(f_b(2:5),[0 1 0],deltab);

fil_b=firpm(N2,f_b,a_b,w_b);

figure(4);

stem(fil_b,'.');

title('impulse response of bandpass filter');xlabel('time sequence') ;ylabel('amplitude');

fft_b=fft(fil_b,1024);

db_b=20*log10(abs(fft_b));

figure(5);

plot([0:length(fft_b)-1]/length(fft_b),db_b);

title(' Log-magnitude response of bandpass filter');xlabel('Normalized frequency');ylabel('Magnitude response (dB)');

axis([0 0.5 -100 10]);

afb=abs(fft_b);

afb_p=zeros(1,1024/2);afb_s=zeros(1,1024/2);

afb_p(fix((wp_b1/(2*pi))*1024)+3:fix((wp_b2/(2*pi))*1024))=afb(fix((wp_b1/(2*pi))*1024)+3:fix((wp_b2/ (2*pi))*1024))-1;

afb_s (fix((ws_b2/(2*pi))*1024)+2:end)=afb(fix((ws_b2/(2*pi))*1024)+2:end/2);

afb_s (1:fix((ws_b1/(2*pi))*1024))=afb(1:fix((ws_b1/(2*pi))*1024));

figure(6); plot([0:1024/2-1]/(1024),afb_p+afb_s);

title('passband and stopband approximation error of bandpass filter');xlabel('Normalized frequency ');ylabel('Passband and stopband ripple');

高通滤波器设计

% Highpass filter: ws=0.4pi, wp=0.55pi, N=34

% peak passband ripple is 0.01, peak stopband ripple is 0.001

ws_h=0.4*pi; wp_h=0.55*pi; N3=34; deltah_1=0.001; deltah_2=0.01;

deltah=[deltah_1,deltah_2];

f_h=[0 ws_h wp_h pi]/pi;

a_h=[0 0 1 1];

[n,f,m,w_h]=firpmord(f_h(2:3),[0 1],deltah)

fil_h=firpm(N3,f_h,a_h,w_h);

figure(7);

stem(fil_h,'.');

title('impulse response of highpass filter');xlabel('time sequence') ;ylabel('amplitude');

fft_h=fft(fil_h,1024);

db_h=20*log10(abs(fft_h));

figure(8);

plot([0:length(fft_h)-1]/length(fft_h),db_h);

title(' Log-magnitude response of highpass filter');xlabel('Normalized frequency');ylabel('Magnitude response (dB)');

axis([0 0.5 -100 10]);

afh=abs(fft_h);

afh_p=zeros(1,1024/2);afh_s=zeros(1,1024/2);

afh_p (fix((wp_h/(2*pi))*1024)+5:end)=afh(fix((wp_h/(2*pi))*1024)+5:end/2)-1;

afh_s (1:fix((ws_h/(2*pi))*1024))=afh(1:fix((ws_h/(2*pi))*1024));

figure(9);

plot([0:1024/2-1]/(1024),afh_p+afh_s);

title('passband and stopband approximation error of highpass filter');xlabel('Normalized frequency ');ylabel('Passband and stopband ripple');

仿真图

1）以下三幅图依次为低通滤波器的冲激响应，对数幅频响应和通带阻带波纹

2）以下三幅图依次为带通滤波器的冲激响应，对数幅频响应和通带阻带波纹

3）以下三幅图依次为高通滤波器的冲激响应，对数幅频响应和通带阻带波纹