数字信号处理5：FIR滤波器设计

Posted 2022-11-22 zuguorui

文章目录

• 1. 滤波器初识
• 2. 最直观的滤波方式：频域滤波
• 3. 傅里叶变换中的加窗
• 4. FIR滤波器设计
• 5. 总结

之前的一系列博客中，详细分解了从卷积到FFT的相关知识，不过那些属于理论，是为了让我们清楚认识到信号处理的本质。本篇博客将会详细讲解数字信号处理最广泛的应用——滤波器。

注意，本章所采用的dB标准为

$$dB=20log\frac{V}{V_{base}}$$

其中，$V_{base}$表示某种作为标准的幅度值，它或者是电平，或者是信号振幅。在这里，我们通常将它认为是原信号的“大小”，这个大小不光可以是信号在时域的幅度，也可以是在频域的频率分量的大小。

1. 滤波器初识

滤波器，如同字面意思是过滤信号的某些频率分量。入门时书上最常说的几种滤波器，比如低通、带通、高通、带阻等滤波器。如下是一个低通滤波器。

它将信号中的中高频分量降低了60dB。按照dB的定义，也就是中高频信号的幅度降低为原来的0.001倍。

实际上，更为广泛的滤波器，并不只有“滤”的作用，它同样还可以对某些频段进行增强，如下滤波器：

它将信号中的中高频部分提高了60dB。按照dB的定义，也就是中高频信号的幅度提高到原来的1000倍。

2. 最直观的滤波方式：频域滤波

回想一下《数字信号处理2：傅里叶变换》博客中讲的傅里叶变换的直观理解那部分。傅里叶变换之后的结果是一组虚数，这组虚数的相位反应的就是原信号中各频率正弦波的相位，而这些虚数的绝对值某种程度上反应了对应频率的信号在原始信号中的大小。所以，我们只要对频域上的这些虚数的绝对值进行调节，再作傅里叶逆变换，即可以调整原信号中的某些频率分量的大小，也就成为了滤波的一种手段。步骤通常如下：

1. 对信号加窗
2. 对加窗之后的结果作傅里叶变换。
3. 根据你的目的对傅里叶变换结果的某些频率上的分量调节绝对值。
4. 将调节完毕后的结果作傅里叶逆变换。

这便是频域滤波。

3. 傅里叶变换中的加窗

在上一部分中提到了对信号进行加窗，为什么要进行加窗操作？加的又是什么窗？

要理解加窗的意义，有两种方式。

首先，回想一下博客《数字信号处理2：傅里叶变换》中对傅里叶变换的推导过程，对于傅里叶变换所做的所有推导，都是基于无限长周期信号进行的。我们假设用来做傅里叶变换的信号是从周期信号中截取的一整个周期。所以，如果这段信号前后进行拼接成我们假想的周期信号后，中间没有突变点（包括斜率突变、断点等情况），那么直接进行傅里叶变换是没有问题的。而如果前后进行拼接后，中间出现了突变点，那么从直觉就可以知道，这个突变点中包含了相当多的频率分量。但单从我们截取的这段信号来说，是不应该有这些多余的频率分量的。

下面举个例子：

构造一个频率为20的周期信号。首先，我们让这段信号有整数个周期。使用matlab。

T = 200;

x = [0:1:T-1];
sinWave = sin(2 * pi / T * 20 * x);
plot(x, sinWave);
title("频率为20的正弦波");
b = fft(sinWave);
plot(x, abs(b));
title("整数个周期正弦波的频谱")

代码中生成的sin波是完整的20个周期，采样点200个，周期是10。绘制出来的图像如下：

可以看到，频谱上仅在周期为20的位置有值。

如果不给整数个周期。比如只给19.5个周期。显然，这样的信号前后拼接成周期信号后，必然是会有突变点的。

x = [0:1:195];
sinWave = sin(2 * pi / T * 20 * x);
plot(x, sinWave);
title("19.5个周期的正弦波");
b = fft(sinWave);
%b = circshift(b', fix(length(sinWave) / 2))';
plot(x, abs(b));
title("非整数个周期正弦波的频谱")

可以看到，这时在频率20周围出现了其他的频率分量，这就是频率泄漏。

实际运用中，一段信号通常有很多频率的正弦波组成，我们不能保证总是能截取到整数个周期，因此，通过加窗的方式，我们使信号两端加一个渐缓的过程，使得拼接成周期信号后，没有突变点。

举一个例子：汉宁窗。它的形状如下：

对于上面提到的非整数个周期的信号，我们对其进行加窗操作：

x = [0:1:195];
window = hanning(length(x))';
plot(window);
title("汉宁窗");
sinWave = sin(2 * pi / T * 20 * x) .* window;
plot(x, sinWave);
title("加窗后的19.5个周期的正弦波");
b = fft(sinWave);
%b = circshift(b', fix(length(sinWave) / 2))';
plot(x, abs(b));
title("加窗后的非整数个周期正弦波的频谱")

之后信号就变成这样

其频谱

可以看到，加窗之后，频率泄漏的情况改善了很多。

从另一方面理解。之前的博客《数字信号处理4：采样定理》中证明了一个重要的结论：时域卷积等于频域相乘，时域相乘等于频域卷积。

因此，将信号做傅里叶变换，其实隐含了一个过程，就是原信号的频谱与窗函数的频谱进行卷积。而如果不想改变原信号频谱，最理想的情况就是让它与一个冲激信号进行卷积，这样原信号的频谱其实只是发生了位移。这样一来，我们就希望能有一个窗函数，它的频谱是一个冲激函数，或者尽量接近冲激函数。

对于一段信号，不加窗等价于加矩形窗。我们不妨看一下矩形窗的频谱：

a=zeros(1, 400);
a(100:299) = ones(1, 200);
b = fft(a);
b = circshift(b', 200)';
plot(a);
title("矩形窗")

plot([-199:1:200], abs(b));
title("矩形窗频谱")

显然，矩形窗有很多旁瓣，与这样的信号在频域上做卷积，会改变原信号的频谱。

而看一下汉宁窗：

a=zeros(1, 400);
a(100:299) = hanning(200)';

plot(a);
title("汉宁窗")
b = fft(a);
b = circshift(b', 200)';

plot([-199:1:200],abs(b));
title("汉宁窗频谱")

可见汉宁窗的频谱相当接近一个冲激函数。

至此，窗函数的意义已经明了：消除傅里叶变换过程中的频率泄漏。

4. FIR滤波器设计

滤波器是一个线性时不变系统，我们使用一个差分方程来表示该系统的信号传输特性：

$$y[n]=-\sum _{k=1}^N{a}_{k}y[n-k]+\sum _{k=0}^{M-1}{b}_{k}x[n-k]$$

如果$a_k\ne 0$，那么当前输出不仅与当前输入有关，还有过去的$M-1$个输入，以及过去的$N$个输出有关。这样的滤波器成为无限冲激响应滤波器（IIR Filter）。如何理解无限这个词，当前输出与过去的输出有关，那么可以预见，第一个输入对后面所有的输出都会有影响。

如果$a_k=0$，那么当前输出仅与当前输入和过去的$M-1$个输入有关。这样的滤波器被成为有限冲激响应滤波器（FIR Filter）。

这两种滤波器各有优缺点，IIR的优点在于计算快，能以较少的阶数达到性能要求。但设计比较复杂，而且难以设计出具有准确频率响应的滤波器，另外，IIR滤波器的相位不可能是线性的。

而FIR滤波器则反了过来，设计简单，能快速设计出具有精确线性相位以及需要的频率响应的滤波器。但是需要较高的阶数才能达到滤波要求。

设$h[n]$是滤波器，$H[k]$是$h[n]$的傅里叶变换。

我们从频域滤波开始，既然从频域调节就可以达到滤波器的目的，那么干脆将频域上调节的系数看做$H[k]$，然后用傅里叶逆变换求出$h[n]$，频域相乘等于时域卷积，使用$h[n]$与信号$x[n]$进行卷积即可。

首先要明确的一个点，为什么滤波器要是线性相位？按照傅里叶变换公式

$$X[k]=|X[k]|e^{-j{a}_k{w}_0}$$

然后有一滤波器$h[n]$，其频谱是

$$H[k]=|H[k]|e^{-j{b}_k{w}_0}$$

当$X[k]$与$H[k]$的每一项相乘。
$$Y[k]=H[k]X[k]=|H[k]||X[k]|e^{-j(a_k+b_k)w_0}$$

我们已经知道傅里叶变换是可逆的。所以对于$X[k]$来说

$$x[n]=\frac{1}{N}\sum _{k=<N>}|X[k]|e^{-j{a}_k{w}_0}{e}^{jk{w}_{0}n}=\frac{1}{N}\sum _{k=<N>}|X[k]|e^{j(kn-a_k)w_0}$$
得出的结果必然是实数。
那么对于$Y[k]$：

$$y[n]=\frac{1}{N}$$