Python 中 FIR 数字滤波器设计与时序信号处理——信号与系统大作业

Author Avatar
akimoto-cris 12月 26, 2017

概要

FIR 即有限脉冲响应(finite impulse response)，它对应着无限脉冲响应。一般在现实生活中，由于数字方法处理信号时，不可能对无限长的信号进行测量和运算，而是取其有限的时间片段进行分析。能够做到的是用窗函数截断出近似的信号进行处理，或者对信号时间片段进行周期延拓等等方法。

FIR滤波器的特点：

FIR 滤波器可以得到严格的线性相位,但它的传递函数的极点固定在原点,只能通过改变零点位置来改变性能

FIR 滤波器的设计方法主要有三种：窗函数法、频率抽样法和切比雪夫逼近方法。窗函数法较为简单易行，本次作业就窗函数设计的角度展开。

窗函数的“窗”指的是频谱上窗型的频率响应，在频域上有限宽度的窗，经过离散傅里叶反变换后，在时域上就是无限长的序列。不同形状的窗函数的滤波性能也不同，不仅取决于级数，还与因果与非因果等因素有关。

作业中的代码和演示在 Jupyter Notebook 中完成。

In this project, with some previous Python programming experiences, I basically learned:

  1. How to complete basic tasks of data processing with fundamental python packages. (numpy, scipy, etc.)
  2. How to design digital FIR filters with different module of FIRs.
  3. The difference between 2 methods of signal filtering encapsuled within scipy.signal package.
  4. How to read, extract and process data from wfdb-python package, which is an open source ECG dataset.

本编

首先调用所需的 Python 库：

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
import wfdb
from matplotlib.font_manager import FontProperties

myfont = FontProperties(fname=r'C:/Windows/Fonts/msyh.ttf')
# set default size of plots
plt.rcParams['figure.figsize'] = (15, 4)

理想窗函数滤波器

理想窗函数的频响如图，频谱在通带和阻带之间直接截断，没有过渡带，如果是一个非因果的双边滤波器，则他的时域响应是一个离散的升余弦函数：
$$h_{ideal}[n] = \frac{\sin(2 \pi f_c)}{\pi n} = 2f_c sinc(2f_c n)$$

而一个因果的单边理想窗函数 LPF 的频响示意图如下：

接下来用理论的时域响应函数通过FFT得到因果和非因果窗函数各自的频响：

#one side signal (karma)
def sinc_n(M, fc):
    return 2 * fc * np.sinc(2 * fc * np.arange(0, M, 1.0))

#two side signal (non-karma)
def sincdouble_n(M, fc):
    return 2 * fc * np.sinc(2 * fc * np.arange(-M, M, 1.0))

def sinc_freqz(b, title, semilogx = False, showpart = False):
    ''' 
    calculate Z transform.
    Cut-out in time domain signal cause side lobes in Frequency Responce'''

    w, h = signal.freqz(b)
    h_dB = 20 * np.log10(abs(h))
    f = w / (np.pi)       #regularize frequncy to be between 0 and 1

    plt.rcParams['figure.figsize'] = (15, 4) 
    if showpart == False: plt.subplot(122)
    # select plot mode: log space or linear space
    if semilogx == True:
        plt.semilogx(f, h_dB)
    else:
        plt.plot(f, h_dB)
        plt.plot([0, 1], [-21, -21], 'k--', LineWidth = 1.0)
        plt.text(0.26, -21, 'MSB of sinc double = -21 dB', fontdict={'size': 16})
    plt.xlabel('$\omega$')
    plt.ylabel('$H (dB)$')

    # if compare between different sample number
    if showpart == False:
        plt.subplot(121)
        plt.stem(b)
        plt.plot(b, '--')
        plt.xlabel('$n$')
        plt.ylabel('$h(n)$')
        plt.title('FIR ' + title)
        plt.show()

    return f, h_dB

# M is width of window in time domain
M = 30
fc = 0.125
# one/two side sinc signal in time domain
b = sinc_n(M, fc)
b1 = sincdouble_n(M, fc)

f, h_dB = sinc_freqz(b, 'rectangle window')
f1, h1_dB = sinc_freqz(b1, 'rectangle window (two-side)')

结论：

  1. 将两种理想窗函数相比较后，两种窗函数都出现了吉布斯现象（由于 FIR 本身特性不可避免）两种理想右上角的子图中，低通旁瓣更加明显；单边窗出现了负谱现象，并且高频段的抑制效果比双边窗的差很多。
  2. 比较单边和双边窗函数的 MSB （最小阻带衰减），单边窗的 MSB 远小于双边窗的，这意味着单边窗的低通滤波性能远不如双边窗。而双边窗的 MSB 增益也只有 -21dB ，这也是不足以实际产生功用的。

进一步探究窗的宽窄对 MSB 参数的影响：

Ms = [8, 16, 32, 64, 90, 160]
cache = dict()
for i, M in enumerate(Ms):
    cache['b' + str(i)] = sincdouble_n(M, fc)
    _, _ = sinc_freqz(cache['b' + str(i)], 'rectangle window, M = '+str(M), showpart=True)
plt.xlim([0.2, 1])
plt.show()

结论：
改变窗的宽窄对理想双边窗函数的 MSB 参数没有影响，MSB 之和零点的位置和多少有关。

其他常用的窗函数

矩形窗，也就是理想窗由于过于突出的缺点——旁瓣较高即 MSB 小导致变换中带入高频干扰和泄露，在实践中往往被其他窗函数代替，以实现更快的衰减，更小的泄露等关键性能。

常见窗函数：

  • Hanning 汉宁窗：
    $$\omega[k] = 0.5 - 0.5\cos \frac{2\pi k}{M}\quad(0 \leq k \leq M)$$
    优点：与矩形窗相比，泄漏、波动都减小,并且提高选择性
  • Hamming 海明窗：
    $$\omega[k] = 0.54 - 0.46\cos \frac{2\pi k}{M}\quad(0 \leq k \leq M)$$
    特点：仅与 Hanning 的权重不同，都是改进的升余弦窗
  • Gaussian 窗：课堂上经常提到，不展开
    -Blackman 布莱克曼窗：
    $$\omega[k] = 0.42 - 0.5\cos \frac{2\pi k}{M}+0.08\cos \frac{4\pi k}{M}\quad(0 \leq k \leq M)$$
    特点：二阶升余弦窗，主瓣宽，旁瓣比较低，但等效噪声带宽比汉宁窗要大一点，波动却小一点。频率识别精度最低，但幅值识别精度最高，有更好的选择性。
  • Kaiser 凯泽窗：
    $$\omega[k] = \frac{I_0[\beta \sqrt{1 - (1 - \frac{2k}{M - 1})^2}]}{I_0 [\beta]}\quad(0 \leq k \leq M - 1)$$
    特性：Kaiser 窗的性能和旁瓣宽度与参数$\beta$有关

利用 scipy 科学计算包中封装好的 signal.firwin 类方法来可视化这些窗函数：

def get_window(wintype, Nx, cutoff = 0.5):   
    ''' Nx is equivalent to M '''

    wins = dict()

    for w in wintype:
        win = signal.firwin(2 * Nx, cutoff, window=w)
        wins[w] = win

        plt.rcParams['figure.figsize'] = (15, 10) 
        plt.subplot(211)
        plt.plot(win, label = w)
        plt.xlabel('$n$')
        plt.ylabel('$w(n)$')
        plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)

        plt.subplot(212)

        win_f, h = signal.freqz(win)
        plt.plot((win_f / (2 * np.pi))[200:], (20 * np.log10(abs(h)))[200:], label = w)
        plt.xlabel('$\omega$')
        plt.ylabel('$H (dB)$')
        plt.title('amplitude response')
        plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)

    return wins

firwins = get_window(["blackman", 
            "hamming", 
            "hann", 
            ('kaiser', 4.5513)], 61)        #beta = 4.5513
plt.show()

结论：

  1. 这些常见窗函数的 MSB 特性都优于理想矩形窗函数（无论是双边还是单边），
  2. 由于是二阶滤波器，给定窗口范围时的 Blackman 窗的旁瓣增益最大；Kaiser 窗的 MSB、旁瓣宽度等形状特点根据 $\beta$ 的改变而不同。

将 Blackman 窗的频响与双边矩形窗对比：

b3 = signal.firwin(len(b1), 0.22)

w3, h3 = signal.freqz(b3)

h3_dB = 20 * np.log10(abs(h3))
f3 = w3 / (np.pi)

plt.rcParams['figure.figsize'] = (15, 4) 
plt.plot(f1, h1_dB, label = 'rect')
plt.plot(f3, h3_dB, label = 'firwin')
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.show()

可以看出前者的性能明显优于后者。

signal.filtfilt 和 signal.lfilter

当我们对上述设计并实现的 FIR 数字滤波器进行比对和选择之后，遴选出的一些滤波器将会用于实际应用。这就需要用到 scipy 库中给出的信号与滤波器相作用的方法，其中有两种不同的滤波函数：signal.filtfilt 和 signal.lfilter。在查阅官方 doc 之后，由于两个 filter 函数的 API 都很类似，没有很好地展示两个实际上的区别，

下面用一个简单的信号来 demo 一下两者的区别，输入是一个低频率和一个高频率叠加形成的信号：

def filtfilt_demo(b,a,x, lim = (0, 400)):
    (l, h) = lim
    y = signal.filtfilt(b, a, x)
    plt.plot(y[l: h])
    plt.show()
    return y

def lfilter_demo(b,a,x, lim = (0, 400)):
    (l, h) = lim
    y = signal.lfilter(b, a, x)
    plt.plot(y[l: h])
    plt.show()
    return y

t = np.linspace(0, 1.0, 2001)
'''
contruct signal with a low freq component and a high freq component
'''
x_high = np.sin(2 * np.pi * 5 * t)
x_low = np.sin(2 * np.pi * 250 * t)
x = x_high + x_low

plt.rcParams['figure.figsize'] = (10,3)
plt.plot(x[0:250])
plt.show()

#build Butterworth filter
nom, denom = signal.butter(8, 0.125, btype='low')    # LPF
_ = filtfilt_demo(nom, denom, x)
_ =lfilter_demo(nom,denom, x)

结论：
对于低通滤波器，用 signal.filtfilt 施加在输入信号上时，输出信号的前几个时刻就已经很好的显现出低频信号的效果，这是因为函数将滤波器一前一后作用在了输入上两次；
而 signal.lfilter 只对输入信号作用一次，因此需要迟疑一段时间才能收敛到理论输出信号上。

对 ECG 时序信号的处理

ECG，即 Electrocardiogram 心电图，是记录人的心脏跳动等特征的一种常见医疗影像数据，由心电仪探测得到的时序 ECG 数据不可避免的会夹杂环境噪声，下面利用开源数据集 wfdb 进行心电图时序信号的处理

首先从实现下载并部署到 python 库中的 wfdb dataset 提取出一列 ECG 采样数据 sampledata/100，采样点设置为 15000 个:

signals, fields=wfdb.srdsamp('sampledata/100', sampfrom=100, sampto=15000)

# Convert one channel to processable data
# Argument "fs_target" stands for resampling freq
data, _ = wfdb.processing.resample_sig(x=signals[:,0], fs=fields['fs'], fs_target=1280)

plt.plot(data[:5000])
plt.show()

打印出其中一段提取出的数据：

可以看到这段记录的心电图曲线叠加了很多细微的高频噪声，但尚未淹没基本的上升下降的规律——每个周期内都有一个很短的高峰值紧跟着两到三个小幅振动。我们只想要保留这些相对噪声而言比较稳定的低频信号，而将噪声过滤掉。

通过傅里叶变换得到 wfdb data 的频谱：

'''select a test freqency to observe the spectrum of the signal
'''
f_test = 0.3

#add a f_test frequency sin wave to the signal
data_prime = data +  0.1 * np.sin(f_test * np.arange(0, len(data), 1.0))

spectrum = np.fft.fft(data_prime)
spectrum_test = np.fft.fft(data)
plt.plot(abs(spectrum), label = 'signal')
plt.plot(abs(spectrum_test), label = 'test')
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.show()

下一步是用 FIR 低通滤波器过滤提取出的 ECG 数据。滤波器使用 Blackman 窗。

firwin = get_window(['blackman', ('kaiser', 6)], 500, 0.2)

d_firwin = np.convolve(firwin['blackman'], data)
d_b3 = np.convolve(b3, data)
lim = 5000
plt.rcParams['figure.figsize'] = (15, 5) 
plt.subplot(311)
plt.plot(data[: lim])
plt.ylim([-1.0, 1.0])
plt.text(lim, 0, 'original')
plt.subplot(312)
plt.plot(d_firwin[: lim])
plt.ylim([-1.0, 1.0])
plt.text(lim, 0, 'firwin')
plt.subplot(313)
plt.plot(d_b3[: lim])
plt.ylim([-1.0, 1.0])
plt.text(lim, 0, 'b3')
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.ylim([-1.0, 1.0])
plt.show()

由于 窗函数的阶数较低，效果不是很理想。下面使用 signal.butter 函数提供的 Butterworth 型滤波器，现在将多阶滤波器的滤波效果进行对比：

nom1, denom1 = signal.butter(8, 0.038, btype='low')
d_filtered = filtfilt_demo(nom1, denom1, data, (0, 2000))

滤波后的输出信号如下图

分析 Butterwoth 的零极点情况：

(r, p, k) = signal.residuez(nom1, denom1)

plt.plot(r.real, r.imag, 'o')
plt.plot(p.real, p.imag, 'x')
circ = plt.Circle((0, 0), radius=1, edgecolor = 'b', facecolor = 'None')
plt.gca().add_patch(circ)
plt.xlim([-2, 2])
plt.ylim([-2, 2])
plt.gca().set_aspect('equal')
plt.show()

零点和极点的个数各有 8 个，此时可以达到良好的滤波效果。

将输入信号与输出信号进行波形的对比（为了可视化效果，在同一张图上显示时将两个波形错开）：

plt.plot(data[:2000], label = 'signal')
plt.plot(d_filtered[:2000] - 0.1, label = 'Butter')
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.show()

可以看出通过 Butterworth 滤波器的输出波形基本还原了输入心电图数据的走势，很好地分离了高频低频分量，同时也保留了脉冲峰值高度（损失了大约 0.2 的增益）

如果将刚才的阶数为 8 的滤波器改成高通滤波器，输出应该会得到只含有高频分量的信号。将参数 btype 设为 "high"：

nom2, denom2 = signal.butter(8, 0.038, btype='high')
d2_filtered = filtfilt_demo(nom2, denom2, data, (0, 2000))

plt.plot(data[:2000], label = 'signal')
plt.plot(d_filtered[:2000] - 0.1, label = 'low')
plt.plot(d2_filtered[:2000], label = 'high')
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.show()

观察波形，发现和 HPF 的输出和原始波形中的大部分噪声的震动部分相同，如果想尝试进一步优化 LPF 输出波形，一个自然的想法是将这两部分做减法，观察将得到什么结果：

k = int(len(data) / 256)
r = len(data) % 256

for i in range(1, k):
    d2_filtered[i * 256 + 1: (i + 1) * 256] = data[0: 255] - np.mean(data[0: 255])
d2_filtered[len(data) - r: len(data) - 1] = data[255 - r + 1: 255] - np.mean(data[0: 255])

mini = min(d2_filtered[0: 750])
maxi = max(d2_filtered[0: 750])

for i in range(len(data)):
    if d2_filtered[i] > maxi: d2_filtered[i] = maxi
    elif d2_filtered[i] < mini: d2_filtered[i] = mini

plt.plot(d2_filtered[:2000], label = 'high')
plt.plot(data[: 2000] - d2_filtered[:2000], label = 'd_filtered_new')
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.show()

结果显示这个尝试并不成功。黄色的曲线显示了做减法之后的输出波形，尽管原有的脉冲峰值基本没有损失，但曲线并没有变得光滑。

链接

源代码 / 英文报告：https://github.com/Akimoto-Cris/SignalSystem_project

参考：
http://blog.csdn.net/shenziheng1/article/details/53366067

wfdb dataset: https://github.com/MIT-LCP/wfdb-python

Powered by Hexo
Theme - Material