IIRフィルタ（Infinite Impulse Response）の原理
前回「音を分離する（その１：FIR）」に引き続き、音データから特定の周波数を残したり、取り除いたりする方法の一つとして、IIRフィルタについて学びます。
IIRフィルタは以下の式によって定義されます。
\( s_1(n) = \sum_{m=0}^I b(m)s_0(n-m) -\sum_{m=1}^J a(m)s_1(n-m) \)
- \( s_1 \)：フィルタ適用後の音データ
- \( s_0 \)：フィルタ適用前の音データ
- \( b \)：フィルタ係数
- \( a \)：フィルタ係数（フィードバック項）
- \( m \)：フィルタ係数列の添え字
- \( n \)：音データ列の添え字
IIRフィルタの特徴はフィードバック項が存在することです。
フィードバック項とは式中の「\( -\sum_{m=1}^J a(m)s_1(n-m) \)」にあたる箇所のことで、IIRフィルタの式はFIRフィルタの式にちょうどこのフィードバック項を加えた形をしています。
IIRフィルタの処理をブロック図にしたものを見てください。
ブロック図の左方はFIRフィルタとまったく同じ構造しており、右方はフィードバック項に該当します。
このフィードバック項は過去に一度計算して出力した\( s_1(n-1) \)という値を、次の\( s_1(n) \)の計算に再利用します。過去\( J \)回分の出力を加味して次の値の計算が行われるため、データの時系列特性が結果により強く現れます。
また、入力音\( s_0 \)の末尾に無音データが続いたとしても、フィードバック項によって何かしらの応答（フィルタ結果の値）が返ってきます。
低域通過フィルタ（Low-Pass-Filter）
低域通過フィルタを適用してカットオフ周波数よりも低い周波数の音だけを残します。カットオフ周波数には600Hzを設定します。
※音量には十分に注意してください。
高域通過フィルタ（High-Pass-Filter）
高域通過フィルタを適用してカットオフ周波数よりも高い周波数の音だけを残します。カットオフ周波数は600Hzです。
※音量には十分に注意してください。
帯域通過フィルタ（Band-Pass-Filter）
帯域通過フィルタを適用して特定の周波数帯（バンド）の音だけを残します。
バンドは600Hz～2500Hzを設定します。
※音量には十分に注意してください。
帯域阻止フィルタ（Band-Eliminate-Filter）
帯域阻止フィルタを適用して特定の周波数帯（バンド）の音だけを取り除きます。バンドは600Hz～2500Hzを設定します。
※音量には十分に注意してください。
Python+scipyによる実装
IIRフィルタはscipyのsignalモジュールによって実装可能です。
以下に低域通過フィルタ（LPF）の例を示します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.io import wavfile
OUTPUT_FILENAME = "output.wav"
# パラメータの設定
RATE = 16000 # 標本化周波数（サンプリングレート）
DURATION = 3 # 録音時間[s]
FREQ = 250 # 波形周波数
AMP = 0.3 # 振幅
FRAMES = 500 # 描画する範囲
N = 10 # IIRフィルタの次元数
CUTOFF = 600 / (RATE / 2) # カットオフ周波数（ナイキスト周波数で正規化）
# 時間軸の作成
time = np.arange(RATE * DURATION) / RATE
# 正弦波
sin_wave = AMP * np.sin(2 * np.pi * FREQ * time) + AMP * np.sin(2 * np.pi * 2000 * time) + AMP * np.sin(2 * np.pi * 3750 * time)
# 周波数軸の作成
freq = np.fft.fftfreq(len(time), 1.0 / RATE)
# フーリエ変換
spectrum = np.fft.fft(sin_wave)
# IIR filter
b, a = signal.iirfilter(N, CUTOFF, btype="lowpass")
w, h = signal.freqs(b, a, 1000)
lpf_wave = signal.lfilter(b, a, sin_wave)
lpf_freq = np.fft.fftfreq(len(lpf_wave), 1.0 / RATE)
lpf_fft = np.fft.fft(lpf_wave)
# 描画設定
fig, axs = plt.subplots(2, 1, figsize=(8, 4))
axs[0].plot(time[: FRAMES], sin_wave[: FRAMES], label="raw")
axs[0].plot(time[: FRAMES], lpf_wave[: FRAMES], label="filtered")
axs[1].plot(freq[: len(freq)//2], np.abs(spectrum)[: len(spectrum)//2], label="raw")
axs[1].plot(freq[: len(freq)//2], np.abs(lpf_fft)[: len(lpf_fft)//2], label="filtered")
axs[0].set_title("wave")
axs[0].legend(loc='upper right')
axs[1].set_title("frequency")
axs[1].legend(loc='upper right')
plt.tight_layout()
plt.show()
lpf_wave = lpf_wave * (2**15-1)
lpf_wave = lpf_wave.astype("int16")
wavfile.write(OUTPUT_FILENAME, RATE, lpf_wave)iirfilter関数でフィルタ係数の配列を作成し、合成波にはlfilter関数で適用します。
定数「N」はフィルターの次元数、「CUTOFF」はカットオフ周波数を指定します。iirfilter関数のオプション「btype」はlowpass, highpass, bandpass, bandstopのいずれかを指定します。
コメント