介紹一下FIR濾波器的使用

1. 濾波的目的

在信號(hào)的實(shí)時(shí)處理中，濾波是十分必要的，因?yàn)樾盘?hào)中難免會(huì)由于各種原因混入噪聲，干擾我們對(duì)信號(hào)進(jìn)行分析。

這里強(qiáng)調(diào)了“實(shí)時(shí)”，是因?yàn)樵谝恍﹫?chǎng)景中，我們可能需要對(duì)信號(hào)進(jìn)行離線的、交互式的分析，此時(shí)僅僅是對(duì)原始信號(hào)進(jìn)行采集即可，無需軟件 進(jìn)行濾波，只需要硬件在采集信號(hào)時(shí)做一下抗混疊濾波即可，這種場(chǎng)景不是本文所關(guān)心的。

在更多的場(chǎng)景下，比如我們的運(yùn)動(dòng)手表、手環(huán)在處理ppg信號(hào)計(jì)算心率時(shí)，是會(huì)進(jìn)行實(shí)時(shí)計(jì)算的，那么就會(huì)進(jìn)行實(shí)時(shí)濾波處理：因?yàn)榇蠖鄶?shù)人的心率頻率在0.7~3.6Hz范圍內(nèi)，至少在這之外的頻率會(huì)在計(jì)算心率之前濾除掉。

2. 信號(hào)模擬

比如我們要處理一個(gè)信號(hào)，但是我們僅僅關(guān)心信號(hào)100Hz以下的頻段，這時(shí)我們就需要一個(gè)低通濾波器了，此時(shí)我們先模擬出一個(gè)包含5Hz和3000Hz頻率成分的信號(hào)，假設(shè)信號(hào)采樣頻率為8192，采樣時(shí)間為1秒，共計(jì)8192個(gè)點(diǎn)。信號(hào)生成和展示的代碼如下：

import numpy as np
from numpy import cos
import matplotlib.pyplot as plt

pi = np.pi
t = np.linspace(0, 1, 8192)

signal = 3 * cos(2 * pi * 5 * t + pi/3) + 7 * cos(2 * np.pi * 3000 * t + 3/8*pi)

plt.figure()
plt.plot(t, signal)
plt.show()

我們生成了這樣的一個(gè)信號(hào)：

生成的信號(hào)

3. FIR濾波器系數(shù)生成

這一步可以使用matlab進(jìn)行輔助，本文僅僅是想要一個(gè)截止頻率為10Hz的FIR低通濾波器，步驟如下：

1. 打開matlab;
2. 點(diǎn)擊"APP";
3. 找到濾波器設(shè)計(jì)工具，并點(diǎn)擊；
4. 選擇響應(yīng)類型、設(shè)計(jì)方法、階數(shù)等。

輔助設(shè)計(jì)界面

5. 點(diǎn)擊“文件” 、“導(dǎo)出”、“系數(shù)文件”，導(dǎo)出系數(shù)文件，我把導(dǎo)出的系數(shù)（FIR低通濾波系數(shù)僅有分子）畫出來后如下圖：

系數(shù)波形圖

3. FIR濾波原理

采用FIR進(jìn)行濾波，從操作上看是進(jìn)行卷積操作，對(duì)上述濾波器的系數(shù)進(jìn)行FFT變換即可窺見一斑：

import numpy as np
import matplotlib.pyplot as plt
# b = [......] 101個(gè)系數(shù)組成的列表，此處省略
delta_f = 1
plt.plot([delta_f * i for i in range(4097)], abs(np.fft.rfft(b, 8192)))  # 單邊FFT
plt.show()

系數(shù)的傅里葉變換

這里得到的就是濾波器的幅頻響應(yīng)曲線，和濾波器輔助設(shè)計(jì)工具中所展示的幅頻響應(yīng)曲線是一致的。

4. 濾波計(jì)算代碼與結(jié)果

把第二步生成的信號(hào)中高于100Hz的頻率成分（即3000Hz的成分）濾除，得到的結(jié)果如下圖所示，在結(jié)果的開頭和結(jié)尾產(chǎn)生了失真，這是因?yàn)樵诰矸e運(yùn)算過程中，在剛開始和結(jié)尾計(jì)算時(shí)有效信息不足進(jìn)行了補(bǔ)零操作導(dǎo)致的。

濾波的結(jié)果

這一步的運(yùn)算代碼如下：

import numpy as np
from numpy import cos
import matplotlib.pyplot as plt
from filter_coeff import b
pi = np.pi

t = np.linspace(0, 1, 8192)
signal = 3 * cos(2 * pi * 5 * t + pi/3) + 7 * cos(2 * np.pi * 3000 * t + 3/8*pi)

fir_len = len(b)

res = []
for i in range(len(signal)):
    print(i)
    temp = 0
    if i < fir_len:
        data = (fir_len-i) * [0] + signal[0:i].tolist()
    elif i > len(signal) - fir_len:
        data = signal[i:].tolist() + (fir_len - len(signal) + i) * [0]
    else:
        data = signal[i:i+fir_len]
    for j in range(fir_len):
        temp += b[j] * data[-j-1]
    res.append(temp)


plt.figure()
plt.plot(t, res)
plt.show()