はじめに

双2次フィルタ (SOS) やバターワース／チェビシェフフィルタなどといったIIRフィルタは，FIRフィルタと比べると少ない計算量で急峻な周波数特性を実現できる便利なフィルタです。

しかし，基本的に位相特性は非線形であり直線位相ではありません。出力波形が大きく崩れてしまう問題があります。

そこで、位相の遅れが発生しないフィルタリング方法としてゼロ位相フィルタリングという手法があります。

この記事では、ゼロ位相フィルタの実装方法と、scipy.signal.filtfilt の利用方法を紹介します。

ゼロ位相フィルタリングの実装方法

ゼロ位相フィルタリングの方法は簡単です。

フィルタリングしようとしている IIRフィルタ係数を b, a とすると、

IIRフィルタを適用する y = iir(b, a, x)
1 の信号を反転する y_inv = y[::-1]
2 の信号に1と同じIIRフィルタを適用する z = iir(b, a, y_inv)
3 を反転する z_inv = z[::-1]

最終的な出力は z_inv です．

すなわち、入力信号に対してIIRフィルタを順方向と逆方向の両方から適用することで位相遅れをキャンセルしています．

scipy.signal.filtfiltを使ったゼロ位相フィルタ

Pythonではscipy.signal.filtfiltを使うことでゼロ位相フィルタリングを簡単に行うことができます．

import matplotlib.pyplot as plt
import numpy as np
import scipy.signal

duration = 0.1
sr = 16000
t  = np.arange(0, sr*duration) / sr
fo = 440
x  = scipy.signal.sawtooth(2*np.pi*fo*t)
x  = np.r_[np.zeros(100), x]
 
# サンプリング周波数の 1/4 の周波数以下を通すバターワースローパスフィルタ
b, a = scipy.signal.butter(8, 0.25, btype='low') 
 
# 通常のIIRフィルタリング
y_lfilter  = scipy.signal.lfilter(b, a, x)
# ゼロ位相フィルタリング
y_filtfilt = scipy.signal.filtfilt(b, a, x)

plt.plot(x, label="original", alpha=0.5)
plt.plot(y_lfilter, label="lfilter", alpha=0.5)
plt.plot(y_filtfilt, label="filtfilt", alpha=0.5)
plt.legend()
plt.show()

通常のIIRフィルタリングであるscipy.signal.lfilterは入力信号と比べると全体の遅延があり、のこぎり波の形状が大きく変形しています．

一方でゼロ位相フィルタ（scipy.signal.filtfilt）では入力信号と比べると重なっていることから全体の遅延がない（群遅延時間ゼロ）が確認できます。また、波形の変形もlfilterよりは抑えられているように見えます．

ただし、ゼロ位相フィルタのデメリットもあります．

逆方向からのフィルタリングを行うため，信号の立ち上がり前に振動してしまう（プリエコー／プリリンギング）
未来の信号に影響するため、リアルタイム処理に不向き*1
2回IIRフィルタリングする必要があるため計算量が増える

なお、振幅応答は元のIIRフィルターの2乗となります．

ylim_min_dB = -120
to_dB = lambda h : 20 * np.log10(np.abs(h)/np.max(np.abs(h)))

w, h_lfilter = scipy.signal.freqz(b, a)

amp_dB_lfilter  = to_dB(h_lfilter)
amp_dB_filtfilt = to_dB(np.abs(h_lfilter)**2)

# 振幅応答のプロット
plt.plot(w, amp_dB_lfilter, alpha=0.7, label="lfilter")
plt.plot(w, amp_dB_filtfilt, alpha=0.7, label="filtfilt")

plt.ylabel('Amplitude [dB]')
plt.xlabel('Frequency [rad/sample]')
plt.ylim(ylim_min_dB, 0)

plt.grid()
plt.legend()
plt.show()

lfilter x 2 でのゼロ位相フィルタリング

ゼロ位相フィルタの実装方法より、lfilterを2回かければよいことが分かりますのでfiltfiltの結果と一致するか確認しました。

z_lfilter = scipy.signal.lfilter(b, a, y_lfilter[::-1])[::-1]

plt.clf()
plt.plot(y_lfilter, label="lfilter", alpha=0.5)
plt.plot(z_lfilter, label="lfilter-lfilter", alpha=0.5)
plt.plot(y_filtfilt, label="filtfilt", alpha=0.5)
plt.legend()
plt.show()