Matlab滤波器设计：FIR滤波器与IIR滤波器设计实现示例

用户1143655

发布于 2022-12-12 08:17:13

6K01

代码可运行

文章被收录于专栏：全栈之殇全栈之殇

运行总次数：1

代码可运行

Matlab滤波器设计：FIR滤波器与IIR滤波器设计实现示例

!! ✨ Matlab版本为R2022b，与以前的版本兼容。本文摘录汇总于：数字滤波器设计实践介绍 - MATLAB & Simulink Example - MathWorks 中国。

本文使用Matlab中的Signal Processing Toolbox中的designfilt函数，并根据频率响应实现如下两种滤波器：

FIR滤波器：有限长单位脉冲响应（Finite Impulse Response）滤波器，又称为非递归型滤波器；
IIR滤波器：无限脉冲响应数字（Infinite Impulse Response）滤波器，又称为递归滤波器

一、FIR滤波器设计

1.1 低通滤波器设计简介

实际FIR滤波器设计通常包括如下三个参数：

Transition width：过渡带宽度；
Peak passband/stopband ripple：最大通带与阻带波纹；
Filter order：滤波器阶数（截断冲激响应长度的选择）。

通常，由于实现理想低通滤波器所需的冲激响应是无限长的，因此无法设计出理想的FIR低通滤波器。理想冲激响应的有限长度逼近会导致滤波器的通带与阻带中均出现波纹，导致带和阻带之间的过渡带宽度非零。

对于滤波器设计，通带/阻带波纹和过渡带宽度作为与理想低通滤波器存在的偏差如下图所示：

FIR滤波器的优点：

其性能非常稳定；
其可以设计成具有线性相位的滤波器。

但是，FIR滤波器仍然可能存在长瞬时响应，对于某些问题会导致计算成本的增加。

1.2 最小阶FIR滤波器设计

最小阶FIR滤波器设计主要包括如下两步：

通过指定通带/阻带的频率、通带波纹、阻带衰减，我们就可以获得最小阶FIR滤波器设计。
设计算法会选择符合设定的FIR滤波器最小滤波器长度。

Matlab代码如下所示，通过designfilt函数实现，并通过fvtool函数查看设计的低通通FIR滤波器的幅值响应：

Fpass = 0.3;  % 通带频率系数
Fstop = 0.43; % 阻带频率系数
Ap = 1;       % 通带波纹
Ast = 30;     % 阻带衰减

% 通过designfilt函数实现FIR低通滤波器的设计
d = designfilt('lowpassfir', ...    % 低通滤波器
    'PassbandFrequency', Fpass, ... % 通带频率
    'StopbandFrequency', Fstop, ... % 阻带频率
    'PassbandRipple', Ap, ...       % 通带波纹
    'StopbandAttenuation', Ast);    % 阻带衰减

% 通过fvtool函数设计的低通通FIR滤波器的幅值响应
hfvt = fvtool(d);

代码执行结果如下所示：

!! 通过filtord函数可以查看生成的滤波器的阶数，如下图所示，设计的滤波器d的阶数目为19：

注意： 另外,fvtool函数会打开滤波器可视化工具，通过该工具箱我们可以查看设计的滤波器如下各种分析结果：

!! designfilt函数默认选择一个等波纹（线性相位）设计算法Equiripple，可以通过命令info(d)查看滤波器的信息得到滤波器设计算法。通常，线性相位波纹滤波器可以取得理想的结果，这是由于对于给定阶数，线性行为滤波器与理想滤波器的最大可能偏差最小。

🚀 注意： 通常，我们还可以使用Kaiser窗获得最小阶FIR滤波器设计。虽然Kaiser窗通常会得到更大的滤波器阶数，但是该算法的计算成本更低，并且不太可能出现收敛的问题。

下面的代码使用Kaiser窗方法设计与上面相同的滤波器，并与等波纹滤波器（Equiripple算法）的幅值响应进行比较：

% 通过Kaiser窗实现FIR低通滤波器的设计
dk = designfilt('lowpassfir', ...   % 低通滤波器
    'PassbandFrequency', Fpass, ... % 通带频率
    'StopbandFrequency', Fstop, ... % 阻带频率
    'PassbandRipple', Ap, ...       % 通带波纹
    'StopbandAttenuation', Ast, ... % 阻带衰减
    'DesignMethod', 'kaiserwin');   % 使用Kaiser窗方法设计滤波器

% 即将滤波器dk添加到绘图结果中
addfilter(hfvt, dk);
% 添加图例
legend(hfvt, 'Equiripple滤波器设计算法', 'Kaiser窗口滤波器设计算法')

代码执行结果如下图所示：

通过filtord(dk)可以查看滤波器阶数为24：

1.3 以赫兹为单位设置滤波器频率参数

Fpass = 350; % 通带频率
Fstop = 400; % 阻带频率
Ap = 1;      % 通带波纹
Ast = 30;    % 阻带衰减
Fs = 2000;   % 采样频率

% 通过designfilt函数实现FIR低通滤波器的设计
d = designfilt('lowpassfir', ...    % 低通滤波器
    'PassbandFrequency', Fpass, ... % 通带频率
    'StopbandFrequency', Fstop, ... % 阻带频率
    'PassbandRipple', Ap, ...       % 通带波纹
    'StopbandAttenuation', Ast, ... % 阻带衰减
    'SampleRate', Fs);              % 采样频率

% 通过fvtool函数设计的低通通FIR滤波器的幅值响应
hfvt = fvtool(d);

代码执行结果如下图所示：

1.4 固定阶数、固定过渡带带宽的FIR滤波器设计

本小节采用等波纹滤波器算法（Equiripple算法）与最小二乘法（Least square算法）两种滤波器设计算法，设计具有固定阶数、固定过渡带带宽的FIR滤波器。

Matlab代码如下所示：

N = 30;      % 滤波器阶数
Fpass = 400; % 通带频率
Fstop = 450; % 阻带频率
Fs = 2000;   % 采样频率

% 使用默认的equiripple滤波器设计算法设计FIR滤波器deq
deq = designfilt('lowpassfir', ...  % 低通FIR滤波器
    'FilterOrder', N, ...           % 滤波器阶数
    'PassbandFrequency', Fpass, ... % 通带频率
    'StopbandFrequency', Fstop, ... % 阻带频率
    'SampleRate', Fs);              % 采样率

% 使用最小二乘法ls滤波器设计算法设计FIR滤波器deq
dls = designfilt('lowpassfir', ...  % 低通FIR滤波器
    'FilterOrder', N, ...           % 滤波器阶数
    'PassbandFrequency', Fpass, ... % 通带频率
    'StopbandFrequency', Fstop, ... % 阻带频率
    'SampleRate', Fs, ...           % 采样率
    'DesignMethod', 'ls');          % 最小二乘法滤波器设计算法

% 滤波器可视化
hfvt = fvtool(deq, dls);
legend(hfvt, 'Equiripple滤波器设计算法', '最小二乘法滤波器设计算法');

代码执行结果如下图所示：

!! ✨ 说明：

等波纹滤波器非常适合与满足特定容差的情况，比如设计固定最小阻带、衰减的滤波器；但是该方法对于想要最小化通带/阻带中的误差能量的问题，处理结果则通常不太理想。
如果我们想尽可能减少某个频带内信号的能量，请使用最小二乘法设计算法。

1.5 固定阶数、固定截止频率的FIR滤波器设计

我们可以采用窗口设计方法来设计具有固定滤波器阶数和截止频率的滤波器，即我们可以通过使用不同窗口（比如，Hamming、Chebyshev窗口等）来控制阻带衰减，并保持阶数不变。

Matlab代码如下所示：

% Hamming窗口FIR滤波器设计
dhamming = designfilt('lowpassfir', ... % 低通FIR滤波器
    'FilterOrder', 50, ...              % 滤波器阶数
    'CutoffFrequency', 100, ...         % 截止频率
    'SampleRate', 1000, ...             % 采样频率
    'Window', 'hamming');               % Hamming窗口函数

% Chebyshev窗口FIR滤波器设计
dchebwin = designfilt('lowpassfir', ... % 低通FIR滤波器
    'FilterOrder', 50, ...              % 滤波器阶数
    'CutoffFrequency', 100, ...         % 截止频率
    'SampleRate', 1000, ...             % 采样频率
    'Window', {'chebwin', 90});         % 旁瓣衰减为90的Chebyshev窗口函数

% 滤波器可视化及其比较
hfvt = fvtool(dhamming, dchebwin);
legend(hfvt, 'Hamming窗口', 'Chebyshev窗口')

代码执行结果如下图所示：

二、IIR滤波器设计

2.1 IIR滤波器简介

FIR滤波器的缺点主要体现在： 其需要很大的滤波器阶数才能满足实际设计，因此增加对计算的需求，此时可以使用IIR滤波器解决该问题。

!! IIR滤波器设计思想： 如果波纹保持不变，滤波器阶数与过滤带宽度成反比。通过反馈，使用很小的滤波器阶数就可以设计满足需求的滤波器。无限脉冲激励响应IIR（Infinite Impulse Response）的内涵为：当冲激施加到滤波器时，输出永远不会衰减到零。

当计算资源非常宝贵时，IIR 滤波器非常有用。然而，稳定的因果 IIR 滤波器无法提供完美的线性相位。在要求相位线性的情况下，避免使用 IIR 设计。使用 IIR 滤波器的另一个重要原因是相对于 FIR 滤波器，IIR 滤波器的群延迟较小，从而瞬时响应更短。

常用的IIR滤波器包括：Butterworth滤波器、Chebyshev I 类滤波器、Chebyshev I 类滤波器以及椭圆滤波器。他们的主要特点及其异同如下表所示：

IIR滤波器	特点
Butterworth滤波器	一种具有最大平坦度的IIR滤波器，但是这会导致过渡带的宽度增加，并需要较大的滤波器阶数才能减小过渡带的宽度。
Chebyshev I类滤波器	通过允许通带波纹，Chebyshev I类滤波器的过渡带宽度小于同阶Butterworth滤波器。Butterworth 和 Chebyshev I 类滤波器都具有最平坦的阻带。对于给定的滤波器阶数，需要在通带波纹和过渡带宽度之间进行权衡。
Chebyshev II类滤波器	其具有最平坦的通带和等波纹阻带。由于实际情况中，通常不需要非常大的衰减，因此可以允许一些阻带波纹，以换来使用较小的阶数，来获得满足需求的过渡带宽。
椭圆滤波器	其通过允许通带和阻带中波纹来泛化Butterworth与Chebyshev滤波器。随着波纹的减小，椭圆滤波器可以逼近任意Butterworth与Chebyshev滤波器的幅值和相位响应。

2.2 IIR滤波器的实现及性能比较

针对示例4，分别使用Butterworth滤波器、Chebyshev I 类滤波器、Chebyshev I 类滤波器和椭圆滤波器四种IIR滤波器设计方法进行设计。并分别比较四种方法的如下三个方面的滤波器性能参数：

滤波器阶数；
滤波器响应；
群延迟比较；

（1）四种滤波器的实现

Matlab代码如下所示：

Fp = 100;   % 通带频率
Fst = 250;  % 阻带频率
Ap = 1;     % 通带波纹
Ast = 60;   % 阻带衰减
Fs = 2e3;   % 采样频率

% Butterworth 滤波器
dbutter = designfilt('lowpassiir', ...  % IIR低通滤波器
    'PassbandFrequency', Fp, ...        % 通带频率
    'StopbandFrequency', Fst, ...       % 阻带频率
    'PassbandRipple', Ap, ...           % 通带波纹
    'StopbandAttenuation', Ast, ...     % 阻带衰减
    'SampleRate', Fs, ...               % 采样频率
    'DesignMethod', 'butter');          % Butterworth滤波器

% Chebyshev I类滤波器
dcheby1 = designfilt('lowpassiir', ...  % IIR低通滤波器
    'PassbandFrequency', Fp, ...        % 通带频率
    'StopbandFrequency', Fst, ...       % 阻带频率
    'PassbandRipple', Ap, ...           % 通带波纹
    'StopbandAttenuation', Ast, ...     % 阻带衰减
    'SampleRate', Fs, ...               % 采样频率
    'DesignMethod', 'cheby1');          % Chebyshev I类滤波器

% Chebyshev II类滤波器
dcheby2 = designfilt('lowpassiir', ...  % IIR低通滤波器
    'PassbandFrequency', Fp, ...        % 通带频率
    'StopbandFrequency', Fst, ...       % 阻带频率
    'PassbandRipple', Ap, ...           % 通带波纹
    'StopbandAttenuation', Ast, ...     % 阻带衰减
    'SampleRate', Fs, ...               % 采样频率
    'DesignMethod', 'cheby2');          % Chebyshev II类滤波器

% 椭圆滤波器
dellip = designfilt('lowpassiir', ...   % IIR低通滤波器
    'PassbandFrequency', Fp, ...        % 通带频率
    'StopbandFrequency', Fst, ...       % 阻带频率
    'PassbandRipple', Ap, ...           % 通带波纹
    'StopbandAttenuation', Ast, ...     % 阻带衰减
    'SampleRate', Fs, ...               % 采样频率
    'DesignMethod', 'ellip');           % 椭圆滤波器

（2）四种滤波器的阶数比较

Matlab代码如下所示：

% 比较四个 IIR 滤波器的阶数
FilterOrders = [filtord(dbutter) filtord(dcheby1) filtord(dcheby2) filtord(dellip)]

代码执行结果如下图所示：

有代码执行结果可以看出：对于相同的设定约束，Butterworth 方法产生最高阶，椭圆方法产生最小阶。

（3）四种滤波器的响应比较

Matlab代码如下所示：

% 比较四个 IIR 滤波器的响应
hfvt = fvtool(dbutter,dcheby1,dcheby2,dellip);
axis([0 1e3 -80 2]);

代码执行结果如下图所示：

（4）四种滤波器的群延迟比较

对于 IIR 滤波器，我们不仅需要考虑波纹/过渡带宽度的权衡，还需要考虑相位失真的程度。 通常，由于无法在整个奈奎斯特区间内都有线性相位。因此，可以查看相位响应离线性有多远，具体实现思路为：观察（理想情况下为常量的）群延迟，查看它的平坦度。

Matlab代码如下所示：

% 比较四种滤波器的群延迟
hfvt = fvtool(dbutter,dcheby1,dcheby2,dellip,'Analysis','grpdelay');
legend(hfvt,'Butterworth滤波器', 'Chebyshev Type I滤波器', 'Chebyshev Type II滤波器','椭圆滤波器')