1.2.3 数字信号处理部分知识重温(3) FIR滤波器的设计

写在前面：本文不作为一篇严谨的理论推导文章，仅为博主学习时的随手小记，个人的主观认知为主，以求日后能在较快时间内理清思路，找到当下的学习状态。文章中会尽量避免出现过多公式及推导，如果能有幸给其他朋友带来一些帮助，不甚荣幸。

感谢以下文章及老师的指导：

【数字信号处理(北京交通大学 陈后金)】 2.1 离散傅里叶变换（DFT） - 4 DFT性质2_哔哩哔哩_bilibili

因为FIR滤波器设计更为简单，在工程中更常见，实现更简单，因此，对滤波器的学习便到此为止了，之后在项目中也主要使用这种类型的滤波器。IIR滤波器日后有机会再学习。

(重要！)引言：

我们想要得到一个理想滤波器的幅频特性（比如在频域表现为一个矩形，在通带幅度为1，阻带为0，无过渡带），但实际上是无法实现的，需要实现这个要求的滤波器响应时域序列是无限长的。（实际上，用户提出的需求大多也是理想型的，无法物理实现）

因此我们以拥有线性相位的FIR滤波器来拟合理想滤波器，得到近似的幅频特性。根据理想滤波器得到FIR滤波器的途径有三种：①理想滤波器时域序列*在一定范围内取值不为0的序列（又叫窗函数WN[k]），从而将理想滤波器序列截断，截断序列即为FIR滤波器时域响应序列。

我们向用户交付的产品实际上是近似满足用户需求，可物理实现的FIR滤波器

• FIR滤波器的引入

当一个函数（序列）进入系统时，它的时域表达式可以表示为信号函数与系统响应函数的卷积。在频域，Z域表示为信号函数与系统响应函数相乘，为了计算的便利，我们通常在频域进行研究。

如下图，离散线性时不变系统在Z域的系统响应函数如下图一第一行，在ai=0的特殊条件下，H(z)可化简为下图一第二行左式，我们参考常见信号Z域变换公式(下图二），可逆变换出系统在时域的响应序列（下图一第二行右式）,序列元素为H(z)分子系数bi，参见表达式，bi总共有M+1个。可见响应序列拥有有限的序列长度，我们称这样的系统为FIR系统。

系统特性：1.稳定（所有极点的模r都小于1）

1. 系统可为线性相位
2. 只有一个零极点，一定在单位圆内

上图一是线性相位的系统工作示意图，线性相位的含义，是对于不同的输入信号，从系统输出时，经历的时延量是一致的。由此可以保证，各类复杂的，有多种简单信号合成的信号经过系统，由于其构成的所有波形时延一致（节奏一致），输出的复杂信号波形不变，也拥有一个相同的时延。

上图二是具体的系统频域响应函数、相位函数、时延等参数的具体分析。

这是一个非线性相位的反例，对于不同信号的延迟不一致，各类复杂的，有多种简单信号合成的信号经过系统，由于其构成的所有波形时延不一致（节奏不一致），输出的复杂信号波形就会失真。

大家可以想象一下趣味运动会中的两人三足运动，每个成员行进节奏不一致，是不是整体来看，就会东倒西歪，跌跌撞撞。

• FIR滤波器设计参数

如下图一左侧，我们可以看到，对一个FIR滤波器进行（z）频域分析，得到的图像如下，我们主要关注的参数见下图一右侧：

（1）G（信号的衰减量或增益量）-单位：dB

负数的分贝值表示信号的衰减，而正数的分贝值表示信号的增益

     G=-20log10(信号幅度)

（2）Ap（通带最大衰减）和As（阻带最小衰减）：

- Ap（通带最大衰减）：Ap表示通带中最大可接受的衰减值（对应但不等于通带幅度的最小值）

- As（阻带最小衰减）：As表示阻带中最小需要达到的衰减值（对应但不等于阻带幅度最大值）

Ap和As通常是根据特定的应用需求和性能要求来确定的。例如，在音频信号处理中，需要保持通带幅度响应平坦，并且需要有效地滤除杂波和噪声。因此，滤波器的设计需要尽可能地实现最小化的Ap和最大化的As，以获得最佳的性能和效果。

计算公式如下：

- 通带最大衰减 Ap = -20log10(1 - δp)，1 - δp即为对应但不等于通带幅度的最小值，其中δp为通带最大允许的波纹系数，Ap 通常取值为0.01或0.05。

- 阻带最小衰减 As = -20log10(δs)，δs即为阻带幅度最大值，为阻带最小需要达到的衰减值，As通常取值为40dB或60dB。

需要注意的是：

1.Ap与As为衰减量，但与信号幅值衰减变化量无关，与当前信号幅值有关 

2.Ap与As是G的特殊情况

      当我们规定好这些参数之后，便需要依靠这些指标推导滤波器的（z域）频域系统函数，主要的方法有：①窗函数法 ②优化设计法 ③频率取样法

      我们设计FIR滤波器，最终向用户交付的就是滤波器的（z）频域系统函数，时域系统函数（IIR不可以）。

• 线性相位FIR滤波器的特性

  ①线性相位的充要条件

    如下图一，FIR滤波器的线性相位指代的是广义线性相位

要满足这样的线性相位，需要满足下图的充要条件：

时域的充要条件和z域的充要条件是等价的，可以相互转换的：时域位移M单位转换的z域就（*z^-M）；h[k]是一串冲激，根据常见函数序列z域转换表（上面已经放过了）：∑h[ki]*δ(t-ki*T)↔∑h[ki]*z^(-1*ki)=H(z^-1）  (h[ki]视作常数)

  从这我们也可以得到：N为滤波器系统时域响应函数（序列）长度，M为FIR滤波器阶数，M=N-1；比如拥有长度为3的时域序列的滤波器阶数为2。

②FIR滤波器时域特性

      时域序列呈偶对称：(属于1,2型)

      时域序列呈奇对称：（属于3,4型）（3型中心点必为0）

    另：时域序列长度为奇数的必为1,3型，长度为偶数的必为2,4型

    ③FIR滤波器频域特性：

      （1,2型）

  中点：序列首尾序列点坐标之差/2

A（Ω）=∑序列点坐标取值*cos(（该序列点坐标与中点坐标差值）Ω)

φ（Ω）=(-M/2)*Ω=-(N-1)/2*Ω= -1*中点坐标*Ω （严格线性相位）

1型：序列长度N为奇数，呈偶对称

A（Ω）关于原点偶对称、π点偶对称

该型滤波器可设计低通、高通、带通、带阻滤波器

2型：序列长度N为偶数，呈偶对称

  A（Ω）关于原点偶对称、π点对称奇对称

  不能设计高通、带阻滤波器！

  如上图，通过分析A（Ω）图像可以发现，由于A（Ω）对π点奇对称，故在π点处A（Ω）=0，而-π—π是系统的工作频率，π点处是系统的最高工作频率（高频部分），在此处幅度函数A（Ω）=0，意味着，2型滤波器不能实现带阻、高通滤波器

        （3、4型）

中点：序列首尾序列点坐标之差/2

A（Ω）=∑序列点坐标取值*sin(（该序列点坐标与中点坐标差值）Ω)

φ（Ω）=(-M/2)*Ω+π/2=-(N-1)/2*Ω+π/2（广义线性相位）

（注：相位出现π/2的原因将结合3型滤波器在下文解释）

3型：序列长度N为奇数，呈奇对称

           在原点处奇对称，在π处奇对称

           仅能制作带通滤波器！！

（欧拉公式：sinx=—j/2*e^jx + j/2*e^—jx，—j=e^（—j*（π/2）），故（—j/2）*（e^—jx — e^—jx）=sinx==>e^（—j*（π/2））*（e^—jx — e^—jx）=2sinx）

4型：序列长度N为偶数，呈奇对称

  在原点处奇对称，在π点处偶对称

  不能设计带阻、低通滤波器！！

      总结：如下图，我们还可以总结出一些关键信息：

        1.滤波器响应序列h[k]是什么类型的对称，其幅频曲线就在原点呈什么类型的对称

2.h[k]长度N为偶数时，幅频曲线在π点的对称性与h[k]的对称性相反；h[k]长度N为奇数时，幅频曲线在π点的对称性与h[k]的对称性相同。

④FIR滤波器零点分布：

 零点情况一：z=1/zi（倒数）

  前文我们提及了：要想满足线性相位，则必须满足的充要条件，在时域表示为时域序列满足奇对称或偶对称；在Z域表示为：

         我们知道，一个系统的频谱，在其零点zi处，H(zi)=0，z是个复数，是不为0的量，故H（1/zi)=0

 零点情况二：z=z*（共轭）（前面我们已经分析过，FIR数字滤波器的频域表达式为H（z）=∑bi*z^-i（i=0,1,2...）（见本文第一张图）=b0+∑bi*z^-i（i=1,2...)，设∑bi*z^-i（i=1,2,3,4...)=a（z）（实数部分）+b（z）*j（虚数部分），假设零点为zi，则a（zi）=-b0，b（zi）=0。相应的，zi取共轭时，实数系数不变，虚数系数变为相反数，a（zi*）=a（zi）=-b0，b（zi*）=-b(zi)=0，可见z*依旧为零点）

          零点情况三：结合情况一、情况二，可知z=(1/zi)*也是零点

综上所述，一个零点，伴随着会出现另外三个零点

零点分布情况①：

  零点分布情况②：

  零点分布情况③：

  零点分布情况④：

总结：

  1.

1. 关于零点的分布情况无需死记硬背，只需记住z=r*e^jθ即可，当我们得知一个零点的确切值时，自然可以通过将r取倒数，θ取相反数等方式找出其他伴随的几个零点

例题0：

四、窗函数法设计FIR滤波器

  1.设计思想：（Hd（e^jΩ）是理想达到的设计指标，通常无法完全达到）

注：用户给定的Hd（e^jΩ）=e^-j*φ1（Ω）*A1（Ω）通常是以频域图像的形式给出，我们需要通过图像自己建立拟合理想滤波器幅度的FIR滤波器幅度函数A2（Ω），同时选定用来拟合的FIR滤波器类型及阶数，相应的得出相位函数φ2（Ω），根据Hd（e^jΩ）=(e^-j*φ2（Ω）)*A2（Ω）得出准确的数学表达式

*注：1.（φ（Ω）不来自用户给定的滤波器频域图像，是因为用户给定的图像不一定能满足线性相位的条件，故用我们自己设计的FIR滤波器的φ2（Ω）去替换）（大概理解，知道怎么操作就行）

     2.A2（Ω）与A1（Ω）可能存在一定程度不同，比如3,4型滤波器要求滤波器幅频特性奇对称，因此在正半轴的幅度表达式是负半轴的幅度表达式的相反数；但1,2型时正负半轴表达式相同，A2（Ω）=A1（Ω）

     3.滤波器若选择1,3型，则滤波器阶数必须为偶数；2,4型滤波器阶数必须为奇数。

     4.Hd（e^jΩ）=e^-j*φ1（Ω）*A1（Ω）=e^-j*φ2（Ω）*A2（Ω）计算的结果始终是理想滤波器的频谱，那我们用e^-j*φ2（Ω）*A2（Ω）替换e^-j*φ1（Ω）*A1（Ω）去表达理想滤波器频谱的意义是什么呢？答案是保证了理想滤波器在一开始就一定满足线性相位的条件，这样才能转化出FIR滤波器

加窗截断hd[k]中的“窗”，全名叫窗函数，符号为WN[k]（在一个区间内取值不为0）。加窗截断实际上就是hd[k]*WN[k]，将hd[k]的序列截断在WN[k]取值不为0的区间内。hd[k]是理想滤波器的时域响应序列，hd[k]*WN[k]才是FIR滤波器的时域响应序列。

    至此，我们终于得到了FIR滤波器。

以下是窗函数WN[k]的频域图像：

    与sa函数有些类似，窗函数频域也有若干个峰，靠近原点的第一个峰是“主瓣”，其余的峰都是“旁瓣”。

主瓣和旁瓣有如下性质：

1. 主瓣越窄/窗函数N越大，过渡带越小
2. 旁瓣越小，在阻带衰减越大

  例子1：

在本题中，我们引入了第一种窗函数称为矩形窗，长度为N=M+1，幅度为1，以下图示是理想滤波器系统响应时域序列加窗之后，得到的逼近拟合理想时域序列的可实现序列，这个才是我们设计的FIR滤波器的系统响应时域序列。

下图给出了不同阶FIR滤波器（不同长度的系统响应时域序列/加不同长度的窗）的幅频特性曲线

下图将幅度转换为了信号的衰减量（G=-20log10（幅度））表示

1. 第一种类型：矩形窗

  注：主瓣宽度是将主瓣正半轴和负半轴两部分算进去了，我们实际可用的只是正半轴这一半（负频率实际都是不存在的）

4.第二种类型：加权窗

1.  实际上就是在矩形窗的基础上*一个可变系数（函数）

1. 与矩形窗相比，由于加权窗的初始和末端序列值比较平滑的衰变，因此加权窗的高频分量小，旁瓣分量小（阻带衰减大）
2. 代价是牺牲了主瓣宽度（主瓣变宽，过渡带增大）
3. 我们对窗函数款式的选择，实质上就是对过渡带宽度和阻带衰减两参数间的权衡，二者不可兼得。在此先放出关键参数，后文再具体分析

• 汉宁窗

 1.通/阻带波动量δp/δs约0.0064，过渡带宽度6.2π/N

    2.主瓣宽度：8π/N

• 海明窗

    1.主瓣宽度与汉宁窗一致，过渡带宽度更大，为7π/N

       2.较汉宁窗，阻带衰减更大

• 布莱克曼窗

1.阻带衰减、过渡带宽度进一步增大

例题2：

附注：1.布莱克曼窗参数：主瓣宽度：12π/N  过渡带宽度：11.4π/N  δp/δs：0.0002  Ap：0.0017dB  As：74dB

2.我们设计的窗函数参数：主瓣宽度<=需求  过渡带宽度<=需求  Ap<=需求  δp/δs<=需求  As>=需求

   3.Ωs-Ωp=过渡带宽度 

   4.Ωs-Ωp=Ωc（将Hd（Ω）IDTFT转为hd[k]时，Ωc作为积分上下限）（FIR的通带阻带截止频率不好控制，为了同时兼顾两方，我们计算时折中，选取两截止频率的中间频率）

本道例题与上道例题的区别：本道题在上道题的基础上，明确地提出了Ap,As，过渡带宽度等指标，因此为了满足这些指标，必须要对窗函数做出筛选；上一题对这些参数没有要求时，默认用矩形窗即可

例题3：

5.第三种类型：可调窗

I。(β)--修正贝塞尔函数，此处不展开

M=N-1，即为滤波器阶数

例子:4：

• 频率取样法设计FIR滤波器（选看）

    窗函数法的操作逻辑，是制作出的可物理实现滤波器时域响应序列尽可能的逼近，模拟理想滤波器时域响应序列(在一定区间的时域序列保持相似或相同)；

频率取样法则是制作出的可物理实现滤波器频域响应幅频特性尽可能的逼近，模拟理想滤波器频域响应幅频特性。（在幅频特性图若干坐标处取值保持相同，类似频域采样）

设计原理：

  如上图，要使用频率取样法设计时域序列长度为N=M+1的序列，阶数为M的FIR滤波器，则必须保证：频域上（0-2π）均匀划分成N等分，在每一份的坐标点上频域取值与理想滤波器保持一致。

    见上图，和前文的方法一样，理想滤波器频谱Hd(e^jΩ)依旧是由用来模拟的FIR滤波器频域幅度函数Ad（Ω）与e^-j*FIR滤波器相位函数相乘来表示

    设计步骤：

这里有一点与之前的设计不同，我们之前滤波器的响应频域范围为-π到π，而在频率取样法中为0到2π。频率取样法中为0到π对应窗函数法的-π到0，频率取样法中π到2π对应窗函数法0到π。因此我们仅需要采样0到π内的点，π到2π的点取值按照所选滤波器的对称特性，将0到π内的点取值映射到π到2π的点即可。

例题5：

注：增大滤波器的阶数，相当于增加在频域采样到的点，因此得到的FIR滤波器频谱会更加接近理想滤波器的样子，但是FIR滤波器的阻带衰耗是不会随滤波器阶数增大而增大。要想增加阻带衰耗，需采用下文的方法。

方法改进：（知道怎么操作就行）

    五、优化设计（选看）

根据一定的误差准则，计算FIR滤波器的系数h[k]，使得误差准则定义的误差ε=[w(Ω)*(A(Ω)-D(Ω))]在0<Ω<π达到最小。