【DSP】数字滤波器设计（所有设计方式大全，含matlab代码分析）

对于离散时间系统(或数字滤波器)而言，其系统函数（传递函数）通常表示为 $H(z)$ 。
当我们把复平面上的 z 变量限制在单位圆 $z = e^{j\omega}$ 上时，便得到该系统在频域的响应函数，即频率响应：
$H\bigl(e^{j\omega}\bigr) = \sum_{n=0}^{N} h(n)\, e^{-j \omega n}$
这里， $h(n)$ 是滤波器的单位脉冲响应。
频率响应可以写成幅度和相位的形式：
$H\bigl(e^{j\omega}\bigr) = |H\bigl(e^{j\omega}\bigr)| \, e^{\,j\,\phi(\omega)}$
其中
$H\bigl(e^{j\omega}\bigr)| = \sqrt{\bigl[\Re(H(e^{j\omega}))\bigr]^2 + \bigl[\Im(H(e^{j\omega}))\bigr]^2}, \quad \phi(\omega) = \arg\bigl(H(e^{j\omega})\bigr)$

1.2 与 s 平面、z 平面的关系

在模拟滤波器中，我们通常考察 $H(s)$ 在 $s = j\omega$ 轴上的响应。
在数字滤波器中，常用 z-变换来描述系统。当 $z = e^{j\omega}$ （即单位圆）时，对应的是离散时间傅里叶变换 (DTFT)。
因此，只要在 z-平面上把 z 取到单位圆，便可得到滤波器的频率响应。

（2）数字滤波器在 z-平面上的极零分布与频率响应

2.1 一般形式

假设数字滤波器的系统函数为有理函数形式，即

$H(z) = \frac{B(z)}{A(z)} = \frac{b_0 + b_1 z^{-1} + \cdots + b_M z^{-M}}{1 + a_1 z^{-1} + \cdots + a_N z^{-N}}$

当我们把 z 替换为 $e^{j\omega}$ 时，就得到频率响应：

$H\bigl(e^{j\omega}\bigr) = \frac{B\bigl(e^{j\omega}\bigr)}{A\bigl(e^{j\omega}\bigr)}$

2.2 以极零形式表达

如果把分子与分母分别因式分解，令

$B(z) = b_0 \prod_{m=1}^{M} (1 - z_m z^{-1}), \quad A(z) = \prod_{n=1}^{N} (1 - p_n z^{-1})$

其中 $z_m$ 表示零点（zero）， $p_n$ 表示极点（pole）。那么，

$H(z) = \frac{b_0 \prod_{m=1}^{M} (1 - z_m z^{-1})}{\prod_{n=1}^{N} (1 - p_n z^{-1})}$

将 $z = e^{j\omega}$ 代入，即得到式 (3-1-20) 所示的数字滤波器频率响应形式：

$H\bigl(e^{j\omega}\bigr) = b_0 \,\frac{\displaystyle\prod_{m=1}^{M} \bigl(1 - z_m e^{-j\omega}\bigr)}{\displaystyle\prod_{n=1}^{N} \bigl(1 - p_n e^{-j\omega}\bigr)}$

（3）数字滤波器的幅度特性与相位特性

3.1 幅度特性

如果记

$H\bigl(e^{j\omega}\bigr) = \frac{\displaystyle \prod_{m=1}^{M} \bigl(e^{-j\theta_m} - z_m\bigr)}{\displaystyle \prod_{n=1}^{N} \bigl(e^{-j\psi_n} - p_n\bigr)}$

在图中的推导里，会将 $z_m$ 和 $p_n$ 写成其极坐标或指数形式，如

$z_m = B_m e^{\,j\,\theta_m}, \quad p_n = A_n e^{\,j\,\psi_n}$

则频率响应的幅度可写成

$|H(e^{j\omega})| = \frac{\prod_{m=1}^{M} |e^{-j\omega} - z_m|}{\prod_{n=1}^{N} |e^{-j\omega} - p_n|}$

在教科书中常见的写法（如式(3-1-22)）为：

$|H(e^{j\omega})| = \frac{\prod_{m=1}^{M} |B_m|}{\prod_{n=1}^{N} |A_n|}$

当然，具体形式要看零极点如何表达距离，常见的推导会把单位圆上的点与零极点间的几何距离写出来，然后再简化成上式(或其变体)。如果零极点都在单位圆内部或外部，就会有相应的简化。

3.2 相位特性

相位特性是频率响应的辐角（argument），在图中对应式 (3-1-23)：
$\phi(\omega) = \arg\bigl(H(e^{j\omega})\bigr) = \sum_{m=1}^{M} \theta_m(\omega) \;-\; \sum_{n=1}^{N} \psi_n(\omega)$
其中 $\theta_m(\omega)$ 和 $\psi_n(\omega)$ 分别是由零点、极点对频率 $\omega$ 产生的相位贡献。
当我们在单位圆上考察时， $e^{j\omega}$ 与每个零点、极点的夹角都会对最终的相位产生累加/累减作用。

3、数字滤波器的分类

数字滤波器可以从多个角度进行分类：

（1）按频率特性分类

低通滤波器：允许低频信号通过，抑制高频信号。
高通滤波器：允许高频信号通过，抑制低频信号。
带通滤波器：仅允许某一频率范围内的信号通过。
带阻滤波器：抑制某一频率范围内的信号，通过其他频率信号。

（2）按实现方式分类

按实现方式的不同,数字滤波器可以分为非递归型(Nonrecursive)数字滤波器和递归型(Recursive)数字滤波器两大类。

2.1 非递归型数字滤波器

定义
非递归型滤波器在实现上不包含反馈回路（feedback），即滤波器输出信号 $y(n)$ 不会通过系统再“反馈”到滤波器输入端参与新的计算。数学形式上，这类滤波器通常可以用有限长卷积或有限长加权和来表示：
$y(n) = \sum_{k=0}^{M} b_k x(n-k)$
其中， $x(n)$ 为输入信号， $y(n)$ 为输出信号， $b_k$ 为滤波器系数，M 为有限阶次。
特点
1. 系统方程中只包含输入信号及其有限次延时项，不包含输出信号的延时项。
2. 相对实现较简单，运算稳定性好，一般不会出现滤波器“发散”现象。
3. 由于没有输出的反馈，若要实现较陡峭的滤波特性，往往需要较高的滤波器阶数（更多系数）。

2.2 递归型数字滤波器

定义
递归型滤波器则在实现上包含反馈回路，即输出 $y(n)$ 的当前值不仅与输入信号相关，还与先前若干个输出值相关。其典型的差分方程形式可以写为：
$y(n) = \sum_{k=0}^{M} b_k x(n-k) \;+\; \sum_{k=1}^{N} a_k y(n-k)$
可以看到，滤波器的输出不但依赖输入信号，也依赖于过去的输出值。
特点
1. 由于包含反馈，滤波器结构相对更复杂。
2. 可以用较低的阶数（即较少系数）实现陡峭的幅频特性，效率更高。
3. 可能存在稳定性问题，需要在设计时确保滤波器极点都在单位圆内。
4. 如果系数设计或实现不当，存在出现振荡甚至发散的风险。

补充说明

非递归型滤波器通常也被称为“直接型”或“有限长单位冲激响应型”结构，但在一些文献中会将“非递归型(FIR)”和“递归型(IIR)”直接对应起来。然而严格来说，“FIR”和“IIR”更多是从滤波器脉冲响应的长度来划分（见下文），而“非递归型”“递归型”是从实现方式（是否存在反馈）来划分。

一般情况下，FIR 滤波器往往通过非递归方式实现，但也可以通过递归结构实现（比如某些多相结构、线性相滤波器的特殊实现）；同理，IIR 滤波器通常通过递归结构实现，但理论上也可以通过特定的级联或并联形式将其做成表面上“非递归”的结构。不过这些都属于比较高级、特殊的实现方式，初学者往往先理解最常见的对应关系即可。

（3）按脉冲响应特性分类

按数字滤波器对脉冲响应的特性来分类,数字滤波器可以分为有限脉冲响应FIR(FiniteImpulse Response)数字滤波器和无限脉冲响应 [IR(Infinite Impulse Response)数字滤波器。

FIR滤波器的传递函数只有零点，不含极点，它的单位脉冲响应 h(k) 只包含有限个非零值，即这种数字滤波器的脉冲响应是时间有限的。

IIR滤波器既有零点又有极点，它的脉冲响应 h(k) 中含有无限多个非零值，即这种滤波器的脉冲响应是无限长时间序列，在一定的时间后可能变小，但不会为零。

一般情况下，FIR滤波器比较适合用非递归型方式来实现，而IIR滤波器比较适合用递归型方式来实现。但无论是FIR滤波器还是IR滤波器，都可用两种方法中的任一种来实现，只是按上述方法更简便而已。

3.1 有限脉冲响应（FIR）滤波器

定义
当数字滤波器对单位脉冲输入的响应（脉冲响应 $h(n)$ ）在有限长度内会完全衰减至零，或仅在有限个采样点上非零、之后保持零，则称该滤波器为有限脉冲响应滤波器。
特点
1. $h(n)$ 在时间域（或序列）上是有限长的，意味着从差分方程角度看，没有依赖过去输出的反馈项。
2. 常见结构是非递归型，因此容易实现线性相位，稳定性好。
3. 对同样的滤波指标，可能需要的滤波器阶数（系数个数）比 IIR 更高，计算量更大。
4. 设计方法较为系统，如窗函数法、频率抽样法、最小平方误差法等。

3.2 无限脉冲响应（IIR）滤波器

定义
当数字滤波器对单位脉冲输入的响应 $h(n)$ 理论上无限延长（即在 $n \to \infty$ 的过程中依然可能非零），则称为无限脉冲响应滤波器。
特点
1. 一般实现上包含反馈，因此是递归型结构居多。
2. 可以在较低阶数下实现较陡峭的滤波特性，计算效率更高。
3. 需要注意滤波器的稳定性，极点必须位于单位圆内，否则输出可能会发散。
4. 设计多借鉴模拟滤波器设计思想，如双线性变换法、脉冲不变法等，从经典模拟滤波器（Butterworth、Chebyshev 等）映射得到数字滤波器。

补充说明

FIR 和 IIR 的核心区别在于脉冲响应长度：FIR 的冲激响应是有限长度，而 IIR 的冲激响应是无限长度。

在实际工程中，FIR 滤波器多用于对相位特性有较高要求、或必须保证严格线性相位的场景；IIR 滤波器则常用于对实时性或运算资源较为敏感、又需要较好幅频特性的场合。

即使 FIR 滤波器的冲激响应是有限长，理论上也可以通过某些方法（例如多相分解或其他技巧）做成表面上“带反馈”结构；同样 IIR 也可通过某些级联/并联等手段改写。但通常最直观的对应仍是：FIR→非递归结构、IIR→递归结构。

4、数字滤波器的构成

一个线性时不变（LTI）的离散时间数字滤波器通常可以用其系统函数 $H(z)$ 表征。对于因果的 IIR（无限长单位冲激响应）滤波器，常见的差分方程形式为：

$y(k) = \sum_{i=0}^{M} b_i\,x(k-i) \;-\; \sum_{j=1}^{N} a_j\,y(k-j)$

其中：

$x(k)$ 为输入信号；
$y(k)$ 为输出信号；
$b_i$ 和 $a_j$ 分别为滤波器的前馈（FIR部分）和反馈（IIR部分）系数。

在 z-域中，上述滤波器的系统函数可写为

$H(z) \;=\; \frac{Y(z)}{X(z)} \;=\; \frac{\sum_{i=0}^{M} b_i\,z^{-i}}{1 + \sum_{j=1}^{N} a_j\,z^{-j}}$

数字滤波器的“结构”或“构成”指的是如何将这些数学关系转化为由延时单元（ $z^{-1}$ ）、加法器和乘法器组成的实际电路或程序实现方式。下文介绍的几种主要结构是最常用的实现方法。

（1）直接构成法

1.1 直接构成法 I

“直接构成法 I”（Direct Form I）是最直观的实现方式，直接根据差分方程中的前馈部分和反馈部分各自使用一条延时线（即多个串联的 $z^{-1}$ ），分别实现：

前馈（FIR）部分：对输入信号 $x(k)$ 进行延时并乘以各自系数 $b_i$ ，然后加和得到某个中间信号；
反馈（IIR）部分：对输出信号 $y(k)$ 的延时量乘以各自系数 $a_j$ ，再与前馈部分相减（或加）得到当前输出 $y(k)$ 。

在图 3-1-4(a)（“直接构成法 I 型”）中可以看到：

左侧有一组延时单元 $z^{-1}$ 依次对输入 $x(k)$ 延时，输出相应的“抽头”信号，再乘以系数 $b_i$ 后相加。
右侧有一组延时单元对输出 $y(k)$ 进行延时，乘以系数 $a_j$ ，与前馈部分合并到输出端。

特点：

思路清晰，直接按照差分方程“前馈—反馈”两个部分分别搭建；
延时单元个数较多，前馈与反馈分别需要 $\max(M,\,N)$ 个左右（具体还看是否进行合并）。

1.2 直接构成法 II

“直接构成法 II”（Direct Form II）常被称作“共用延时线”或“转置结构”。它的关键在于合并前馈与反馈部分所需的延时单元，从而减少硬件或存储单元的需求。

通过一些代数变换（可以参考教材中 3-1-3、3-1-4 的推导），把前馈和反馈的延时操作合并到同一条延时线中；
一般可把滤波器差分方程转换为一种“状态空间”形式，然后将中间状态与输出、输入的乘积相加得到结果。

在图 3-1-4(b)（“直接构成法 II 型”）中可以看到：

仅有一条延时线（若为最高阶 $\max(M,N)$ 需要相应数量的串联延时）；
系数 $b_i$ 和 $a_i$ 分布在这条延时线的不同位置，通过加法器与输入、输出相结合，完成滤波功能。

特点：

相比直接构成法 I，所需的延时单元数通常更少（为 $\max(M, N)$ 个，而 Direct Form I 可能需要 M+N 个）；
但是对系数量化误差和有限字长效应相对更敏感，在实际实现中需要权衡。

（2）间接构成法

直接构成法简单、直观，但也有局限性，特别是在高阶数字滤波器中，各系数的变化对其频率响应影响比较大，在精度达不到要求的情况下，很难保证滤波器的性能。所以，在实际中经常将高阶数字滤波器分解成一系列低阶数字滤波器，然后再按一定的规则将它们组合起来。

2.1 级联构成法

“级联构成法”也称为“串联构成法”，是将一个高阶数字滤波器的系统函数 $H(z)$ 因式分解成若干个低阶（通常是一阶或二阶）子滤波器的乘积，然后将这些子滤波器串联起来。

理论依据

若

$H(z) = \prod_{i=1}^{L} H_i(z)$

则可以分别实现每个 $H_i(z)$ ，再将它们的输出依次作为下一个的输入，从而得到整个滤波器的输出。常见做法是：

对 IIR 滤波器而言，将高阶多项式按极点-零点配对，拆成若干二阶段（二阶节），有时还会有一阶段；
每个二阶段用直接构成法 I 或 II 等方式实现；
通过级联连接各段，形成最终的滤波器。

优势与特点

数值稳定性：将高阶滤波器拆分为多个二阶或一阶滤波器，可以在每个子滤波器处更好地控制极点位置、量化误差和稳定性；
灵活性：可以针对每个二阶节选用最优的结构（例如 Direct Form II 或其他形式），也更易于系数调整和实现；
模块化设计：每个二阶节都是一个可复用的基本模块，有助于硬件或软件的模块化实现。

2.2 并联构成法

“并联构成法”是将数字滤波器的系统函数 $H(z)$ 分解成若干子滤波器的和，然后将这些子滤波器并联，最后把并联各支路的输出相加。常用做法是部分分式分解（Partial Fraction Expansion）：

$H(z) = \sum_{i=1}^{L} H_i(z)$

每个 $H_i(z)$ 可能是低阶的传递函数（同样是一阶或二阶），分别用简单的直接构成法实现，然后把这些输出相加，就得到滤波器最终输出。

特点

对于有理函数型滤波器（即 IIR），可以将分子分母做部分分式分解，分解成若干形如
$\frac{B_i(z)}{A_i(z)}$
的子滤波器；
这种方式适合某些特定系数分布或特殊滤波器的实现，也可以缓解某些系数量化问题；
与级联类似，也是一种模块化思路，只是从乘法分解变成了加法分解。

（3）格形结构

格形滤波器（Lattice Filter） 是一种通过级联“格状单元”实现滤波功能的结构形式。每个“格状单元”由加法器、延迟单元和乘法器构成，且具有“前向”和“后向”两个信号通道（或称正向预测误差、反向预测误差），两路信号通过一系列反射系数（Reflection Coefficients）或称“格系数”进行组合。格形结构最初多用于线性预测编码（LPC）和自适应滤波器中，因其数值稳定性好、可逐级更新系数而广受关注。

概括来说，格形滤波器有以下几个关键点：

模块化设计：高阶滤波器通过多个低阶“格单元”逐级构成。
反射系数：每个格单元用一个或一对系数来控制信号的前向与后向传输。
数值稳定性：在固定点实现中，格形结构对量化误差更不敏感；在自适应算法中，格系数的更新往往也比较稳定。

在数字滤波器中，格形结构可分为 FIR 全零点、IIR 全极点 以及 IIR 零极点 三大类。下面分别说明其原理与实现方式。

3.1 FIR型全零点系统的格形结构

系统函数与基本结构

对于一个 FIR 型全零点系统，假设其系统函数为

$H(z) = 1 + \sum_{i=1}^{M} b(i)\,z^{-i}$

它的格形实现如图 3-1-7 所示（即题主给出的示意图）。从结构上看，每一级格单元都由以下组成：

一路正向通道（前向信号）；
一路反向通道（后向信号）；
乘法器（系数 $k_i$ 或类似符号）；
加法器和减法器。

在实现中，可以将滤波器系数 $b(i)$ 转化为格系数 $k_1, k_2, \dots, k_M$ 。该转换有一定的数学推导过程（如线性预测或最小二乘逼近时的求解），但一旦确定了这些格系数，滤波器就可以用模块化的格形单元来级联实现。

特点与分析

有限脉冲响应（FIR）：由于没有反馈回路，系统的单位脉冲响应在有限长度后即为零。
线性相位可控：FIR 滤波器易于设计成线性相位，而格形结构在实现上往往可以减少系数敏感度。
数值稳定性：格形结构中的每一级都是较低阶的滤波单元，不会产生大范围的累积量化误差；在硬件实现中有利于降低对系数量化精度的要求。

3.2 IIR型全极点系统的格形结构

系统函数与基本结构

若滤波器是 IIR 型全极点系统，其传递函数可写为

$H(z) = \frac{1}{A(z)} = \frac{1}{1 + \sum_{i=1}^{M} a(i)\,z^{-i}}$

图 3-1-8 显示了其格形结构示意图。与 FIR 格形类似，这种结构依旧保留了“前向通道”和“后向通道”的思路，但此时每个格单元会有反馈路径以实现极点特性。

推导思路

极点实现：IIR 滤波器具有极点，意味着滤波器差分方程中存在对输出的反馈项。格形单元通过“前向-后向”误差和反射系数的组合，等效实现了对输出信号的递归反馈。
分级分解：高阶多项式 $A(z)$ 可以通过一系列低阶分式或反射系数的组合来表达。这一点与“级联”或“并联”结构的思路类似，但格形结构具有更好的数值稳定性。

特点与分析

适用于递归系统：能在较低阶数下实现较陡峭的滤波特性；
数值稳定性优于直接型结构：因为每个格单元只负责部分极点参数，量化误差不会大范围扩散；
自适应滤波优势：在自适应滤波中，更新格系数比更新直接型滤波器的所有系数更稳定，容易收敛。

3.3 IIR型零极点系统的格形结构

系统函数与结构示意

对于同时包含零点和极点的 IIR 系统，传递函数一般可以写成

$H(z) = \frac{B(z)}{A(z)} = \frac{1 + \sum_{i=1}^{M} b(i)\,z^{-i}}{1 + \sum_{i=1}^{M} a(i)\,z^{-i}}$

图 3-1-9（题主给出的示意）则展现了零极点系统的格形结构。在这里，每一级格单元不仅包含前向、后向误差的运算，还要同时对分子多项式 $B(z)$ 和分母多项式 $A(z)$ 的系数进行格形分解。

推导原理

双格形（Two-Lattice）或并行格形：可以把零极点系统看作全零点系统与全极点系统的组合，或者将其分解为两个子系统，再在格形级联中予以综合。
反射系数的含义更复杂：此时每一级的系数既要反映分子部分（零点）的贡献，也要反映分母部分（极点）的贡献，需要有相应的推导公式来确定这些格系数。

特点与分析

灵活性：可在同一结构中实现对零点与极点的联合控制；
性能更高：对特定滤波指标（如较窄带滤波、高衰减带要求等），同时调节零点与极点有更高的设计自由度；
实现复杂度略高：相较于全零点或全极点结构，多了额外的参数耦合，需要更仔细的数值分析。

3.4 格形滤波器的优点与应用

数值稳定性
- 在固定点或有限字长运算环境中，直接型结构容易出现系数量化和极点漂移问题，而格形结构由于模块化、分布式反馈，对量化误差不敏感，更易保持系统稳定。
适合自适应算法
- 传统自适应滤波器如果直接更新所有系数，容易在高阶情形下出现收敛慢或数值不稳定；格形结构可将滤波器分解成多个低阶模块，逐级更新格系数，更容易保证收敛性和稳定性。
模块化设计
- 每一级格单元实现相对简单，滤波器阶数的增加只需要在后端再级联一个格单元，具有很好的扩展性。
易于实现线性预测和语音编码
- 在语音信号处理中，线性预测编码（LPC）模型常使用格形结构来实现预测滤波器；这也是格形结构在数字信号处理领域广受青睐的一个重要原因。
实现多种滤波器类型
- 无论是全零点（FIR）、全极点（IIR）还是零极点混合（IIR）滤波器，都可以用格形结构来实现，满足不同频率特性需求。

二、特殊的数字滤波器

《见白纸笔记10》

1、微分器

2、积分器

3、匹配滤波

三、数字滤波器设计的基本思想

《见白纸笔记：12、数字滤波器设计入门》

四、典型模拟低通滤波器

《参考白纸笔记：13、模拟滤波器途径设计数字滤波器的前半部分》

1、巴特沃斯模拟低通滤波器

（1）巴特沃斯滤波器的幅频特性

（2）巴特沃斯滤波器的传递函数和极点

2、切比雪夫Ⅰ型和Ⅱ型模拟低通滤波器

（1）切比雪夫滤波器的幅频特性

（2）切比雪夫滤波器的传递函数和极点

3、椭圆形模拟低通滤波器

4、模拟滤波器设计的matlab函数

在这里讨论模拟滤波器并不是要构成一个模拟滤波器，而是在模拟滤波器的基础上通过某种变换，转换成数字滤波器(将在以后几节中讨论)，使相应的数字滤波器也有模拟滤波器的某些特性。

在MATLAB中常用的模拟滤波器设计函数如下：

（1）模拟滤波器阶数选择函数

下面结合四个典型函数（buttord、cheb1ord、cheb2ord、ellipord），详细说明模拟滤波器阶数选择函数的作用、调用格式以及各参数含义。它们都是 MATLAB 中用于根据滤波器指标自动计算模拟滤波器所需的最小阶数 n 以及相应的截止频率（或转折频率）Wn 的函数。需要注意的是，这些函数在调用时通常通过添加 's' 参数来表明设计的是“模拟（s 域）”滤波器，而非离散（z 域）滤波器。

1.1 `buttord`（巴特沃斯滤波器阶数选择函数）

函数名： buttord
功能： 计算满足指定幅度指标的巴特沃斯模拟滤波器的最小阶数以及截止频率。
调用格式：

其中：
- Wp：通带边缘频率（analog 设计时单位为 rad/s 或某些情况下 Hz，但在 MATLAB 里通常以 rad/s 为默认，具体与其他函数调用上下文有关）。
- Ws：阻带边缘频率（同上，单位与 Wp 一致）。
- Rp：通带内最大衰减（以 dB 为单位），或称通带纹波（Ripple）。
- Rs：阻带最小衰减（以 dB 为单位），或称阻带纹波。
- 's'：表明此滤波器设计在模拟域（s 域）。
- 返回值 n 是满足给定指标的最小滤波器阶数；Wn 是相应的截止频率（有时也称为 3 dB 频率或转折频率），后续可用于 butter 函数设计滤波器系数。
特性：
巴特沃斯滤波器在通带和阻带都呈单调特性，没有通带波纹或阻带波纹；相对于其他滤波器类型，它在过渡带会更“缓和”，但往往需要更高的阶数才能达到同样的衰减指标。

1.2 `cheb1ord`（切比雪夫 I 型滤波器阶数选择函数）

函数名： cheb1ord
功能： 计算满足指定幅度指标的切比雪夫 I 型（Chebyshev Type I）模拟滤波器的最小阶数以及截止频率。
调用格式：

参数含义与 buttord 类似：
- Wp、Ws、Rp、Rs、's' 与上面相同含义。
- 返回值 n、Wn 的意义也相同：滤波器阶数和相应的截止频率。
特性：
- 切比雪夫 I 型滤波器在通带呈等波纹（equiripple）特性，阻带单调下降。
- 通带纹波大小由 Rp 控制；阻带只需保证在 Ws 开始后衰减到 Rs dB 即可。
- 与巴特沃斯相比，在相同通带/阻带要求下，切比雪夫 I 型往往能用更低的阶数实现相同的幅度衰减指标，但缺点是通带会出现波纹（不再是平坦响应）。

1.3 `cheb2ord`（切比雪夫 II 型滤波器阶数选择函数）

函数名： cheb2ord
功能： 计算满足指定幅度指标的切比雪夫 II 型（Chebyshev Type II）模拟滤波器的最小阶数以及截止频率。
调用格式：
- 参数含义：与前面类似，Wp、Ws、Rp、Rs、's' 的意义不变。
- 返回值 n、Wn 同样表示滤波器阶数与截止频率。
特性：
- 切比雪夫 II 型滤波器（又称逆切比雪夫滤波器）在阻带呈等波纹特性，通带单调。
- 通带纹波可以做得很小甚至接近 0（由 Rp 控制），但阻带会有波纹。
- 在某些应用场合，对通带平坦度要求较高时，会选用切比雪夫 II 型。

1.4 `ellipord`（椭圆型滤波器阶数选择函数）

函数名： ellipord
功能： 计算满足指定幅度指标的椭圆（Cauer 或 elliptic）滤波器的最小阶数以及截止频率。
调用格式：
- 参数意义：与上面相同。
- 返回值：n、Wn 依然是滤波器阶数与截止频率。
特性：
- 椭圆滤波器在通带和阻带均呈等波纹（equiripple）。
- 在给定相同 Rp 和 Rs 指标下，椭圆滤波器通常所需的阶数是最小的，因此可以用最少的级数实现最严格的幅度衰减要求。
- 代价是响应里既有通带波纹，也有阻带波纹；相位特性也更为复杂。

1.5 参数说明与设计流程

通带边缘频率 Wp 和 阻带边缘频率 Ws
- 设计模拟滤波器时，一般以 rad/s 计；若有具体需求，也可转换为 Hz（但 MATLAB 里若使用 's'，多数场景下直接用 rad/s）。
- 若是低通滤波器，则 Wp 是低通的通带截止，Ws 是低通的阻带起始；高通、带通、带阻滤波器的频率含义则相应变化。
通带允许衰减 Rp（通带纹波）
- 以 dB 为单位，指通带最大允许的衰减量（或最大纹波深度）。
- 如 Rp=1 dB，表示通带范围内信号的衰减不会超过 1 dB。
阻带最小衰减 Rs
- 以 dB 为单位，指在阻带中至少需要达到的衰减水平。
- 如 Rs=40 dB，表示阻带范围内信号至少要衰减到 -40 dB 以下。
输出 n
- 这四个函数的第一个输出都是滤波器阶数，即最小的满足设计要求的整数阶数。
输出 Wn
- 为根据滤波器类型和指标计算得到的实际使用截止频率或转折频率。
- 该频率可直接带入对应的滤波器设计函数（如 butter, cheby1, cheby2, ellip 等）去求滤波器系数。
典型设计流程
- 第一步：根据规格（Wp、Ws、Rp、Rs）调用上述阶数选择函数，得到滤波器阶数 n 及截止频率 Wn。
- 第二步：使用与之对应的滤波器设计函数（butter, cheby1, cheby2, ellip），将得到的 n 和 Wn 代入，例如：
  $[b, a] = \text{butter}(n, Wn, 's') \quad\text{or}\quad [b, a] = \text{cheby1}(n, Rp, Wn, 's') \quad \dots$
- 第三步：得到滤波器传递函数 $H(s) = \frac{b(s)}{a(s)}$ 的系数。如果后续要数字实现，还需进行双线性变换（bilinear）或其他方法变换到 z 域。

（2）设计模拟滤波器的函数

主要围绕 MATLAB 中用于设计模拟滤波器的四个函数：butter, cheby1, cheby2, ellip，以及它们的常用参数和调用方式进行说明。

2.1 `butter`（巴特沃斯模拟滤波器）

函数名: butter
功能: 设计巴特沃斯模拟滤波器（巴特沃斯滤波器在通带内幅度响应尽可能平坦，无波纹）。
调用格式:
```
[b, a] = butter(n, Wn, 's') 
[b, a] = butter(n, Wn, 'ftype', 's')
```
其中：
- n 为滤波器阶数(正整数)。
- Wn 为截止频率或截止频带(模拟域中的角频率，单位通常为 rad/s)，当是一个数时，设计低通或高通；当是形如 [W1 W2] 的向量时，设计带通或带阻。
- 's' 表示在 s 域(模拟域) 下设计滤波器。
- 'ftype' 可选参数，用于指定滤波器类型，如 'high'（高通）、'bandpass'（带通）、'stop'（带阻）等。

2.2 `cheby1`（切比雪夫 I 型模拟滤波器）

函数名: cheby1
功能: 设计切比雪夫 I 型模拟滤波器（通带等波纹，通带有波纹，阻带则为单调衰减）。
调用格式:
```
[b, a] = cheby1(n, Rp, Wn, 's') 
[b, a] = cheby1(n, Rp, Wn, 'ftype', 's')
```
其中：
- n 为滤波器阶数。
- Rp 为通带允许的最大衰减(或通带波纹)，单位 dB。
- Wn 为截止频率或截止频带(模拟域角频率)。
- 's' 指明在 s 域设计滤波器。
- 'ftype' 同上，用于指定滤波器类型（高通、带通、带阻等）。

2.3 `cheby2`（切比雪夫 II 型模拟滤波器）

函数名: cheby2
功能: 设计切比雪夫 II 型模拟滤波器（阻带等波纹，阻带有波纹，通带则为单调）。
调用格式:
```
[b, a] = cheby2(n, Rs, Wn, 's') 
[b, a] = cheby2(n, Rs, Wn, 'ftype', 's')
```
其中：
- n 为滤波器阶数。
- Rs 为阻带最小衰减(或阻带波纹)，单位 dB。
- Wn 为截止频率或截止频带(模拟域角频率)。
- 's' 指明在 s 域设计滤波器。
- 'ftype' 用于指定滤波器类型。

2.4 `ellip`（椭圆模拟滤波器）

函数名: ellip
功能: 设计椭圆型(椭圆函数)模拟滤波器（通带、阻带均等波纹，综合性能最好但实现较复杂）。
调用格式:
```
[b, a] = ellip(n, Rp, Rs, Wn, 's') 
[b, a] = ellip(n, Rp, Rs, Wn, 'ftype', 's')
```
其中：
- n 为滤波器阶数。
- Rp 为通带波纹(单位 dB)。
- Rs 为阻带最小衰减(单位 dB)。
- Wn 为截止频率或截止频带(模拟域角频率)。
- 's' 表示在 s 域设计滤波器。
- 'ftype' 同上，用于指定滤波器类型。

2.5 参数与设计要点

n (滤波器阶数)
- 所有这四类模拟滤波器设计函数都需要滤波器阶数 n。阶数越高，理论上滤波器的过渡带越陡峭，但实现也越复杂。
- 一般在设计时，通过对通带/阻带指标(如最大通带衰减、最小阻带衰减等)进行计算或使用内置函数求解，可以确定所需的最低阶数 n。
Wn (截止频率 / 截止频带)
- 在 模拟滤波器 设计时，Wn 通常指 模拟角频率（rad/s），有时也会用归一化后的频率(具体取决于设计需求和 MATLAB 内部约定)。
- 若只给出一个数值 Wn，则默认设计低通或高通(需结合 ftype 判断)。
- 若给出 [W1 W2]，则设计带通或带阻。如：
  - 'bandpass' 时，通带在 [W1, W2] 之间。
  - 'stop'(或 'bandstop') 时，阻带在 [W1, W2] 之间。
Rp (通带波纹 / 通带衰减)
- 在切比雪夫 I 型、椭圆滤波器中，Rp 表示通带最大衰减(或波纹)，单位 dB。
- Rp 越小，通带越平坦，但相应地阶数可能需要提高。
Rs (阻带波纹 / 阻带衰减)
- 在切比雪夫 II 型、椭圆滤波器中，Rs 表示阻带最小衰减(或阻带波纹)，单位 dB。
- Rs 越大，阻带越好（衰减越大），但滤波器阶数往往也随之增大。
's'
- 这四个函数都可以通过末尾加 's' 参数来指定 在模拟域(s 域) 进行滤波器设计。
- 若不加 's'，默认是数字滤波器设计(在 z 域)；加了 's' 则是模拟滤波器设计。
'ftype'
- ftype 用于指定滤波器类型，如 'low'（低通，默认）、'high'（高通）、'bandpass'（带通）、'stop'（带阻）等。
- 当 Wn 为一个标量时，'high' 表示高通滤波器；当 Wn 为 [W1 W2]，'stop' 表示带阻滤波器，'bandpass' 表示带通滤波器，等等。

2.6 返回值[b,a]

设计滤波器时返回的 [b, a] 通常代表传递函数的分子和分母系数。具体来说：

b 是传递函数分子多项式的系数向量。
a 是传递函数分母多项式的系数向量。

传递函数 H(s)（或 H(z)）的标准形式为：

$H(s) = \frac{b_0 s^n + b_1 s^{n-1} + \cdots + b_n}{a_0 s^n + a_1 s^{n-1} + \cdots + a_n}$

其中，[b, a] 就分别包含了 b0,b1,…,bn 和 a0,a1,…,an 的数值。这样一来，你可以利用这些系数直接构建滤波器的传递函数，或用于后续的滤波器分析、仿真与实现。

（3）模拟低通滤波器原型设计函数

3.1 `buttap`（巴特沃斯滤波器原型）

函数：buttap
功能：巴特沃斯模拟低通滤波器原型
调用格式：
[z,p,k]=buttap(n)

3.2 `cheblap`（切比雪夫 I 型滤波器原型）

函数：cheblap
功能：切比雪夫 I 型模拟低通滤波器原型
调用格式：
[z,p,k]=cheblap(n,Rp)

3.3 `cheb2ap`（切比雪夫 II 型滤波器原型）

函数：cheb2ap
功能：切比雪夫 II 型模拟低通滤波器原型
调用格式：
[z,p,k]=cheb2ap(n,Rs)

3.4 `ellipap`（椭圆滤波器原型）

函数：ellipap
功能：椭圆模拟低通滤波器原型
调用格式：
[z,p,k]=ellipap(n,Rp,Rs)

说明：

以上 4 个函数中，n 是由模拟滤波器的阶数决定的次数；Rp、Rs 分别是通带内的波纹系数和阻带内的波纹系数，单位都是 dB。
通过以上 4 个函数可以求出原型滤波器的传递函数。一般表示为：
$H(s) = k\cdot \frac{(s - z_1)(s - z_2)\cdots(s - z_m)}{(s - p_1)(s - p_2)\cdots(s - p_n)}$
这样得到的滤波器原型具有对应的零点和极点，可进一步用于设计相应的模拟低通滤波器。
可以通过zp2tf()函数将零点极点增益形式转换为传递函数形式，格式为 [b,a]=zp2tf(z,p,k);

（4）频率变换函数

4.1 `lp2bp`（低通到带通变换）

函数：lp2bp
功能：低通到带通模拟滤波器变换

调用格式：

[bt, at] = lp2bp(b, a, Wo, Bw) 
[At, Bt, Ct, Dt] = lp2bp(A, B, C, D, Wo, Bw)

4.2 `lp2hp`（低通到高通变换）

函数：lp2hp
功能：低通到高通模拟滤波器变换

调用格式：

[bt, at] = lp2hp(b, a, Wo) 
[At, Bt, Ct, Dt] = lp2hp(A, B, C, D, Wo)

4.3 `lp2bs`（低通到带阻变换）

函数：lp2bs
功能：低通到带阻模拟滤波器变换

调用格式：

[bt, at] = lp2bs(b, a, Wo, Bw) 
[At, Bt, Ct, Dt] = lp2bs(A, B, C, D, Wo, Bw)

4.4 `lp2lp`（低通到低通变换）

函数：lp2lp
功能：低通到低通模拟滤波器变换

调用格式：

[bt, at] = lp2lp(b, a, Wo) 
[At, Bt, Ct, Dt] = lp2lp(A, B, C, D, Wo)

说明：

以上各个函数将截止频率为 1 rad/s 的模拟低通滤波器原型变换成带通、高通或带阻滤波器，也可以再次变换为新的低通滤波器。
Wo 表示带通滤波器的中心角频率，Bw 表示带通滤波器的带宽；对于带阻滤波器，同样用 Wo 表示中心频率，Bw 表示带宽。
若假设带通滤波器的低端频率为 W1、高端频率为 W2，则可令
$W_0 = \sqrt{W_1 \cdot W_2}, \quad B_w = W_2 - W_1$
调用中各参数含义如下：
- b, a：滤波器传递函数分子、分母多项式系数；
- A, B, C, D：滤波器的状态空间表示；
- Wo：变换后滤波器的中心(或截止)角频率；
- Bw：带宽(带通或带阻滤波器时需要)。

模拟低通滤波器可表示为

$H(s) = \frac{b_0 + b_1 s + b_2 s^2 + \cdots}{a_0 + a_1 s + a_2 s^2 + \cdots}$

经过频率变换后，可得到新的滤波器函数(或对应的状态空间方程)。若以连续状态方程表示低通滤波器，则变换后滤波器也可表示为：

$\dot{x} = A x + B u$

$y = C x + D u$

（5）`freqs`（滤波器频率分析函数）

函数：freqs
功能：如果给定模拟滤波器的分子、分母系数 b,a，即可求出该滤波器在给定频率点处的幅值和相角响应。

调用格式：

[H, w] = freqs(b, a) 
[H, w] = freqs(b, a, n) 
H = freqs(b, a, w)

说明：

在调用格式 [H, w] = freqs(b, a, n) 中，freqs 函数会在 0 到某个默认最大频率范围内选取 n 个等间隔的模拟频率点（若不指定，则默认取 200 个频率点），并计算对应的幅频和相频响应。
若调用 [H, w] = freqs(b, a)，则使用默认的频率采样点；若调用 H = freqs(b, a, w)，则可自定义一个频率向量 w（单位通常是 rad/s），freqs 会在所给定的频率向量上计算滤波器的响应。

（6）`residue`（模拟系统留数计算函数）

下面给出对图中“模拟系统留数计算函数”的详细说明与分析，包含函数功能、调用方式及其在模拟系统分析中的意义

6.1 函数名称与功能

函数名称: residue
功能: 用于计算模拟系统传递函数（或一般的有理函数）在部分分式展开（Partial Fraction Expansion）下的留数 (Residues) 和极点 (Poles)，以及可能存在的多项式项 (Direct Term)。

换言之，给定一个模拟系统的传递函数

$H(s) = \frac{b_0 + b_1 s + b_2 s^2 + \cdots + b_M s^M}{a_0 + a_1 s + a_2 s^2 + \cdots + a_N s^N}$

我们希望将其写成如下形式的部分分式之和：

$H(s) = \sum_{i=1}^{N} \frac{R_i}{s - p_i} + K(s)$

其中， $R_i$ 为对应极点 $p_i$ 的留数， $K(s)$ 为当分子次数 $\text{M} \ge \text{N}$ 时所产生的多项式项（又称多项式余项或直接项），即将 $H(s)$ 做多项式除法后剩余的那部分。

在 MATLAB 中，residue 函数就能完成这一部分分式展开，并返回每个极点对应的留数，以及多项式项的系数。

6.2 调用格式

通常，residue 的调用方式是：

$[r, p, k] = \text{residue}(B, A)$

其中：

B、A：分别为传递函数分子和分母多项式系数向量。
- 例如，若
  $H(s) = \frac{b_0 + b_1 s + b_2 s^2 + \cdots + b_M s^M}{a_0 + a_1 s + a_2 s^2 + \cdots + a_N s^N}$
  则 $B = [b_0, b_1, b_2, ..., b_M]$ ， $A = [a_0, a_1, a_2, ..., a_N]$ 。
r：函数执行后返回的留数向量 $[R_1, R_2, \ldots, R_k]$ 。
p：对应的极点向量 $[p_1, p_2, \ldots, p_k]$ 。
k：若多项式阶数 $\mathrm{M} \ge \mathrm{N}$ ，则 k 返回多项式项 $K(s)$ 的系数；否则为一个空向量 []。当分子多项式次数小于分母多项式次数时， $H(s)$ 本身就是一个严格的真有理函数，不会出现多项式项。

注意：MATLAB 中 residue 函数是既可以用来处理连续系统(s 域)的有理函数，也可以处理离散系统(z 域)的有理函数，调用方式和返回值含义相同，只是变量符号（s 或 z）不同。在模拟系统分析中，默认指的是 s 域表达式。

6.3 传递函数形式与部分分式展开

传递函数的一般形式

一个模拟系统常用的传递函数形式为：

$H(s) = \frac{b_0 + b_1 s + b_2 s^2 + \cdots + b_M s^M}{a_0 + a_1 s + a_2 s^2 + \cdots + a_N s^N}$

若 $\mathrm{M} < \mathrm{N}$ ，则是一个真有理函数；若 $\mathrm{M} \ge \mathrm{N}$ ，则可将其做多项式长除以分母多项式，得到

$H(s) = Q(s) + \frac{\text{others}}{a_0 + a_1 s + \cdots + a_N s^N}$

其中 $Q(s)$ 为多项式（即 residue 中的 k 部分）。

部分分式展开的意义

将 $H(s)$ 做部分分式展开后，可以将传递函数表示成各个极点（简单极点或高重数极点）对应的分式之和：

$H(s) = \sum_{i=1}^{k} \frac{R_i}{(s - p_i)^m} + \dots + K(s)$

对于简单极点而言，通常是

$H(s) = \sum_{i=1}^{k} \frac{R_i}{s - p_i} + K(s)$

$R_i$ ：极点 $p_i$ 对应的留数（系数）。
$p_i$ ：系统的极点，满足分母为零的值，即 $a_0 + a_1 p_i + \cdots + a_N p_i^N = 0$ 。
$K(s)$ ：可能的多项式余项。

在控制系统或信号处理分析中，留数和极点信息常用于：

逆拉普拉斯变换：快速得到系统响应的时域表达式。
网络函数分解：有助于理解电路或系统的部分响应。
稳定性分析：通过查看极点位置判断系统稳定性。
实现结构：根据分式展开，可以设计或实现由若干标准部分组成的滤波器或控制器。

6.4 使用示例（简要）

假设有一个 3 阶的模拟系统传递函数：

$H(s) = \frac{s^2 + 3s + 2}{s^3 + 2s^2 + 5s + 2}$

可在 MATLAB 中编写：

B = [1 3 2]; % 分子系数 b0 + b1*s + b2*s^2 
A = [1 2 5 2]; % 分母系数 a0 + a1*s + a2*s^2 + a3*s^3 

[r, p, k] = residue(B, A);

r 返回每个极点对应的留数；
p 返回系统极点；
k 若分子阶数 < 分母阶数，可能为空；若分子阶数 ≥ 分母阶数，则 k 会是多项式系数。

在得到 r、p、k 后，可以根据

$H(s) = \sum_{i} \frac{r_i}{s - p_i} +$ 多项式项

来重建系统，也可以用于做逆拉普拉斯变换或其他分析。

6.5 主要用途与注意事项

系统分解与逆变换
- residue 常用来做部分分式展开，是执行逆拉普拉斯变换(或逆 Z 变换)时的重要步骤。
- 对于模拟系统，若想得到时域响应 h(t) ，可以先做部分分式展开，然后对每项做逆拉普拉斯变换。
极点与稳定性
- p 是系统的极点。若所有极点实部都为负，则系统是稳定的（连续系统判断标准）。若有极点实部为正，则系统不稳定；若有极点实部为零，则需进一步分析临界稳定性。
多项式项 (k)
- 当分子多项式次数 M≥N 时，residue 会返回非空的 k，即多项式余项。可以把它看作系统的前馈项或多项式项。
- 若要做逆拉普拉斯变换，需要先把多项式项转成 $K_0 + K_1 s + \dots$ 这样的形式，再分别做变换。
数值精度
- 在高阶系统或系数数量级差异很大时，数值计算可能导致极点计算不准确，需要结合 vpa、sym 等符号工具或进行适当的数值处理。
注意变量域
- residue 本质上不区分 s 域还是 z 域，主要看你输入的多项式系数。如果是模拟系统(连续域)，则 B 和 A 对应 s 的多项式；若是数字系统(离散域)，则对应 z 的多项式。

6.6 小结

residue 是 MATLAB 中常用的一个函数，用于部分分式展开和获取系统极点与留数。在模拟滤波器或模拟控制系统的分析中，通过 residue(B, A) 可以轻松获得系统的零极点信息（尤其是极点和对应留数），并可进一步进行系统响应分析、稳定性分析或逆拉普拉斯变换等操作。

5、案例

《见案例程序，篇幅原因这里不过多分析》

案例3.1、巴特沃斯、切比雪夫Ⅰ型、切比雪夫Ⅱ型和椭圆型滤波器的相同和不同之处

案例3.2、设计模拟滤波器的几种编程方法的相同和不同之处

案例 3.3、在频带变换的模拟滤波器设计中，怎样计算Wn和Bs

五、利用模拟滤波器设计数字滤波器

《参考白纸笔记：13、模拟滤波器途径设计数字滤波器的后半部分》

利用模拟滤波器设计数字滤波器的方法主要基于以下思路：先设计出一个满足要求的模拟滤波器（原型滤波器），再通过一定的映射方法将其转换为数字滤波器。

1、基本原理

模拟滤波器设计技术成熟、设计公式齐全，可以根据频率响应指标（如截止频率、通带/阻带波纹等）设计出所需的滤波器原型。将模拟滤波器转换为数字滤波器的核心问题在于如何将连续时间系统（s 域）的特性映射到离散时间系统（z 域）。这个过程必须保证映射后的系统满足如下要求：

稳定性：模拟滤波器所有极点位于左半平面（Re(s)<0），转换后对应的数字滤波器极点应位于单位圆内（|z|<1）。
频率响应：尽可能保持模拟滤波器设计时所期望的幅频特性或相位特性。
时间响应：在一些方法中，还要求数字滤波器的脉冲响应与模拟滤波器的脉冲响应尽可能一致（在采样时刻）。

2、常见的映射方法：冲击不变法（脉冲不变法）

基本原理：
直接对模拟滤波器的脉冲响应 $h_a(t)$ 进行采样，令

$h[n] = T \cdot h_a(nT)$

再通过 Z 变换得到数字滤波器的传递函数。

主要步骤：

求模拟滤波器脉冲响应：对设计好的模拟传递函数 $H(s)$ 求逆拉普拉斯变换，得到 $h_a(t)$ 。
采样：按照采样周期 T 采样，得到离散序列 $h[n]$ 。
Z 变换：对 $h[n]$ 做 Z 变换，得到数字滤波器的传递函数 $H(z)$ 。

优缺点：

优点：能保持模拟滤波器在采样点上的时域响应；适用于脉冲响应主要集中在低频部分的情况。
缺点：频谱搬移造成频谱泄露；难以确定 $h_a(t)$ 从无限长到有限长的截断。因此，这种方法不好用，实际用的少。

3、常见的映射方法：双线性变换法

基本原理：
直接在频域从模拟到数字的映射，采用双线性变换公式：

$s = \frac{2}{T} \cdot \frac{z-1}{z+1}$

其中 T 为采样周期。此变换将 s 域的左半平面映射到 z 域的单位圆内，保证稳定性，同时将无穷远点映射到 z = -1。

$H(z)=H_{a}(\frac{2}{T}\cdot \frac{z-1}{z+1})$

因此，模拟频率 ω 与数字频率 Ω 之间的对应关系为：

$\boxed{\omega=\frac{2}{T}\tan\left(\frac{\Omega}{2}\right)}$

或者反过来写数字频率与模拟频率的关系：

$\boxed{\Omega=2\arctan\left(\frac{\omega T}{2}\right)}$

主要步骤：

设计模拟原型：例如使用 butter, cheby1, ellip 等函数设计出模拟低通滤波器原型。
预畸变频率：由于双线性变换引起频率非线性变换（频率扭曲），需要进行预畸变（Pre-warping），即计算模拟滤波器的截止频率 $\omega_c$ 对应数字滤波器的截止频率 $\Omega_c$ 的关系：
$\omega_c = \frac{2}{T}\tan\left(\frac{\Omega_c T}{2}\right)$
变换公式代入：将模拟传递函数中的 s 用双线性变换公式替换，得到数字传递函数 $H(z)$ 。

优缺点：

优点：稳定性有保证；对于低通、高通、带通、带阻滤波器都有广泛应用；可利用 MATLAB 中现成函数（如 bilinear）实现。
缺点：频率响应存在非线性畸变（频率扭曲），需要预畸变来补偿，且带外响应会有一定变化。

3.5、双线性映射方法设计数字滤波器步骤

（1）目标：数字滤波器的设计目标。

（2）根据模拟频率 ω 与数字频率 Ω 之间的对应关系将设计目标转为模拟的

（3）找到好的模拟滤波器 $H_a(s)$ （几种经典模型或更多其他好的模型）

（4）带入 $s = \frac{2}{T} \cdot \frac{z-1}{z+1}$ 到 $H_a(s)$ 得到： $H_d(z)$

（5）最后由 $H_d(z)$ 得到 $h_d(k)$

4、常见的映射方法：匹配 z 变换法（Matched z-Transform / Pole-Zero Mapping）

基本原理：
直接将模拟滤波器的极点和零点通过映射关系转换到 z 域。例如，对于每个模拟极点 $p_i$ 和零点 $z_i$ ，采用映射公式：

$z = e^{p_i T}$

从而得到数字滤波器的极点和零点。

主要步骤：

极点零点分解：将模拟传递函数 $H(s)$ 表示为极点和零点形式。
极点映射：对每个模拟极点 $p_i$ 使用 $z_i = e^{p_i T}$ 得到对应的数字极点；同理，对于零点进行相同的映射（注意：有时模拟系统可能没有零点，数字系统也需要适当调整）。
增益调整：计算数字滤波器的增益，使得在某些参考频率处数字滤波器的响应与模拟滤波器一致。

优缺点：

优点：能直接反映模拟滤波器的极点分布，通常在低频范围内能较好地匹配时域特性；对某些应用（如信号整形）较有优势。
缺点：不保证整体频率响应匹配，可能在某些频段引入失真；零点映射处理较为复杂，需保证数字滤波器的零点结构合理。

六、陷波器与全通滤波器

1、陷波器

在信号测量中，有时信号会被一些干扰的正弦信号所淹没，所以要通过一些窄带滤波器把这些正弦信号滤除。陷波器的目的就是构成窄带陷波，用来滤除这些单个的正弦干扰信号，而其他频率成分都能通过。

这里讨论二阶数字陷波器，假设要滤除的频率为 $\omega_0$ ，所以滤波器的幅度特性在 $\omega = \omega_0$ 处为零，在其他频率上接近常数，这就构成了滤除单频干扰的滤波器。

令零点 $z = e^{\pm j\omega_0}$ 使滤波器的幅度特性在 $\omega = \pm \omega_0$ 处为零。为了使幅度离开 $\omega = \pm \omega_0$ 后迅速上升到一个常数，将两个极点放在非常靠近零点的地方，极点为 $z = re^{\pm j\omega_0}$ ，传递函数为：

$H(z) = \frac{b_0 (z - e^{j\omega_0})(z - e^{-j\omega_0})}{(z - re^{j\omega_0})(z - re^{-j\omega_0})} = \frac{b_0 [1 - 2\cos(\omega_0)z^{-1} + z^{-2}]}{1 - 2r\cos(\omega_0)z^{-1} + r^2z^{-2}} \quad (3-5-1)$

其中 $0 < r \leq 1$ ，极点的位置是 $z = re^{\pm j\omega_0}$ ，半径 $r$ 与陷波器的带宽有关。陷波器在 -3 dB 处的频带定义为带宽 $Bw$ ，则 $Bw$ 和增益系数 $b_0$ 为：

$Bw = (1 - r)f_s, \quad b_0 = \frac{\left | 1 - 2r\cos\omega_0 + r^2 \right |}{2(1 - \cos\omega_0)} \quad (3-5-2)$

图中展示了 Bw = 0.01, 0.05, 0.1 三种情况下陷波器的幅度响应。

分析：

若 r 较小，滤波器的阻带宽 Bw 较宽， $\omega = \pm \omega_0$ 附近的频率成分会有显著影响。
当 $r \to 1$ 时，陷波器带宽 Bw 减小，同时陷波深度也减小。

因此，设计陷波器时须根据实际需求选择合适的 r 值。

2、全通滤波器

（1）概念

全通滤波器是一种幅度响应恒定为1的滤波器（即对所有频率都不衰减），但具有可调的相位响应。

其幅度特性：

$|H(e^{j\omega})| = 1, \quad 0 \leq \omega \leq 2\pi$

全通滤波器的功能不在于“放大或抑制”信号的幅度，而是“调整相位”。
关键用途：相位补偿、群延迟控制、相位均衡。

（2）基本结构与推导

$H_{ap}(z) = \frac{\sum\limits_{i=0}^{N} a_i z^{-N+i}}{\sum\limits_{i=0}^{N} a_i z^{-i}}, \quad a_0 = 1$

分子与分母系数是对称的反转：称为“反射对称结构”。
若定义分母为 $A(z)$ ，则有 $H_{ap}(z) = z^{-N} \frac{A(z^{-1})}{A(z)}$ 。

（3）极点零点对称关系

若有极点 $p_i$ ，则必须有对应的零点 $z_i = \frac{1}{p_i^*}$ 。

$H_{ap}(z) = \prod_{i=1}^N \frac{z^{-1} - p_i^*}{1 - p_i z^{-1}}$

若滤波器是实系数，则共轭对称。
若要求滤波器稳定（极点在单位圆内），则所有极点 $|p_i| < 1$ 。

（4）推广形式

包括实数极点（ $\alpha_i$ ）和复数极点（ $\beta_i$ ）两部分：

$H_{ap}(z) = \prod_{i=1}^{K_R} \frac{z^{-1} - \alpha_i}{1 - \alpha_i z^{-1}} \cdot \prod_{i=1}^{K_C} \frac{(z^{-1} - \beta_i^*)(z^{-1} - \beta_i)}{(1 - \beta_i^* z^{-1})(1 - \beta_i z^{-1})}$

$\alpha_i$ ：实极点
$\beta_i$ ：复极点
$K_R$ 、 $K_C$ ：分别为实极点数和复极点对数

（5）极点零点分布图分析

(a) 实数极点零点：

零点位于极点的倒数处，即 $z = \alpha^{-1}$ 。

(b) 复数极点零点：

极点在单位圆内，如 $\beta = r e^{j\theta}$ ，对应零点 $1/\beta^* = r^{-1} e^{j\theta}$ 在单位圆外。
极点与零点关于单位圆共轭对称。

（6）相位响应与群延迟

对于只有一个极点的情况：

$H_{ap}(e^{j\omega}) = e^{-j\omega} \cdot \frac{1 - r e^{j(\omega - \theta)}}{1 - r e^{-j(\omega - \theta)}}$

相位响应：

$\varphi_{ap}(\omega) = -\omega - 2 \arctan \left( \frac{r \sin(\omega - \theta)}{1 - r \cos(\omega - \theta)} \right)$

群延迟：

$\tau_{\varphi}(\omega) = - \frac{d\varphi_{ap}(\omega)}{d\omega} = \frac{1 - r^2}{1 + r^2 - 2r \cos(\omega - \theta)} \quad > 0$

群延迟总为正值，表示信号被延迟。
在 IIR 滤波器中常用来补偿相位非线性。

七、IIR数字滤波器设计的matlab函数

在IIR数字滤波器设计中有 <1>直接的设计函数，即在已知数字滤波器的指标后调用函数直接设计得到滤波器的系数，或者 <2>先设计模拟滤波器，再通过转换得到数字滤波器的系数。以下分别给出它们的MATLAB函数。不论用哪一种方法设计的数字滤波器都以巴特沃斯、切比雪夫工型、切比雪夫Ⅱ型和椭圆型模拟滤波器为原型，并把它们映射为数字滤波器。

1、IIR数字滤波器的阶数选择函数

IIR数字滤波器设计中与阶数选择有关的MATLAB函数有以下几种：

（1）buttord（计算巴特沃斯数字滤波器的阶数）

其调用格式为：

[n, Wn] = buttord(Wp, Ws, Rp, Rs)

其中，Wp 是通带截止频率，Ws 是阻带截止频率，Rp 是通带最大衰减（dB），Rs 是阻带最小衰减（dB），n 是计算出的滤波器阶数，Wn 是滤波器的归一化截止频率。

（2）cheb1ord（计算切比雪夫I型数字滤波器的阶数）

其调用格式为：

[n, Wn] = cheb1ord(Wp, Ws, Rp, Rs)

与 buttord 相似，Wp、Ws、Rp 和 Rs 分别代表通带和阻带的截止频率以及对应的最大衰减。

（3）cheb2ord（计算切比雪夫II型数字滤波器的阶数）

其调用格式为：

[n, Wn] = cheb2ord(Wp, Ws, Rp, Rs)

与 cheb1ord 类似，但适用于切比雪夫II型滤波器，参数含义相同。

（4）ellipord（计算椭圆型数字滤波器的阶数）

其调用格式为：

[n, Wn] = ellipord(Wp, Ws, Rp, Rs)

该函数用于设计椭圆型数字滤波器，Wp、Ws、Rp 和 Rs 分别代表通带和阻带的截止频率及其对应的衰减。

这些函数的共同点是都需要输入滤波器的通带和阻带的截止频率，以及对应的衰减参数，然后它们会根据输入的要求计算并返回合适的滤波器阶数和归一化截止频率。

以上4个函数与在模拟滤波器设计中一样,只是在输入参数中不再带有's'。输出参数中的 n 是求出滤波器最小的阶数，Wn是等效低通滤波器的截止频率；Wp和Ws分别是通带和阻带的频率(截止频率)，是归一化的角频率，单位为rad/pi，数值在0~1之间。当Wp>Ws时，为高通滤波器；当Wp和Ws为二元矢量时，为带通或带阻滤波器，这时求出的Wn也为二元矢量。

2、模拟滤波器数字化

（1）impinvar（脉冲响应不变法）

impinvar：用于模拟滤波器的脉冲响应不变法（impulse invariance method）将模拟滤波器转换为数字滤波器。其调用格式为：

[bz, az] = impinvar(b, a, Fs)

其中，b 和 a 是模拟滤波器的系统系数，bz 和 az 是数字滤波器的系数，Fs 是采样频率。该方法通过保持模拟滤波器的脉冲响应与数字滤波器相同，进行数字化处理。

说明：b 和 a 是模拟滤波器的系统系数，bz 和 az 是数字滤波器的系统系数。Fs 是采样频率，它定义了数字滤波器的采样频率，通常与模拟滤波器的脉冲响应频率保持一致。使用 impinvar(b, a) 时，采样频率默认设置为 1 Hz。

（2）bilinear（双线性映射法）

bilinear：用于双线性 Z 变换方法将模拟滤波器转换为数字滤波器。其调用格式为：

[bz, az] = bilinear(b, a, Fs) 
[zd, pd, kd] = bilinear(z, p, k, Fs)

其中，b 和 a 是模拟滤波器的系统系数，bz 和 az 是数字滤波器的系数，Fs 是采样频率。第二种形式适用于双线性 Z 变换的极点零点转换，z、p 和 k 是模拟滤波器的零点、极点和增益，zd、pd 和 kd 是对应的数字滤波器的零点、极点和增益。

说明：该方法通过双线性 Z 变换将模拟滤波器的零点和极点映射到数字域，转换后数字滤波器的系数由 bz 和 az 给出。Fs 是采样频率，定义了数字滤波器的采样率。

这两个函数常用于将模拟滤波器（如巴特沃斯、切比雪夫滤波器等）转换为数字滤波器，在实际设计中非常有用。

3、直接设计数字滤波器的函数

IIR数字滤波器设计中与直接设计数字滤波器有关的 MATLAB 函数如下所示。这些函数均是根据设定的滤波器阶数和性能指标（如通带/阻带衰减）直接生成数字滤波器系数：

（1）`butter`（巴特沃斯滤波器设计函数）

功能：设计巴特沃斯数字滤波器。

调用格式：

[b, a] = butter(n, Wn) 
[b, a] = butter(n, Wn, 'ftype')

n：滤波器阶数
Wn：归一化截止频率
'ftype' 可选值：'high'（高通），'stop'（带阻）

（2）`cheby1`（切比雪夫Ⅰ型滤波器设计函数）

功能：设计通带等波纹（等波纹在通带内）的切比雪夫Ⅰ型数字滤波器。

调用格式：

[b, a] = cheby1(n, Rp, Wn) 
[b, a] = cheby1(n, Rp, Wn, 'ftype')

Rp：通带最大衰减（dB）

（3）`cheby2`（切比雪夫Ⅱ型滤波器设计函数）

功能：设计阻带等波纹的切比雪夫Ⅱ型数字滤波器。

调用格式：

[b, a] = cheby2(n, Rs, Wn) 
[b, a] = cheby2(n, Rs, Wn, 'ftype')

Rs：阻带最小衰减（dB）

（4）`ellip`（椭圆型滤波器设计函数）

功能：设计具有通带与阻带同时等波纹的椭圆型数字滤波器（最小阶数实现指定指标）。

调用格式：

[b, a] = ellip(n, Rp, Rs, Wn) 
[b, a] = ellip(n, Rp, Rs, Wn, 'ftype')

通用说明

这些函数的参数 Wn 可以是一个标量（低通/高通滤波器）或一个二元素向量 [W1, W2]（用于带通/带阻滤波器）。
如果使用 'ftype' 参数：
- 'high'：设计高通滤波器，截止频率为 Wn
- 'stop'：设计带阻滤波器，截止频率为 [W1, W2]
所有频率参数 Wn 为归一化频率，其取值范围为 0 ~ 1（对应于 0 ~ π radians/sample）。

4、滤波器特性分析函数

（1）`freqz`（求得数字滤波器的幅值响应和相位响应）

功能：已知数字滤波器的系统系数 b 和 a，求得数字滤波器的幅值响应和相位响应。

调用格式：

[H, w] = freqz(b, a) 
[H, w] = freqz(b, a, n) 
[H, f] = freqz(b, a, n, Fs) 
[H, w] = freqz(b, a, n, 'whole') 
[H, f] = freqz(b, a, n, 'whole', Fs) 
H = freqs(b, a, w) 
freqz(b, a)

说明：
在调用格式[H,w]=freqz(b,a)中，freqz函数自定义计算512个数字频率，响应曲线在H中，对应的数字角频率在矢量w中，w在 0~ $\pi$ 区间内，表示归一化角频率。参数 b 和 a 是数字滤波器的系数。设置参数 n 后，即表示频率矢量和 H 都有 n 个元素；当带有采样频率Fs后，输出频率矢量f，将在0~Fs/2的区间内。设置参数 'whole' 后频率将在0~2 $\pi$ 区间内，或在 0~Fs 的区间内。当不带任何输出参数时，freqz(b,a)将直接显示出该滤波器的幅频响应曲线与相频响应曲线。

（2）`fvtool`（显示数字滤波器的响应曲线，支持多个滤波器一起显示）

功能：显示数字滤波器的响应曲线；支持多个滤波器的频率响应一起显示。

调用格式：

fvtool(b, a) 
fvtool(Hd) 
fvtool(b1, a1, b2, a2, ..., bn, an) 
fvtool(Hd1, Hd2, ...)

说明：
- 输人参数b和a是数字滤波器系数；b1,a1,b2,a2,…,bn,an是多个滤波器的系数(要求多个滤波器的响应曲线在一张图上画出)；Hd是滤波器系数集合；Hd1，Hd2，…是多个滤波器的系数集合。从调用格式可以看出，多个滤波器的响应曲线通过本函数一次就能同时画在一张图上。
- 该函数通过图形界面显示滤波器的频率响应，适用于多个滤波器的比较。

（3）`grpdelay`（群延迟分析）

功能：已知数字滤波器系数 b 和 a，输出数字滤波器群延迟随频率变化的曲线。

调用格式：

[gd, w] = grpdelay(b, a) 
[gd, w] = grpdelay(b, a, n) 
[gd, f] = grpdelay(b, a, n, Fs) 
[gd, w] = grpdelay(b, a, 'whole') 
[gd, f] = grpdelay(b, a, 'whole', Fs) 
gd = grpdelay(b, a, w) 
grpdelay(b, a)

说明：
- 群延迟表示滤波器相位响应的导数，反映了信号各频率成分的时延变化。
- grpdelay函数类似于freqz函数，自定义计算512个数字频率，设置n后就有n个数字频率。延迟量曲线在gd中，对应的数字角频率在矢量w中，w在0~π区间内，表示归一化角频率。w和gd的长度或为512个(默认值)或为n个。参数b和a是数字滤波器的系数，对应的群延迟响应为
  $\tau _{g}(\omega )=-\frac{d\varphi (\omega )}{d\omega }$
  带有采样频率Fs时，输出频率矢量f将在0~Fs/2的区间内。设置参数whole'后频率将在 $0\sim 2\pi$ 区间内，或在0~Fs区间内。不带任何输出参数时同freqz(b,a)，将直接显示出该滤波器的群延迟随频率变化的曲线。

（4）`impz`（生成数字滤波器的脉冲响应）

功能：已知数字滤波器系数 b 和 a，求出数字滤波器的脉冲响应。

调用格式：

[h, t] = impz(b, a) 
[h, t] = impz(b, a, n) 
[h, t] = impz(b, a, n, Fs)
impz(b, a, n, Fs)

说明：

b和a为数字滤波器系数，n为脉冲响应样点总数，Fs为采样频率，默认值为1，h为滤波器单位脉冲响应向量，t为[0:n一1]’，是h对应的时间向量。当函数输出默认时，绘制滤波器脉冲响应图，当n为默认时,函数自动选择n值。

（5）`zplane`（绘制系统的零极点图）

功能：已知数字滤波器的零点 z 和极点 p，或滤波器系数 b 和 a，绘制系统的零极点图。
调用格式：
```
zplane(z, p) 
zplane(b, a)
```
说明：
- z 和 p 是系统的零点和极点，b 和 a 是滤波器的系统系数。
- 该函数在复平面上绘制零极点图，用于分析系统的稳定性和响应特性。函数在Z平面绘出零点和极点。极点用x表示,零点用○表示。

5、滤波使用函数

（1）`filter`（数字滤波器的应用）

功能：实现数字滤波器对数据的滤波。

调用格式：

y = filter(b, a, x) 
[y, zf] = filter(b, a, x, zi)

说明：
- b 和 a 为数字滤波器系数；x为滤波器的离散输入，y为滤波器的输出，y为与x具有相同大小的向量，zi和 zf 表示x在分段滤波时的初始状态和最终状态

（2）`filtfilt`（零相位数字滤波）

功能：实现数字滤波器对数据进行零相位滤波。
调用格式：
```
y = filtfilt(b, a, x)
```
说明：
- b 和 a 是数字滤波器的系统系数，x 是滤波器的输入，y 是滤波器的输出。
- 与 filter 函数不同，filtfilt 函数执行 双向滤波，即向前和向后滤波，确保输出信号的相位不变（零相位失真）。

6、数字陷波器和全相位数字滤波器的相关函数

（1）`iirnotch`（数字陷波器设计 Notch Filter）

函数名称: iirnotch
功能: 设计数字陷波器。
调用格式:
```
[b, a] = iirnotch(Wo, Bw)
```
- Wo 是陷波器的中心频率。
- Bw 是陷波器的带宽。
- b 和 a 是数字滤波器的系数。

（2）`iirgrpdely`（全相位数字滤波器设计 All-pass Digital Filter）

函数名称: iirgrpdely
功能: 设计一个在已知频率范围内接近规定群延迟的全相位数字滤波器
调用格式：
```
[b, a] = iirgrpdely(N, F, Edges, Gd)
```
- N 是滤波器的阶数（必须是偶数）。
- F 和 Gd 是指定在某一频率区间内群延迟的指标。
- Edges 是指定滤波器的边缘频率。
- b 和 a 是数字滤波器的系数。

这些函数主要用于设计陷波滤波器和全相位数字滤波器，以满足特定的频率特性或群延迟要求。在MATLAB中，可以利用这些函数设计相应的IIR滤波器，调节参数以获得所需的性能。

7、`residuez` 和 `residuesz`（数字系统留数的计算函数）

数字系统留数的计算函数：

函数名称: residuez 和 residuesz
功能: 计算数字系统的留数和极点，或者计算数字系统系数。
调用格式:
```
[r, p, k] = residuez(B, A) 
[B, A] = residuesz(r, p, k)
```
- B 和 A 是数字滤波器的系数（分子和分母）。
- r 是系统的留数。
- p 是系统的极点。
- k 是系统的直接增益。

解释：

residuez 函数用于计算数字系统的留数（r）、极点（p）和直接增益（k）。
residuesz 函数则是反向操作，当已经知道数字系统的留数（r）、极点（p）和增益（k）时，能够恢复原始的数字滤波器系数 B 和 A。

8、`freqz_m`（分析滤波器的频率特性：幅，相，群延迟）

一些有用的MATLAB函数。

有一些MATLAB函数不是MATLAB自带的,而是由其他一些书本自编的,但对IIR和FIR滤波器设计和应用很有用。

这些函数已包含在案例附带程序包中的basic_tbxsp 子目录中。

函数名称: freqz_m

功能: 分析滤波器的幅度响应的频率特性、相位响应及群延迟与频率的关系。
调用格式:
```
[db, mag, pha, grd, w] = freqz_m(b, a)
```
- db 是幅度响应的分贝值。
- mag 是幅度响应的线性值。
- pha 是相位响应。
- grd 是群延迟响应。
- w 是一组统一的角频率。

说明：

freqz_m 是一个用于分析滤波器频率响应的函数，特别用于计算滤波器的幅度响应、相位响应和群延迟等特性。
db 表示滤波器的幅度响应的分贝值，mag 是其线性值，pha 表示相位响应，grd 是群延迟。
w 是频率范围，它是与 freqz 函数类似的，但在实际应用时，这个函数可以提供更多的滤波器特性分析参数。

七、IIR滤波器设计的案例

《见案例程序，篇幅原因这里不过多分析》

案例3.4、用留数求得脉冲不变法数字滤波器与调用 impinvar 函数得到的是否一样

案例3.5、在调用bilinear函数时为何有的Fs处用实际频率值，有的却用Fs=1

案例3.6、为什么不能用 impinvar 函数

案例3.7、为什么滤波器的输出会溢出或没有数值

案例3.8、用 bilinear 函数时,如果Wp和Ws都没有先做预畸会有什么结果

案例3.9、如何把任意S系统转换为Z系统

案例3.10、把滤波器的滤波过程用差分方程的运算来完成

案例3.11、滤波函数 filter 的调用格式为[y,zf]=filter(b,a,x,zi),其中的zi和zf有何作用

案例3.12、如何使用数字陷波器滤除工频信号

案例3.13、如何设计数字全通滤波器对 IIR滤波器进行相位补偿

案例3.14、为什么零相位滤波在起始和结束两端都受瞬态效应的影响

八、线性相位与FIR系统的相位特性

《参考白纸笔记：12、FIR系统线性相位的h(n)的4种形式》

九、FIR型数字滤波器的窗函数设计法

1、理想数字滤波器的单位脉冲响应

在设计FIR滤波器时，通常首先考虑理想滤波器的单位脉冲响应。理想滤波器在频域上具有完美的通带与阻带划分，但其时域单位脉冲响应往往是无限长的，不能直接实现，需要通过截断获得有限长的响应。下面分别讨论低通、高通、带通、带阻滤波器的理想单位脉冲响应。

（1）理想数字低通滤波器

理想低通滤波器的频率响应定义为：

$H_d(e^{j\omega}) = \begin{cases} 1, & |\omega| \leq \omega_c \\ 0, & \omega_c < |\omega| \leq \pi \end{cases}$

其中， $\omega_c$ 为截止频率。

通过傅里叶反变换可得其时域响应：

$h_d[n] = \frac{1}{2\pi}\int_{-\omega_c}^{\omega_c} e^{j\omega n} d\omega = \frac{\sin(\omega_c n)}{\pi n}, \quad n\neq 0$

而在 $n=0$ 时，由极限处理可得：

$h_d[0] = \frac{\omega_c}{\pi}$

因此，理想低通滤波器的单位脉冲响应为无限长的sinc函数。

（2）理想数字高通滤波器

理想高通滤波器的频率响应为：

$H_d(e^{j\omega}) = \begin{cases} 0, & |\omega| \leq \omega_c \\ 1, & \omega_c < |\omega| \leq \pi \end{cases}$

利用频谱平移的方法，可以将高通滤波器视为低通滤波器经过调制（频谱镜像）。其单位脉冲响应可以写为：

$h_d[n] = \delta[n] - \frac{\sin(\omega_c n)}{\pi n},\quad n\neq 0$

而在 $n=0$ 时：

$h_d[0] = 1 - \frac{\omega_c}{\pi}$

这样，同样是无限长的序列，但中心有一个直流脉冲，后续项为sinc函数的负值。

（3）理想数字带通滤波器

理想带通滤波器的频率响应在两个截止频率 $\omega_{c1}$ 和 $\omega_{c2}$ 之间为1，其他为0，即：

$H_d(e^{j\omega}) = \begin{cases} 1, & \omega_{c1} \leq |\omega| \leq \omega_{c2} \\ 0, & \text{others} \end{cases}$

其单位脉冲响应为两个低通滤波器响应的差分，即：

$h_d[n] = \frac{1}{\pi n} \Big[\sin(\omega_{c2} n) - \sin(\omega_{c1} n)\Big], \quad n\neq 0$

在 $n=0$ 时：

$h_d[0] = \frac{\omega_{c2}-\omega_{c1}}{\pi}$

（4）理想数字带阻滤波器

理想带阻滤波器的频率响应则在带阻区间为0，其他区间为1，可表示为：

$H_d(e^{j\omega}) = \begin{cases} 1, & |\omega| \leq \omega_{c1} \text{ or } |\omega| \geq \omega_{c2} \\ 0, & \omega_{c1} < |\omega| < \omega_{c2} \end{cases}$

利用补集的思想，其单位脉冲响应为：

$h_d[n] = \delta[n] - \frac{1}{\pi n}\Big[\sin(\omega_{c2} n)-\sin(\omega_{c1} n)\Big],\quad n\neq 0$

在 $n=0$ 时：

$h_d[0] = 1 - \frac{\omega_{c2}-\omega_{c1}}{\pi}$

总结：理想滤波器的单位脉冲响应虽然在频域满足理想化要求，但由于时域响应无限长，因此必须通过窗函数截断，使其可实现化。

2、FIR型数字滤波器的矩形窗设计法

FIR数字滤波器设计的基本思想是先得到理想滤波器的无限长脉冲响应，然后使用窗函数截断以获得有限长度的FIR滤波器。矩形窗是最简单的一种窗函数，其设计过程和缺陷是后续加窗方法改进的基础。

（1）矩形窗设计过程

矩形窗设计法的基本步骤包括以下几个部分：

1.1 矩形窗截断过程

从理想滤波器的单位脉冲响应 $h_d[n]$ （无限长）开始，通过一个长度为 N 的矩形窗 $w[n]$ 将其截断。矩形窗定义为：

$w[n] = \begin{cases} 1, & 0 \leq n \leq N-1 \\ 0, & \text{others} \end{cases}$

截断过程即计算：

$h[n] = h_d[n] \cdot w[n]$

在实际设计中，由于要求滤波器因果，常常需要对 $h_d[n]$ 进行中心化截取，再平移为因果滤波器。

1.2 因果性处理过程

由于理想单位脉冲响应 $h_d[n]$ 通常是关于 $n=0$ 对称的（或延迟 N/2 ），需要对截断后的序列进行时移，使得系统成为因果系统。令：

$h_{\text{causal}}[n] = h[n - \alpha]$

其中 $\alpha$ 通常选取为 $\alpha = \frac{N-1}{2}$ （当 N 为奇数时，滤波器关于中心对称；当 N 为偶数时，则需做适当处理）。这样，滤波器从 $n=0$ 开始，满足实际系统对因果性的要求。

1.3 计算频率响应

滤波器的频率响应 $H(e^{j\omega})$ 为截断后的序列 $h[n]$ 的离散时间傅里叶变换（DTFT）：

$H(e^{j\omega}) = \sum_{n=0}^{N-1} h[n] e^{-j\omega n}$

对于矩形窗截断的情况，实际频率响应可看作理想滤波器响应与矩形窗傅里叶变换的卷积，从而引入旁瓣效应和通带波动。

（2）矩形窗截断的频谱

利用矩形窗截断必然会导致频谱泄漏现象。矩形窗的频谱（窗函数的DTFT）为：

$W(e^{j\omega}) = \sum_{n=0}^{N-1} e^{-j\omega n} = e^{-j\omega\frac{N-1}{2}} \frac{\sin\left(\frac{N\omega}{2}\right)}{\sin\left(\frac{\omega}{2}\right)}$

从上式可以看出：

主瓣宽度：主瓣宽度大约为 $\frac{4\pi}{N}$ （具体定义中常以零交叉点距离确定）。
旁瓣衰减：矩形窗的旁瓣衰减率较低（约为20 dB左右），这会导致滤波器的阻带衰减性能较差。

矩形窗的这种特性决定了使用矩形窗设计的FIR滤波器在阻带中会有较高的旁瓣泄漏，可能无法满足严格的阻带衰减要求。

（3）窗函数指标

在窗函数设计中，有几个关键指标用来描述窗函数的性能：

主瓣宽度（Mainlobe Width）：主瓣宽度决定了滤波器的过渡带宽度。主瓣越窄，滤波器过渡带越陡峭。
旁瓣衰减（Sidelobe Attenuation）：旁瓣的最大衰减值，表示滤波器在阻带的泄漏水平。旁瓣衰减越高，滤波器在阻带的抑制效果越好。
能量分布：能量集中在主瓣内的窗函数往往具有更好的频域特性。

矩形窗虽然实现简单，但其主瓣较窄（可以实现较陡的过渡带）却牺牲了旁瓣衰减，因而在设计上存在一定折中。

3、窗函数设计法

由于矩形窗存在旁瓣泄漏较高的问题，后续发展出多种改进型窗函数。之前介绍过的任何一个窗,其参数都是固定的，都不能同时实现窗函数的几个指标的优化。因此，需要根据应用不同，选择不同的窗函数。而凯泽窗可通过调节窗函数的参数在各指标之间进行折中。下面重点讨论这一种重要的窗函数——凯泽窗（Kaiser窗），以及加窗后频谱的特性。

（1）凯泽窗函数

凯泽窗是一种参数化窗函数，其表达式为：

$w[n] = \frac{I_0\left(\beta \sqrt{1-\left(\frac{2n}{N-1}-1\right)^2}\right)}{I_0(\beta)},\quad 0\le n\le N-1$

其中 $I_0(\cdot)$ 表示零阶改进贝塞尔函数，参数 β 用来调控窗函数的旁瓣衰减和主瓣宽度：

β 越大，旁瓣衰减越高，但主瓣宽度也会增大，导致过渡带变宽；
β 较小时，窗函数更接近矩形窗，主瓣较窄，但旁瓣较高。

当 β =0时,凯泽窗相当于矩形窗;当 β =5.44时,凯泽窗相当于海明窗;当 β =8.5时,凯泽窗相当于布莱克曼窗。选择不同的值即可得到不同的频域特征。Kaiser 经过大量数值实验得到了一组设计数字滤波器凯泽窗参数估计公式,即估算 N 和 β 值。

N 的计算公式：

$N \approx \frac{A_s - 7.95}{14.36 \times \Delta f} + 2 \quad \text{(3-9-11)}$

其中，Δf 是通带宽度，计算公式为： $\Delta f = \frac{(\omega_s - \omega_p)}{2 \pi}$

β 的计算公式如下：

$\beta = \begin{cases} 0.1102(A_s - 8.7), & A_s > 50 \\ 0.5842 (A_s - 21)^{0.4} + 0.07886(A_s - 21), & 21 \leq A_s \leq 50 \\ 0, & A_s \leq 21 \end{cases}$

其中，As 阻带衰减（单位为dB）。

凯泽窗的优势在于其可以在主瓣宽度与旁瓣衰减之间实现灵活调控。利用凯泽窗设计FIR滤波器时，先计算理想滤波器的单位脉冲响应，再乘以凯泽窗，从而得到加窗后的有限长序列：

$h[n] = h_d[n] \cdot w[n]$

这种方法在许多工程应用中非常流行，因为它既能较好地控制旁瓣衰减，又能兼顾过渡带陡峭性。

（2）加窗频谱的特性

通过加窗后的滤波器频谱可以理解为理想滤波器频率响应与窗函数频谱的卷积：

$H(e^{j\omega}) = H_d(e^{j\omega}) \ast W(e^{j\omega})$

这带来的主要特性有：

旁瓣泄漏的降低：使用凯泽窗或其他改进窗函数（如汉明窗、布莱克曼窗等）可以显著提高旁瓣衰减，从而减少频谱泄漏，使得阻带性能更好。
过渡带的平滑化：由于窗函数本身频谱的主瓣宽度会叠加到滤波器的主瓣上，加窗后滤波器的过渡带宽度通常会比理想滤波器更宽，这是一种从理想到可实现化的折中。理想滤波器具有无限陡峭的过渡带，但实际设计中无法实现，窗函数的选择决定了实际过渡带的陡峭程度。
相位特性保持线性：大多数窗函数设计方法（包括凯泽窗）保留了理想滤波器的对称性，因而设计出的FIR滤波器具有线性相位特性，这是数字信号处理中非常重要的性质。

通过合理选择窗函数的参数（如凯泽窗中的 β 和滤波器的长度 N），可以在阻带衰减和过渡带宽度之间达到设计要求的平衡。

十、FIR型数字滤波器的频率采样设计法

1、预期频率特性的设置方法

在FIR滤波器设计中，首先需要明确预期的频率响应，也就是设计目标。这个目标通常涉及以下几个方面：

通带要求：在通带内，要求滤波器具有平坦的幅度响应（或满足特定幅度响应，如均衡器要求）和线性相位特性。
阻带要求：在阻带内，需要有足够的衰减以抑制不需要的频率成分，从而减小信号的干扰。
过渡带要求：通带和阻带之间的过渡带宽度及斜率，也是设计时需要平衡的重要指标。

数学描述

假设我们设计的FIR滤波器长度为 N，其理想（或预期）频率响应记为 $H_d(e^{j\omega})$ 。由于数字系统中频率是周期性的，通常我们只需要在区间 $[0, 2\pi)$ （或 $[-\pi,\pi]$ ）描述滤波器的频率响应。

为了利用频率采样法，我们在频率域上选取 N 个均匀分布的采样点，定义：

$\omega_k = \frac{2\pi k}{N},\quad k=0,1,\dots,N-1$

在这些采样点上的预期响应为：

$H[k] = H_d\left(e^{j\omega_k}\right)$

如果滤波器要求线性相位，那么滤波器的时域冲激响应 $h[n]$ 必须满足对称或反对称条件。对于对称型FIR滤波器，通常有：

$h[n] = h[N-1-n],\quad n=0,1,\dots,N-1$

这也是频率采样法设计中常采用的一种约束，确保设计出的滤波器具备线性相位特性。

2、频率采样法的设置过程

频率采样法设计FIR滤波器的基本思想是在频域上直接确定滤波器在一组离散频率点的取值，然后通过反离散傅里叶变换（IDFT）获得时域滤波器系数。具体步骤如下：

（1）频域采样及预期值确定

首先，在频域上均匀选择 N 个点：

$\omega_k = \frac{2\pi k}{N},\quad k=0,1,\dots,N-1$

并根据设计目标确定这些点的理想频率响应值 $H[k]$ 。例如，对于一个低通滤波器：

$H[k] = \begin{cases} 1, & \omega_k \in \text{passband} \\ 0, & \omega_k \in \text{stopband} \end{cases}$

注意在过渡带上可以采用渐变设计或直接设置为理想值，然后通过后续处理平衡波动。

（2）反离散傅里叶变换（IDFT）

得到频率采样值后，利用IDFT计算时域冲激响应：

$h[n] = \frac{1}{N}\sum_{k=0}^{N-1} H[k] e^{j\omega_k n},\quad n=0,1,\dots,N-1$

这一公式将频率域的离散采样点转换为时域系数。由此得到的 $h[n]$ 是严格满足 $H[k]$ 的有限长时域实现，同时也继承了线性相位的特性（前提是 $H[k]$ 具有共轭对称性或满足相应条件）。

（3）设计过程中的数学推导与注意事项

对称性条件：如果要求线性相位，则应确保频率响应采样值满足：
$H[k] = H^*[N-k],\quad k=1,\dots,N-1$
这样反变换后的 $h[n]$ 为实数且满足 $h[n]=h[N-1-n]$ 。
插值与过渡带：由于采样点有限，在过渡带区域的处理需要精心设计。直接取理想值可能引起振铃效应（Gibbs现象），因此在部分设计中会采用加权修正。
时域窗截断：与窗函数法类似，频率采样法在获得 $h[n]$ 后，如果希望设计更长的滤波器以获得更好的频率分辨率，也可以看作是一种“窗函数”截断方法的变形。此时，可以对 $h[n]$ 进行进一步加窗处理以改善旁瓣特性。

3、频率采样法的改进

虽然传统的频率采样法能直接通过IDFT得到满足预期频率响应的滤波器，但在边界频率处和过渡带区域往往会出现问题。为了解决这些问题，主要有以下两种改进方法。

（1）边界频率处采样值的修正

在频率采样法中，边界频率处的采样点（例如通带和阻带交界处）可能因为不连续性导致时域滤波器出现较大的波动。修正方法主要包括：

过渡平滑处理：在边界附近设计一段平滑过渡区域，使得 $H[k]$ 在边界处从1逐渐过渡到0，而非突然跃变。数学上可采用加权函数（例如余弦窗）使得边界附近的采样值满足：
$H[k] = \frac{1}{2}\left[1 + \cos\left(\pi \frac{\omega_k - \omega_b}{\Delta \omega}\right)\right]$
其中 $\omega_b$ 为边界频率， $\Delta \omega$ 为过渡带宽度。这样做可以减少IDFT后时域冲激响应的振铃现象。
采样点密度调整：在关键频率区域适当增加采样点的密度，再经过插值或加权平均，获得更平滑的边界采样值。理论上，这种方法等效于在频率上做非均匀采样，然后再映射到均匀采样点上。
数学修正公式：设原理想响应在边界处存在跳跃，通过修正函数 $f(\omega)$ 使得
$\tilde{H}_d(e^{j\omega}) = H_d(e^{j\omega}) \cdot f(\omega)$
其中 $f(\omega)$ 在边界处逐渐过渡（例如利用分段多项式或余弦函数），从而使得边界处采样值满足连续性要求。然后再以修正后的采样值进行IDFT计算。

（2）频率采样法和窗函数法的结合

另一种改进思路是将频率采样法与窗函数法结合，以利用两者的优点：

频率采样与加窗：首先按频率采样法得到初步时域冲激响应 $h[n]$ ，然后对 $h[n]$ 应用一个合适的窗函数 $w[n]$ （如凯泽窗、汉明窗等），得到：
$h_w[n] = h[n] \cdot w[n]$
这样既保留了频率采样法精确定义的频率响应，又通过窗函数抑制了时域截断带来的旁瓣泄漏问题。
组合设计：在频率采样过程中，可以考虑在设计频率响应采样值时就预留“窗函数效应”的余地，即在设定 $H[k]$ 时，将窗函数的频域影响考虑进去。例如，设计时采用加窗频谱与理想响应的卷积模型，对 $H[k]$ 进行预补偿，使得经过IDFT与后续加窗后更接近预期响应。
数学证明：假设理想频率响应为 $H_d(e^{j\omega})$ ，窗函数频谱为 $W(e^{j\omega})$ ，两者在时域卷积关系为：
$H(e^{j\omega}) = H_d(e^{j\omega}) \ast W(e^{j\omega})$
在频率采样法中直接取 $H[k]=H_d(e^{j\omega_k})$ 然后进行IDFT，效果相当于得到一个带有矩形窗效应的时域响应。通过预先对 $H[k]$ 进行修正，再乘以 $w[n]$ ，可以使得实际频率响应更接近理想设计目标，同时降低由于边界效应和截断引起的旁瓣振铃现象。
实际应用：在实际工程设计中，往往采用试凑法（trial-and-error）结合仿真优化滤波器响应。设计者可以先使用频率采样法得到初步滤波器，再通过调整窗函数参数（如窗长、窗函数类型、窗函数参数等）来获得更优的频率响应。

4、总结

本文详细分析了FIR数字滤波器的频率采样设计法，并给出以下几点要点：

预期频率特性的设计方法
- 明确通带、阻带及过渡带要求。
- 在频域上选取均匀采样点 $\omega_k = \frac{2\pi k}{N}$ ，并确定采样点上的理想响应 $H[k] = H_d(e^{j\omega_k})$
- 对于要求线性相位的设计，必须满足采样值的共轭对称性条件。
频率采样法的设计过程
- 通过IDFT将频域采样值转换为时域冲激响应：
  $h[n] = \frac{1}{N}\sum_{k=0}^{N-1} H[k] e^{j\omega_k n}$
- 强调了在设计过程中对称性条件和过渡带插值的必要性，以及可能存在的振铃现象（Gibbs现象）。
频率采样法的改进
- 边界频率处采样值的修正：采用平滑过渡（例如余弦平滑）、调整采样密度以及数学修正函数等方法，使得边界处的采样值连续，从而改善时域响应的振铃现象。
- 频率采样法与窗函数法的结合：先利用频率采样法获得初步时域响应，再对其加窗以抑制旁瓣泄漏，从而在频域上获得更接近预期设计目标的滤波器响应。

十一、最优等波纹FIR滤波器的设计

《参考白纸笔记：12、数字滤波器设计入门》里面的思路比这部分的思路要更有递进性。

1、最小最大问题的设计

最优等波纹FIR滤波器设计的基本思想是：在满足通带和阻带要求的前提下，使得滤波器的频率响应与理想响应之间的误差在各个频率处均匀分布，即最小化最大误差。这就是所谓的最小最大（minimax）设计问题，也称为切比雪夫逼近问题。

（1）问题描述

假设理想滤波器的幅度响应为 $D(\omega)$ （例如，对于低通滤波器，通带理想响应为1，阻带为0），而设计得到的滤波器幅度响应为 $H(\omega)$ 。定义误差函数为：

$E(\omega)=W(\omega)\bigl[D(\omega)-H(\omega)\bigr]$

其中 $W(\omega)$ 为加权函数，用以在不同频率区间内赋予不同的重要性。最小最大设计目标即为使得全频带的绝对误差的最大值最小：

$\min_{h[n]}\, \max_{\omega} \left|E(\omega)\right|$

（2）数学推导

在FIR滤波器设计中，滤波器的频率响应为

$H(\omega) = \sum_{n=0}^{N-1} h[n] e^{-j\omega n}$

最小最大问题可以理解为：在所有满足线性相位（或其他约束）的滤波器中，寻找一组 $h[n]$ ，使得误差函数在整个设计区间内的最大绝对值 $\|E(\omega)\|_\infty$ 达到最小。

利用切比雪夫逼近理论可以证明，最优解具有等波纹特性，即在预定的频带上，误差函数 $E(\omega)$ 的极值点数量达到上限，并且各极值值（绝对值）相等。换句话说，最优滤波器设计会使得误差在通带和阻带内呈现交替的极大和极小，且波纹幅值相同。

若记最优解的最大误差为 $\delta$ ，那么在最优条件下，对于一系列选定的频率点 $\omega_k$ 有：

$E(\omega_k) = \pm \delta, \quad k=1,2,\ldots,M$

其中 M 是极值点的个数。利用对偶理论和凸优化方法，可以将这一无限维问题转化为有限维的线性规划或凸优化问题来求解。

2、对极值数目的限制

最优等波纹滤波器设计中的一个关键理论是极值点定理。该定理指出，在最优解中，误差函数 $E(\omega)$ 必须在设计区间内出现至少 $L+2$ 个交替变化的极值点，其中 L 与滤波器的自由参数数目有关。

（1）自由参数与极值点个数

对于一个具有线性相位的FIR滤波器，其冲激响应通常具有对称性。如果滤波器长度 NNN 为奇数，则冲激响应满足：

$h[n] = h[N-1-n], \quad n=0,1,\dots,N-1$

这时自由参数个数为 $\frac{N+1}{2}$ ；若 N 为偶数，自由参数个数为 $\frac{N}{2}$ （或者稍有不同的计数方式，取决于具体约束）。设自由参数数目为 L+1 ，则理论上最优设计中误差函数 $E(\omega)$ 至少应出现 L+2 个交替极值。

（2）等波纹特性与交替定理

交替定理（Alternation Theorem）说明，如果在预定的设计区间内存在 L+2 个交替极值点，并且这些极值点的误差幅值相等，则当前解必定是最优解。反之，如果解不是最优的，则必定存在超过 L+2 个交替极值。

数学上，设误差函数在一组频率点 $\{\omega_0,\omega_1,\ldots,\omega_{L+1}\}$ 处满足：

$E(\omega_k) = (-1)^k \delta, \quad k=0,1,\dots,L+1$

这一条件不仅是最优解的必要条件，也可以视为最优化问题的约束之一。利用这一交替性质，可以对最优滤波器设计问题构造一组方程，通过这些方程求解最优滤波器系数 $h[n]$ 和最优误差 $\delta$ 。

3、Parks-McClellan算法

Parks-McClellan算法（也称为Remez交换算法）是一种实现最优等波纹FIR滤波器设计的数值方法。该算法利用交替定理的理论基础，将最小最大设计问题转化为寻找具有等波纹性质的误差函数，从而达到全局最优。

（1）算法基本思想

Parks-McClellan算法主要步骤如下：

初始选取极值点：在设计区间内均匀选取 L+2 个频率点作为初始的交替点。
构造插值方程：利用这些极值点构造一组关于滤波器系数和最优误差 δ\deltaδ 的线性方程组。对于线性相位FIR滤波器，可写为：
$D(\omega_k) - \sum_{n=0}^{L} a[n] \cos(\omega_k n) = (-1)^k \delta,\quad k=0,1,\dots,L+1$
其中 $a[n]$ 是由原始系数 $h[n]$ 通过对称性关系提取的独立参数。
求解方程组：求解上述线性方程组得到滤波器系数和当前的误差幅值 $\delta$ 。
交换步骤：利用当前解计算整个设计区间内的误差函数，寻找新的极值点。如果新选取的极值点与原先设定的不一致，则更新极值点集，并返回步骤2进行迭代。算法不断交换极值点，直到满足收敛条件，即误差函数在所有交替点处波纹幅度相等且交替顺序保持不变。

3.2 数学推导

设设计区间为 $\omega \in \Omega$ ，理想响应 $D(\omega)$ 和加权函数 $W(\omega)$ 定义如前，令设计参数为 $a[n]$ （经过对称性约束后独立的系数个数为 L+1 ）。在交替点 $\omega_k$ 处，我们有插值方程：

$D(\omega_k)W(\omega_k) - \sum_{n=0}^{L} a[n] \cos(\omega_k n) = (-1)^k \delta,\quad k=0,1,\dots,L+1$

构成 L+2 个方程，未知数为 L+1 个滤波器参数 $a[n]$ 以及 $\delta$ 共 L+2 个变量。该方程组通常写成矩阵形式：

$\begin{bmatrix} \cos(0\cdot\omega_0) & \cos(1\cdot\omega_0) & \cdots & \cos(L\cdot\omega_0) & 1 \\ \cos(0\cdot\omega_1) & \cos(1\cdot\omega_1) & \cdots & \cos(L\cdot\omega_1) & -1 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ \cos(0\cdot\omega_{L+1}) & \cos(1\cdot\omega_{L+1}) & \cdots & \cos(L\cdot\omega_{L+1}) & (-1)^{L+1} \end{bmatrix} \begin{bmatrix} a[0] \\ a[1] \\ \vdots \\ a[L] \\ \delta \end{bmatrix} = \begin{bmatrix} D(\omega_0)W(\omega_0) \\ D(\omega_1)W(\omega_1) \\ \vdots \\ D(\omega_{L+1})W(\omega_{L+1}) \end{bmatrix}$

求解上述方程组可得到滤波器参数和当前的极值误差 $\delta$ 。

3.3 交换与收敛

在求解了参数后，利用当前滤波器设计计算整个频率区间的误差函数：

$E(\omega)=W(\omega)\Bigl[D(\omega)-\sum_{n=0}^{L} a[n] \cos(\omega n)\Bigr]$

算法在整个区间上寻找新的交替极值点集。如果新选取的交替点与原来的不同，则用新点集重新构造方程组，重复迭代，直到交替极值点的分布稳定，即所有极值点处的误差交替并且幅度趋于一致。交替定理保证了当满足 L+2 个交替极值条件时，解达到最优。

十二、FIR滤波器设计中的matlab函数

在MATLAB中自带了FIR滤波器设计的窗函数法、频率采样法和最佳等波纹法相应的函数。

1、与窗函数相关的函数

（1）窗函数（窗函数设计）

`1.1 boxcar`（矩形窗 Boxcar Window）

函数调用：wind = boxcar(N)
说明：N为窗函数的长度，通常与FIR滤波器阶数相关，且窗函数长度是M+1。

`1.2 hanning`（汉宁窗 Hanning Window）

函数调用：wind = hanning(N)
说明：汉宁窗是一种常见的窗函数，适用于多种滤波器设计。

`1.3 hamming`（海明窗 Hamming Window）

函数调用：wind = hamming(N)
说明：海明窗是另一个广泛使用的窗函数，类似于汉宁窗，但有稍微不同的频率响应。

`1.4 blackman`（布莱克曼窗 Blackman Window）

函数调用：wind = blackman(N)
说明：布莱克曼窗用于抑制旁瓣，适用于需要较好频率响应的滤波器设计。

`1.5 kaiser`（凯泽窗 Kaiser Window）

函数调用：wind = kaiser(N, beta)
说明：凯泽窗具有可调节的参数beta，能在不同的设计要求下调整窗函数的形状，适用于需要高精度滤波器设计。

（2）`kaiserord`（求解凯泽窗相关的FIR滤波器参数）

函数调用：
- [n, Wn, beta, ftype] = kaiserord(f, a, dev)
- [n, Wn, beta, ftype] = kaiserord(f, a, dev, fs)
  说明：
- f是频率数组，a是幅度数组，dev是最大允许的误差。
- f 和 a 是频率和滤波器幅值矢量，当设计低通滤波器时，f=[fp fs]，a=[1 0]；当设计带通滤波器时，f=[fs1 fpl fp2 fs2]，a=[0 1 0]，对高通和带通滤波器可按低通和带通给出。f 的长度是2*length(a)-2。
- dev也是矢量,包含有通带的波纹和阻带的衰减(线性值，如果给的是分贝值的话需要转化为线性值)，其长度与 a 等长。fs 是采样频率(单位为Hz)，若没给fs ，则 f 中的各值必须是归一化的频率值。输出参数中 n 是滤波器的阶数，Wn是滤波器的频率参数，beta是kaiser窗函数中的参数，ftype是滤波器的类型,在fir1函数调用时将会用到。

窗函数使用总结：

在MATLAB中，窗函数常用于设计FIR滤波器。窗函数的长度通常是滤波器阶数加1（M+1），与滤波器设计的特性密切相关。
每种窗函数（如短形窗、汉宁窗、凯泽窗等）都有其特定的适用场景和性能表现，选择合适的窗函数对于滤波器设计至关重要。

2、设计FIR的函数

（1）`fir1` （用窗函数法设计FIR滤波器）

fir1 函数

功能：用窗函数法设计FIR滤波器，求得滤波器系数 B 序列。
调用格式：
- B = fir1(n, Wn)：设计一个低通FIR滤波器。
- B = fir1(n, Wn, 'ftype')：根据 ftype 设计不同类型的FIR滤波器（如'low'、'high'、'bandpass'等）。
- B = fir1(n, Wn, Window)：使用指定的窗函数Window设计FIR滤波器。
- B = fir1(n, Wn, 'ftype', Window)：使用指定的窗函数Window设计FIR滤波器。
说明：
firl函数是以经典方法实现加窗线性相位FIR数字滤波器设计,它可设计出低通带通、高通和带阻滤波器。
通过 B=fir1(n,Wn)形式可得到n阶FIR低通滤波器,n是波器阶数,Wn是滤波器带宽。滤波器系数包含在B中,可表示为

$b(z) = b(1) + b(2)z^{-1} + b(3)z^{-2} + \cdots + b(n+1)z^{-n}$ ，这构成一个用海明窗函数截止频率为Wn的线性相位滤波器,0≤Wn≤1.目Wn=1时对应于 0.5fs。
当Wn=[W1 W2] 时，fir1 函数这种形式可得到带通滤波器,其通带为 W1<ω<W2
通过 B=fir1(n,Wn,'ftype') 形式可设计高通和带阻滤波器，由ftype 决定
B=firl(n,Wn,Window) 利用列矢量 Window中指定某种窗函数进行滤波器设计(Window可以是任一种)，Window长度为n+1，如果不指定 Window参数，则firl函数在默认时采用海明窗。

（2）`fir2`（用频率采样法设计FIR滤波器）

fir2 函数

功能：使用频率采样法设计FIR滤波器。
调用格式：
```
B = fir2(n, f, a, window)
```
其中：
- n：FIR滤波器的阶数。
- f：频率点的数组，表示设计滤波器的频率范围。
- a：对应的幅度响应数组，表示频率点f上对应的幅度值。
- window：可选参数，表示窗口函数，用于优化滤波器性能。
说明：
- 设计低通滤波器：当f = [fp fs]时，a = [1 0]，表示设计低通滤波器，其中fp是通带截止频率，fs是阻带开始频率。
- 设计高通滤波器：当f = [fs fp]时，a = [0 1]，表示设计高通滤波器，其中fs是通带开始频率，fp是阻带截止频率。
- 设计带通滤波器：当f = [fs fp f1 f2]时，a = [0 1 0]，设计带通滤波器，fs和f1为带通频带的边界频率，f2和fp为带阻频带的边界频率。
- 归一化频率：f的值需要归一化，即频率值应该在[0, 1]之间，0表示0 Hz，1表示Nyquist频率（采样频率的一半）。
窗口函数：
- 可以通过window参数选择不同的窗函数（如hann、hamming、blackman、kaiser等）来优化FIR滤波器的性能。
- 默认情况下，fir2函数使用海明窗（Hamming window）。
应用场景：
- 适用于各种类型的滤波器设计（如低通、高通、带通、带阻滤波器），通过频率采样法进行设计。
- 适合于需要精确频率响应的场合。

这个函数在设计具有特定频率响应的FIR滤波器时非常有用，特别是当你需要指定不同频率点的幅度响应时。

（3）`firpmord`（计算等纹波FIR滤波器的阶数）

firpmord 函数

功能：计算FIR滤波器的阶数。
调用格式：
```
[n, fo, ao, w] = firpmord(f, a, dev)
[n, fo, ao, w] = firpmord(f, a, dev, fs)
```
- f：频率点数组，定义滤波器的频率响应。
- a：幅度响应数组，定义频率点f的幅度值。
- dev：误差值，表示滤波器与理想滤波器之间的偏差的矢量,包含有通带的波纹和阻带的衰减(线性值),其长度与a等长。
- fs：采样频率（单位为Hz）。若没给fs，则 f 中的各值必须是归一化的值。
说明：
- f和a分别是频率和幅度数组，f的长度为2 * length(a) - 2，f的值必须归一化到[0, 1]。
- dev是期望误差，表示实际滤波器与理想滤波器之间的误差范围。fs是采样频率，若未给定fs，则f的值需要进行归一化。
- n是计算出来的滤波器阶数，fo是频率点，ao是对应的幅度值，w是权重值。
应用场景：
- 用于根据给定的频率响应和期望误差估算FIR滤波器的阶数。
- 适用于设计复杂滤波器并根据设计要求调整阶数。

（4）`firpm`（等纹波FIR滤波器设计）

firpm 函数（Parks-McClellan法设计）

功能：使用Parks-McClellan方法设计FIR滤波器。
调用格式：
```
b = firpm(n, f, a) 
b = firpm(n, f, a, w) 
b = firpm(n, f, a, 'ftype') 
b = firpm(n, f, a, w, 'ftype')
```
- n：FIR滤波器的阶数。
- f：频率点数组，定义滤波器的频率响应。
- a：幅度响应数组，定义对应频率点的幅度值。
- w：权重数组，用于调整各个频率响应的误差。
- 'ftype'：滤波器类型，指定滤波器的类型（如低通、高通、带通等）。
说明：
- 在MATLAB新版本中，firpm函数替代了旧版本中的remez函数。f和a分别是频率和幅度的数组，f的长度为2 * length(a) - 2，而f的频率必须归一化到[0, 1]范围内。
- w是权重数组，表示每个频率点的误差偏差。可以通过指定w来控制滤波器的误差。
- ftype是滤波器类型，可以选择'low'（低通滤波器）、'high'（高通滤波器）、'bandpass'（带通滤波器）、'stop'（带阻滤波器）等类型。
应用场景：
- 适用于需要通过最小化误差设计FIR滤波器的情况，特别是当设计带通、带阻或其他复杂类型的滤波器时。
- 适用于优化滤波器性能并准确控制频率响应的场合。
- 这些函数是用于FIR滤波器设计中的高级方法，特别是在复杂的滤波器设计中，使用Parks-McClellan方法进行优化设计时非常有效。

3、另一些有用的与FIR设计有关的函数

有一些MATLAB的函数不是MATLAB自身带的，而是由其他一些书本编制的，在案例程序中，对于FIR滤波器设计和应用很有用，介绍如下。注：这些函数已包含在案例程序包的basictbxsp 子目录中。

（1）`ideal_lp`（理想低通滤波器）

功能：此函数用于生成一个理想低通FIR滤波器的冲激响应。
调用格式：
```
hd = ideal_lp(wc, M)
```
- wc：归一化截止频率，公式为 wc = fc * π / (0.5 * fs)，其中 fc 是截止频率（单位为Hz），fs 是采样频率。
- M：滤波器的长度（即滤波器的阶数）。
说明：该函数返回一个理想低通滤波器的冲激响应。wc是归一化的截止频率，取值范围是 0 到 π。
举例：
- 低通滤波器：hd = ideal_lp(wc, M)
- 高通滤波器：hd = ideal_lp(pi, M)
- 带通滤波器：hd = ideal_lp(wc1, M) + ideal_lp(wc2, M)，其中 wc1 和 wc2 是带通滤波器的低、高截止频率。

（2）`hr_type1`（第1类线性相位FIR滤波器的振幅响应）

功能：该函数用于计算第1类线性相位FIR滤波器的振幅响应。
调用格式：
```
[Hr, w, a, L] = hr_type1(h)
```
- h 是FIR滤波器的冲激响应序列。
- 输出 Hr 是滤波器的振幅响应，w 是频率点，a 和 L 为相关的参数。

（3）`hr_type2`（第2类线性相位FIR滤波器的振幅响应）

功能：该函数用于计算第2类线性相位FIR滤波器的振幅响应。
调用格式：
```
[Hr, w, b] = hr_type2(h)
```
- h 是FIR滤波器的冲激响应序列。
- 输出 Hr 是滤波器的振幅响应，w 是频率点，b 为相关的参数。

（4）`hr_type3`（第3类线性相位FIR滤波器的振幅响应）

功能：该函数用于计算第3类线性相位FIR滤波器的振幅响应。
调用格式：
```
[Hr, w, c] = hr_type3(h)
```
- h 是FIR滤波器的冲激响应序列。
- 输出 Hr 是滤波器的振幅响应，w 是频率点，c 为相关的参数。

（5）`hr_type4`（第4类线性相位FIR滤波器的振幅响应）

功能：该函数用于计算第4类线性相位FIR滤波器的振幅响应。
调用格式：
```
[Hr, w, L] = hr_type4(h)
```
- h 是FIR滤波器的冲激响应序列。
- 输出 Hr 是滤波器的振幅响应，w 是频率点，L 为相关的参数。
说明：该函数的输入参数 h 是FIR滤波器的冲激响应，输出 Hr 是滤波器的振幅响应。

（6）`ampl_res`（FIR滤波器的振幅响应）

功能：该函数用于计算FIR滤波器的振幅响应。
调用格式：
```
[Hr, w, P, L, type] = ampl_res(h)
```
- h 是FIR滤波器的冲激响应序列。
- 输出 Hr 是振幅响应，w 是频率点，P 是频率分量，L 是滤波器的长度，type 是滤波器类型。
说明：该函数还可以根据提供的冲激响应序列 h，计算其对应的振幅响应。

这些函数帮助设计和分析FIR滤波器的频率响应、振幅响应等特性，尤其在进行线性相位滤波器设计时非常有用。

十三、FIR滤波器设计案例

案例3.15、在窗函数法设计FIR中如何选择窗函数和阶数N

案例3.16、用 ideal_lp 函数和 fir1 函数设计的滤波器是否相同

案例3.17、用凯泽窗设计FIR滤波器的优点

案例3.18、为什么FIR滤波器不适用于设计数字陷波器

案例3.19、通过FIR滤波器的输出，延迟量如何校正

案例3.20、通过fir2函数设计任何响应的FIR滤波器

案例3.21、通过 firpm函数设计的FIR 滤波器为什么达不到指标要求

案例3.22、如何设计多频带的FIR滤波器

案例3.23、如何用FIR滤波器设计数字微分器

案例3.24、如何用FIR滤波器设计数字希尔伯特变换器

十四、用FDATool设计数字滤波器

1、IIR滤波器设计

2、FIR滤波器设计

3、SOS系数的进一步说明

SOS称为二阶分割滤波器系数，它不同于我们以前讲述的一般的滤波器系数，这种系数形式只使用在IIR滤波器中。
在3.1节中曾介绍过滤波器的级联形式(如图3-1-5所示)。由FDAToI(与3.15节中介绍的fdesign+design方法相同)设计的IIR滤波器常表示为如图3-14-22所示的形式，即每一个滤波器是由L个二阶滤波器级联而成，其中每一个Hi(z)(i=1,2,…,L)都是二阶滤波器,可表示为

$H(z) = \frac{b_{i0} + b_{i1} z^{-1} + b_{i2} z^{-2}}{a_{i0} + a_{i1} z^{-1} + a_{i2} z^{-2}}, \quad a_{i0} = 1$

式中， $a_{i0}$ 恒为1.

4、案例

案例3.25、如何把SOS或Hd转变为滤波器的系数

十五、用fdesign和design设计数字滤波器

1、使用fdesign和design设计数字滤波器的函数介绍

在MATLAB中除了3.6节、3.12节和3.14节介绍的MATLAB函数和工具设计IIR和FIR数字滤波器外，还有一对函数fdesign和design也可以用来设计数字滤波器。这一对函数既可用来设计IIR滤波器，也可以用来设计FIR滤波器。它们和FDATooI一样，但不是以图形形式设计，而是以编程的形式。

函数：

低通滤波器设计函数：d = fdesign.lowpass

高通滤波器设计函数：d = fdesign.highpass

带通滤波器设计函数：d = fdesign.bandpass

带阻滤波器设计函数：d = fdesign.bandstop

带陷滤波器设计函数：d = fdesign.notch

注：函数的参数不再详细分析

功能：

给出滤波器参数指标集合d

滤波器设计：

以上由fdesign给出的 d 只是滤波器参数的集合，还不是滤波器的系数，要通过design函数才能把滤波器参数集合d转换成滤波器的系数集合h。

名称：design
功能：计算出滤波器系数集合

        调用格式：
                h= design(d)
                h = design(d,designmethod)
        说明：输人参数d是由fdesign函数给出的参数集合;designmethod是给出IIR或FIR 滤波器具体某一种设计方法,具体方法列于表3-15-6中。

在应用designmethod的具体值时都要用单引号括起来。输出参数h是滤波器系数的集合，是一个结构型数组。

2、案例

案例3.26、为什么在使用 design函数时常会出现“invalid design method’。

案例3.27、用 fdesign+design 的方法与前几节介绍的经典方法设计的滤波器是否相同

案例3.28、用 fdesign+design 方法有什么优点

十六、三分之一倍频程滤波器

1、三分之一倍频程滤波器的基本概念

三分之一倍频程滤波器是将整个频谱按照对数（即频率比而非绝对差值）均匀划分成多个频带，每个频带的宽度正好是一个倍频程（即频率翻倍）宽度的三分之一。下面我将从几个角度来解释这一概念：

（1）什么是倍频程？

倍频程（Octave）：指的是在一个频带中，最高频率是最低频率的两倍。比如，一个倍频程从100 Hz到200 Hz；下一个倍频程从200 Hz到400 Hz。
这种划分方式符合人耳对频率的感知（即对数感知），因为我们对频率变化的敏感度不是线性的，而是对比例更敏感。

（2）三分之一倍频程是什么意思？

三分之一倍频程：就是把一个倍频程再平均分成三个相等的部分。换句话说，每个频带的频率范围不是固定的绝对值，而是满足上限频率与下限频率之比相同。
举个例子，如果有一个中心频率为 $f_c$ 的滤波器，它的下界和上界可以通过下面的公式确定：
$f_l = \frac{f_c}{2^{1/6}},\quad f_u = f_c \cdot 2^{1/6}$
这里的 $2^{1/6}$ 大约为1.122。这意味着，从下界到中心频率，再到上界，频率每增加一段都是乘以1.122，从而在对数尺度上均匀分布。

（3）为什么使用三分之一倍频程滤波器？

更高的分辨率：相比直接使用倍频程滤波器（即频率翻倍），三分之一倍频程滤波器提供了更细致的频率划分，能更准确地分析信号中各个频段的能量分布。
广泛应用于音频和声学测量：例如在噪声监测、音乐调音等场景中，细分的频段能帮助更好地理解和处理不同频率成分。

（4）实际例子

假设中心频率为1000 Hz：

下界： $1000 / 2^{1/6} \approx 1000 / 1.122 \approx 891$ Hz
上界： $1000 \times 2^{1/6} \approx 1000 \times 1.122 \approx 1122$ Hz

这样，滤波器主要对891 Hz到1122 Hz之间的信号起响应，其他频率的信号则被衰减。

总的来说，三分之一倍频程滤波器通过在对数尺度上对频谱进行细分，使得每个滤波器的带宽（频率范围）与中心频率呈固定比例关系，从而实现更加精细和均匀的频谱分析。

2、数学模型与参数计算

（1）带宽与品质因数

带宽（BW）：对于三分之一倍频程滤波器，带宽可以表示为：
BW=fu−fl=fc(21/6−2−1/6)BW = f_u - f_l = f_c \left(2^{1/6} - 2^{-1/6}\right)BW=fu−fl=fc(21/6−2−1/6)
品质因数（Q）：品质因数用于衡量滤波器的选择性。对于理想的带通滤波器，
Q=fcBW=fcfc(21/6−2−1/6)=121/6−2−1/6Q = \frac{f_c}{BW} = \frac{f_c}{f_c \left(2^{1/6} - 2^{-1/6}\right)} = \frac{1}{2^{1/6} - 2^{-1/6}}Q=BWfc=fc(21/6−2−1/6)fc=21/6−2−1/61
这个值是固定的，与中心频率无关。计算后大约为4.32左右。高Q值表示滤波器的带宽较窄，选择性较高。

（2）频率响应

设计滤波器时，我们通常关注幅频响应和相位响应：

幅频响应：描述在各个频率上滤波器对信号幅度的衰减情况。理想的三分之一倍频程滤波器在其中心频率处幅值最高，并随着频率远离中心逐渐衰减。
相位响应：反映滤波器对信号各频率分量的相位延迟。对于一些应用（如高保真音频处理），线性相位是十分重要的，但通常线性相位FIR滤波器的计算量较大。

3、滤波器设计与实现方法

（1）FIR与IIR滤波器设计

FIR滤波器：有限脉冲响应滤波器的设计方法包括窗函数法、频率抽样法、最小二乘法等。FIR滤波器可实现严格的线性相位，但往往需要较高的阶数以达到理想的滤波器特性。对于三分之一倍频程滤波器，如果要求相位特性好，可以考虑使用FIR设计。
IIR滤波器：无限脉冲响应滤波器采用递归结构，常见设计方法有Butterworth、Chebyshev或椭圆滤波器设计。这类滤波器阶数通常较低，计算效率高，但相位响应通常不是严格线性。常用于实时处理或对相位要求不严格的场合。