音を分離する（その１：FIR）

Python
2023.05.20

FIRフィルタ（Finite Impulse Response）の原理

今回は音データから特定の周波数を残したり、取り除いたりすることで音データを加工する方法について学びます。
その第１弾として以前投稿した記事「音を響かせる」でも登場したFIRフィルタについて掘り下げていきます。

FIRフィルタは以下の式によって定義されるデジタルフィルタです。
これは「畳み込み」と呼ばれる計算の定義とまったく同じです。

\( s_1(n) = \sum_{m=0}^J b(m)s_0(n-m) \)

  • \( s_1 \)：フィルタ適用後の音データ
  • \( s_0 \)：フィルタ適用前の音データ
  • \( b \)：フィルタ係数
  • \( m \)：フィルタ係数列の添え字
  • \( n \)：音データ列の添え字

下図はこの式を理解しやすくするために、ブロック図にしたものです。

ブロック図は「加算器」「乗算器」「遅延器」の基本三要素と呼ばれる3つのパーツからできています。

  1. 入力された音声\( S0(0) \)は最上段の乗算器と遅延器に渡される。
  2. ２段目の遅延器は音声\( S0(0) \)を記録し、初期値０を２段目の乗算器と遅延器に渡す。
  3. 最下段の遅延器は同様に２段目の遅延器から受け取った０を記録し、初期値０を最下段の乗算器に渡す。
  4. 各乗算器は受け取った値に各フィルタ係数値\( b(n) \)を乗算し、後続の加算器に渡す。
  5. 加算器は受け取った値を加算し、次の加算器または出力に値を渡す。
  6. 以降、1～5の手順を繰り返す。
    ただし、音声\( S0 \)の添え字は１つずつ増加し、遅延器は初期値ではなく直前に記録した値を渡すようになります。

以前は音データにディレイやリバーブといった残響効果かけるために使用しましたが、FIRフィルタは特定の周波数成分を残したり、取り除いたりすることができます。

  • 低域通過フィルタ（Low-Pass-Filter）
    低周波成分をだけを残し、高周波成分を減衰または除去するフィルタです。
  • 高域通過フィルタ（High-Pass-Filter）
    高周波成分をだけを残し、低周波成分を減衰または除去するフィルタです。
  • 帯域通過フィルタ（Band-Pass-Filter）
    特定の周波数帯域（Band）だけを残し、それ以外の成分を減衰または除去するフィルタ
  • 帯域阻止フィルタ（Band-Eliminate-Filter）
    特定の周波数帯域（Band）だけを減衰または除去し、それ以外の成分を残すフィルタ

この記事では250Hz、2000Hz、3750Hzの正弦波の合成波に対して上記の各フィルタを適用していきましょう。

※高い周波数の音が含まれます。音量には十分に注意してください。

低域通過フィルタ（Low-Pass-Filter）

低域通過フィルタを適用すると、カットオフ周波数よりも低い周波数の音を残すようにはたらきます。
今回はカットオフ周波数を600Hzに設定し、一番低い250Hzの音だけが通過するようにします。フィルタ適用後の波形（左上）は波形端が少し歪んでしまっていますが、およそ250Hzの音データだけが残っています。また、同じ波形をフーリエ変換するとちょうど250Hzにピークをもつグラフが得られます。（右上）
フィルタ係数（左下）はフーリエ変換する（右下）と、カットオフ周波数である600Hzよりも低い値はおよそ1、高い値がおよそ0になっています。

※音量には十分に注意してください。

高域通過フィルタ（High-Pass-Filter）

高域通過フィルタを適用すると、カットオフ周波数よりも高い周波数の音を残すようにはたらきます。
先ほどと同様にカットオフ周波数を600Hzに設定し、それよりも高い周波数成分だけが残るようにしました。

※高い周波数の音が含まれます。音量には十分に注意してください。

帯域通過フィルタ（Band-Pass-Filter）

帯域通過フィルタを適用すると。指定した帯域の音を残すようにはたらきます。帯域には600Hz～2500Hzを指定しました。

※高い周波数の音が含まれます。音量には十分に注意してください。

帯域阻止フィルタ（Band-Eliminate-Filter）

帯域阻止フィルタを適用すると、指定した帯域の音を取り除くようにはたらきます。帯域は先ほど同様に600Hz～2500Hzを指定しました。

※高い周波数の音が含まれます。音量には十分に注意してください。

Python+scipyによる実装

FIRフィルタはscipyのsignalモジュールによって実装可能です。
firwin関数またはfirwin2関数でフィルタ係数の配列を作成し、合成波にlfilter関数で適用します。

以下は低域通過フィルタの例です。

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        # 描画する範囲
NUMTAPS = 255       # FIRフィルタの長さ（タップ数）
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)

# FIR filter
lpf_b = signal.firwin(NUMTAPS, CUTOFF, pass_zero=True)
lpf_b_freq = np.fft.fftfreq(len(lpf_b), 1.0 / RATE)
lpf_b_fft = np.fft.fft(lpf_b)
lpf_wave = signal.lfilter(lpf_b, 1, 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, 2, figsize=(8, 4))
axs[0, 0].plot(time[: FRAMES], lpf_wave[: FRAMES])
axs[0, 1].plot(freq[: len(freq)//2], np.abs(lpf_fft)[: len(lpf_fft)//2])
axs[1, 0].plot(lpf_b)
axs[1, 1].plot(lpf_b_freq[: len(lpf_b_freq)//2], np.abs(lpf_b_fft)[: len(lpf_b_fft)//2])

axs[0, 0].set_title("filtered wave")
axs[0, 1].set_title("filtered wave frequency")
axs[1, 0].set_title("filter cofficient")
axs[1, 1].set_title("filter cofficient frequency")
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)

firwin関数のオプション「pass_zero」は0HzからCUTOFFで指定した値の周波数までを通過させるか否か設定することができます。
つまり、「pass_zero＝False」とするだけで上記の例は高域通過フィルタに変えることができます。

また、CUTOFFには２つの要素を持つリストを与えて帯域（バンド）を指定することができます。帯域を指定すれば帯域通過フィルタまたは帯域阻止フィルタを作成することができます。

コメント