FIRフィルタ

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

今回からしばらくディジタルフィルタの実験をいろいろやろうと思います。

ディジタルフィルタは、

• 乗算器
• 加算器
• 遅延器

の3つの要素の組み合わせによって構成されます。この組み合わせ方によって、

• FIR（Finite Impulse Response）フィルタ
• IIR（Infinite Impulse Response）フィルタ

の2種類にわけられるとのこと。今回は、FIRフィルタを実装してみます。

FIRフィルタ

FIRフィルタの定義式は、

となり、畳み込みの定義式と同じです。x(n)は入力信号、y(n)は出力信号、b(i)は乗算器のフィルタ係数、Nは遅延器の数（フィルタ係数の数はN+1）です。式は難しそうだけどシグマをほぐしてみるとけっこう簡単。

b(0)	x(n)	最新の入力信号x(n)にフィルタ係数b(0)をかける
b(1)	x(n-1)	1個前の入力信号x(n-1)にフィルタ係数b(1)をかける
b(2)	x(n-2)	2個前の入力信号x(n-2)にフィルタ係数b(2)をかける
...	...	...
b(N)	x(n-N)	N個前の入力信号x(n-N)にフィルタ係数b(N)をかける

これの和をとったものが最新の出力信号値y(n)になります。つまり、最新の出力信号y(n)を得るためには、最新の入力信号x(n)だけでなく、そのN個前までの信号も重み付けして考慮に入れるってわけですね。そして、i個前の入力信号x(n - i)を得るために遅延器が必要になってくる。なるほど。

実験

ここでは、簡単なFIRフィルタを作ってみます。まず、前に作った（2011/6/7）スクリプトで250Hz、2000Hz、3750Hzの混合正弦波を作ります。一部、抜き出すとこんな感じ。

data = createCombinedWave(0.25, [250, 2000, 3750], 8000, 1.0)
play(data, 8000, 16)
save(data, 8000, 16, "data/combined_sine.wav")

波形のグラフとFFTした結果は下のようになる。

音声は、

波形は2000Hzや3750Hzという高周波成分があるためとげとげしい感じです。聴いた感じでもぴーという高い音が混ざっています。FFTの結果は250, 2000, 3750がちゃんとたっていることがわかります。

ここで、

のフィルタをためしてみます。フィルタ係数bを [0.5, 0.5] としたFIRフィルタです。この式は、n番目の出力信号を得るのに、n番目の入力信号とその1つ前のn-1番目の入力信号の平均をとっているので移動平均をとる式だとわかります。音声波形に対して移動平均をとることで波形がなめらかになり、高周波成分が減衰するとのこと。さっそくスクリプトを書いて試してみます。

FIRフィルタの式に忠実に実装したPythonスクリプト

#coding:utf-8
import wave
import struct
from pylab import *

def fir(x, b):
    """FIRフィルタ
    x: 入力信号
    b: フィルタ係数"""
    y = [0.0] * len(x)  # フィルタの出力信号
    N = len(b) - 1      # フィルタ係数の数
    for n in range(len(x)):
        for i in range(N+1):
            if n - i >= 0:
                y[n] += b[i] * x[n - i]
    return y

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("data/combined_sine.wav", "r")
    fs = wf.getnframes()

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

    # FIRフィルタをかける
    b = [0.5, 0.5]
    y = fir(x, b)

    # 正規化前のバイナリデータに戻す
    y = [int(v * 32767.0) for v in y]
    y = struct.pack("h" * len(y), *y)

    # 音声を保存
    save(y, fs, 16, "fir_combined_sine.wav")

FIRフィルタをかけた波形とFFTの結果は下のようになります。

音声は、

おぉぉぉ、波形のとげとげが丸くなり、甲高いぴーってのがなくなりました。FIRフィルタを適用すると入力信号の周波数特性が変化することがわかります。今回、試したようにb = [0.5, 0.5] のFIRフィルタだと高周波の周波数成分が減衰し、低周波成分だけ通す、いわゆるローパスフィルター（LPF）になっています（もっとよいLPFの作り方はあとで）。

今度は、フィルタ係数の数をもっと増やして移動平均を取ってみます。移動平均の幅を長くするともっとなめらかになるはずだけど・・・フィルタ係数を

b = [0.1] * 10

としてみました。連続する10個のサンプル値で移動平均をとっていることになります。結果は、

確かに高周波成分がさらに減衰し、波形ももっと滑らかになることが確認できました。音声は、

高音のピーがさらに弱まっています。

今回は、フィルタ係数bを決め打ちしましたが、任意の周波数特性に変換するFIRフィルタを作るためには、フィルタ係数bを求める何らかの方法が必要になります。フィルタ係数bを求めることが、フィルタの設計にあたるようです。次回は、FIRフィルタの代表例であるローパスフィルタ（LPF）、ハイパスフィルタ（HPF）、バンドパスフィルタ（BPF）の3つのFIRフィルタを設計してみたいと思います。

参考

• Finite impulse response - Wikipedia
• C言語ではじめる音のプログラミング―サウンドエフェクトの信号処理
• 初心者用 畳み込み（たたみこみ）解説 - こっちは連続なので積分だけど
• やる夫で学ぶディジタル信号処理