SciPyのFIRフィルタの使い方 - 人工知能に関する断創録

Pythonで音声信号処理（2011/05/14）

今まで、5回に渡ってFIRフィルタの例であるローパスフィルタ、ハイパスフィルタ、バンドパスフィルタ、バンドストップフィルタを実装してきました。内部で何をやっているか勉強になって大変面白かったんだけど、当然ながらSciPyにはすでに実装されているので今回はそちらを使ってみます。使うのはたった2つの関数だけ。使い方も簡単です。

以下のスクリプトは、ローパス、ハイパス、バンドパス、バンドストップを全部まとめています。メイン関数にあるお好みのフィルタのコメントをはずせばそれぞれ試せます。

#coding:utf-8
import struct
import wave
import numpy as np
import scipy.signal
from pylab import *

"""SciPyのFIRフィルタ関数を使うサンプル"""

def fft(b, y, fs):
    """フィルタ係数bとフィルタされた信号yのFFTを求める"""
    b = list(b)
    y = list(y)

    N = 512    # FFTのサンプル数

    # 最低でもN点ないとFFTできないので0.0を追加
    for i in range(N):
        b.append(0.0)
        y.append(0.0)

    # フィルタ係数のFFT
    B = np.fft.fft(b[0:N])
    freqList = np.fft.fftfreq(N, d=1.0/fs)
    spectrum = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in B]

    # フィルタ係数の波形領域
    subplot(221)
    plot(range(0, N), b[0:N])
    axis([0, N, -0.5, 0.5])
    xlabel("time [sample]")
    ylabel("amplitude")

    # フィルタ係数の周波数領域
    subplot(223)
    n = len(freqList) / 2
    plot(freqList[:n], spectrum[:n], linestyle='-')
    axis([0, fs/2, 0, 1.2])
    xlabel("frequency [Hz]")
    ylabel("spectrum")

    # フィルタされた波形のFFT
    Y = np.fft.fft(y[0:N])
    freqList = np.fft.fftfreq(N, d=1.0/fs)
    spectrum = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in Y]

    # 波形を描画
    subplot(222)
    plot(range(0, N), y[0:N])
    axis([0, N, -1.0, 1.0])
    xlabel("time [sample]")
    ylabel("amplitude")

    # 振幅スペクトルを描画
    subplot(224)
    n = len(freqList) / 2
    plot(freqList[:n], spectrum[:n], linestyle='-')
    axis([0, fs/2, 0, 10])
    xlabel("frequency [Hz]")
    ylabel("spectrum")

    show()

def save(data, fs, bit, filename):
    """波形データをWAVEファイルへ出力"""
    wf = wave.open(filename, "w")
    wf.setnchannels(1)
    wf.setsampwidth(bit / 8)
    wf.setframerate(fs)
    wf.writeframes(data)
    wf.close()

if __name__ == '__main__':
    wf = wave.open("whitenoise.wav", "r")
    fs = wf.getframerate()

    x = wf.readframes(wf.getnframes())
    x = frombuffer(x, dtype="int16") / 32768.0

    nyq = fs / 2.0  # ナイキスト周波数

    # フィルタの設計
    # ナイキスト周波数が1になるように正規化
    fe1 = 1000.0 / nyq      # カットオフ周波数1
    fe2 = 3000.0 / nyq      # カットオフ周波数2
    numtaps = 255           # フィルタ係数（タップ）の数（要奇数）

    b = scipy.signal.firwin(numtaps, fe1)                         # Low-pass
#    b = scipy.signal.firwin(numtaps, fe2, pass_zero=False)        # High-pass
#    b = scipy.signal.firwin(numtaps, [fe1, fe2], pass_zero=False) # Band-pass
#    b = scipy.signal.firwin(numtaps, [fe1, fe2])                  # Band-stop

    # FIRフィルタをかける
    y = scipy.signal.lfilter(b, 1, x)

    # フィルタ係数とフィルタされた信号のFFTを見る
    fft(b, y, fs)

    # 音声バイナリに戻して保存
    y = [int(v * 32767.0) for v in y]
    y = struct.pack("h" * len(y), *y)
    save(y, fs, 16, "whitenoise_filtered.wav")

SciPyの信号処理関係の関数は、scipy.signalの下にまとまっています。FIRフィルタを設計する関数は、firwin()という1つにまとめられていて引数を変えることでローパス、ハイパス、バンドパス、バンドストップが作れました。firwin()の返り値は、フィルタ係数です。

SciPyの場合、今まで実装してきたのと微妙に言い回しや作法が違います。今までエッジ周波数といってきたのは、SciPyではカットオフ周波数、フィルタ係数の数（N+1）はタップ数（numtaps）と書かれています（いろいろ検索してみたらこっちのほうが一般的みたいだ）。あと、今までの実装では式の導出の関係でサンプリング周波数fsが1.0になるように正規化してきたけれど、SciPyはナイキスト周波数fs/2が1.0になるような正規化が必要。

pass_zeroという引数は、0Hz（一番左）を通すか通さないかを意味します。これによってローパスとハイパス、バンドパスとバンドストップが区別できます。たとえば、ハイパスフィルタはカットオフ周波数より左は通さない（0に減衰）ので一番左にある0Hzは通さない、つまり、pass_zero=Falseです。省略するとTrueになります。ここら辺は最初はうまい実装だなと思ったのですが、けっこうややこしいので素直にlowpass()、highpass()とか作ればいいのにと思えてきました。

遷移帯域幅の指定はwidthという引数がそれっぽいけどよくわからない。あと打ち切りに使う窓関数はwindowという引数がある。デフォルトはハミング窓。違うのはそれくらいかな？

実際に信号にフィルタをかける関数は、lfilter()です。引数にフィルタ係数と入力信号を与えるとフィルタをかけた信号が返ってきます。第2引数は重みのようでどう使うのかいまいち理解できなかったけど上の4つのフィルタをつくるときは1を指定します。

実行結果です。まず、ローパスフィルタ。左上がフィルタ係数の波形領域b(i)、左下がその周波数領域B(f)、右上がフィルタをかけたホワイトノイズの波形領域y(n)、右下がその周波数領域Y(f)です。左下を見ればどんなフィルタかわかりますね。

f:id:aidiary:20111102224258p:plain