人工知能に関する断創録

このブログでは人工知能のさまざまな分野について調査したことをまとめています(更新停止: 2019年12月31日)

ローパスフィルタ

Pythonで音声信号処理(2011/05/14)

FIRフィルタ(2011/10/23)の続きです。今回は、FIRフィルタの代表例であるローパスフィルタ(LPF)を実装していきます。フィルタの設計には、Z変換という技術が必要になります。元の空間では難しい問題をZ変換で別の空間に移す、そこなら簡単に解けて元の空間に戻すといつの間にか複雑な問題が解けている!っておなじみのパターンですね。FIRフィルタの場合、時間領域を周波数領域に移します。時間領域の複雑な畳み込み演算が周波数領域では簡単な掛け算になってしまいます

Z変換

x(n)のZ変換は、

f:id:aidiary:20111028222633p:plain

y(n)のZ変換は、

f:id:aidiary:20111028222634p:plain

で与えられます。元のnの関数をzの関数に変換します。この式がどっからどう出てきてどういう意味があるのか、zとは何なのかと理解に苦しんでいたのですが、今の段階ではあまり考えないことにしよう。この変換を使うとうまくいく、以上(笑)おいおい理解していきたいと思います。

時間領域のFIRフィルタの式は、

f:id:aidiary:20111028222635p:plain

でした。これをY(z)の式にぶち込んでがりがり計算を続け、途中で z \Rightarrow \exp(j 2 \pi f) という置き換えをしてzの関数をさらにfの関数に変換すると

f:id:aidiary:20111028225538p:plain

となります。X(f), Y(f), B(f)は、それぞれx(n), y(n), b(i)のフーリエ変換です。ただし、途中で周波数fをサンプリング周波数fsが1になるように正規化しています。詳しい導出方法は参考文献の書籍を見てください。丁寧に導出されていて非常にわかりやすいです。

この式は、入力信号の周波数特性にフィルタ係数の周波数特性をかけたものが出力信号の周波数特性になっています。時間領域で複雑な畳み込み演算だったFIRフィルタが、周波数領域では簡単な掛け算になってしまいました!フィルタ係数は実際は信号ではないけど、仮に信号とみなしてフーリエ変換してしまえっていうのがミソですね。

試しに前回(2011/10/23)使った移動平均のフィルタ係数 b = [0.5, 0.5] をフーリエ変換して周波数特性を求めてみます。さすがに2つしか値がないとFFTできないので下のように0.0を適当に補います(フーリエ変換の窓幅だけ広げればOK)。

b = [0.5, 0.5]
for i in range(1000):
    b.append(0.0)
fft(b, fs=8000)

フーリエ変換は、高速フーリエ変換(2011/6/18)を参照。結果は、

f:id:aidiary:20111028224921p:plain

となります。上は波形b(i)で下が周波数特性B(f)。b(i)は、本当は波形じゃないけどね。b(i)はすごく見づらいけど最初の2サンプルだけ0.5がたっています。ほかは0.0。問題は下のB(f)で高域成分が減衰している様子がわかります。つまり、入力信号の周波数特性にフィルタ係数の周波数特性をかけあわせると出力信号の周波数が減衰されるって仕組みです。上のB(f)を見るとわかりますが、減衰がすごくなだらかになっています。これを特定のエッジ周波数feのところで急激に減衰するもっと高精度なローパスフィルタを作ろうってのが今回のテーマです(前置き長っ)。

ローパスフィルタの設計

ローパスフィルタってのは、下のようにエッジ周波数feより小さい周波数領域はそのまま通すけど、feより大きい部分の周波数は0にしてしまうようなフィルタです。

f:id:aidiary:20111028231259p:plain

式で書くと

f:id:aidiary:20111028231859p:plain

となります。ナイキスト周波数fs/2は、fs=1と正規化しているので1/2となります。また、FFTの結果は、0を中心として左右対称になるので負の部分も考慮します。先の b = [0.5, 0.5] はなめらかに減衰してましたけど、ローパスフィルタは、エッジ周波数feを境に急激に減衰します。このような特徴を満たすフィルタ係数b(i)を求めるのがローパスフィルタの設計にあたります。

フィルタは周波数を変化させるので周波数領域のB(f)で考えた方が圧倒的にシンプルです。ところが、実際にほしいのは周波数領域のB(f)ではなく、時間領域のフィルタ係数b(i)です。b(i)を得るには、B(f)を逆フーリエ変換してやればよいですね。フーリエ逆変換すると

f:id:aidiary:20111028232250p:plain

となり、sinc(x)というシンク関数が出てきます(シンク関数は2種類の定義がありますが、πがない方です)。b = [0.5, 0.5] に比べれば圧倒的に複雑です。しかも、フィルタ係数の数は無限大・・・これでは畳み込みを計算できないので窓関数を使って有限に打ち切ります。無限大の係数があれば、エッジ周波数feで垂直に減衰します(理想)が、有限で打ち切ると少し滑らかになってしまいます(現実)。フィルタ係数の数(N+1)と滑らか具合=遷移大域幅(δ)の関係は、

f:id:aidiary:20111028232641p:plain

になるそうです。3.1とかどっから出てきたんでしょうね?ここまでくればやっとプログラミングできるので実験してみます。

ローパスフィルタの実験

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

def sinc(x):
    if x == 0.0: return 1.0
    else: return np.sin(x) / x

def fft(x, fs):
    start = 0
    N = 512    # FFTのサンプル数

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

    X = np.fft.fft(x[start:start+N])
    freqList = np.fft.fftfreq(N, d=1.0/fs)

    amplitudeSpectrum = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in X]

    # 波形を描画
    subplot(211)
    plot(range(start, start+N), x[start:start+N])
    axis([start, start+N, -0.5, 0.5])
    xlabel("time [sample]")
    ylabel("amplitude")

    # 振幅スペクトルを描画
    subplot(212)
    n = len(freqList) / 2  # FFTの結果は半分まで見ればOK
    plot(freqList[:n], amplitudeSpectrum[:n], linestyle='-')
    axis([0, fs/2, 0, 2])
    xlabel("frequency [Hz]")
    ylabel("amplitude spectrum")

    show()

def createLPF(fe, delta):
    """ローパスフィルタを設計、fe:エッジ周波数、delta:遷移帯域幅"""
    # 遷移帯域幅を満たすフィルタ係数の数を計算
    # N+1が奇数になるように調整が必要
    N = round(3.1 / delta) - 1
    if (N + 1) % 2 == 0: N += 1
    N = int(N)

    # フィルタ係数を求める
    b = []
    for i in range(-N/2, N/2 + 1):
        b.append(2.0 * fe * sinc(2.0 * math.pi * fe * i))

    # ハニング窓関数をかける(窓関数法)
    hanningWindow = np.hanning(N + 1)
    for i in range(len(b)):
        b[i] *= hanningWindow[i]

    return b

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

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

    # LPFを設計
    fe = 1000.0 / fs        # 正規化したエッジ周波数
    delta = 100.0 / fs      # 正規化した遷移帯域幅
    b = createLPF(fe, delta)

    # フィルタをかける
    y = fir(x, b)
    fft(b, fs)

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

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

窓関数は、ハニング窓(2011/7/16)を使っています(ハミングじゃだめなのかな?)。エッジ周波数を1000Hz、遷移帯域幅を100Hzとしています(いろいろ値を変えると面白いかも)。遷移帯域幅が100Hzの場合、フィルタ係数は248個になります。ここで、NumPyのsinc関数を使うとはまります。NumPyのsinc関数はπがついていますが、ここで使うのはπがついていないシンク関数です。NumPyにはなさそうなので自分で書きました。実行するとフィルタ係数とフィルタ係数の周波数特性が表示されます。

f:id:aidiary:20111028234612p:plain

フィルタ係数b(i)がちゃんとシンク関数の形になっており、B(f)がローパスになっていることも確認できます。入力に与えた音声はホワイトノイズ

です。Audacityというソフトで簡単に作れます。ホワイトノイズの波形と周波数特性は下のようになります。

f:id:aidiary:20111029000132p:plain

一方、ローパスフィルタを通した音声

の波形と周波数特性は、

f:id:aidiary:20111029000212p:plain

となります。エッジ周波数1000Hz以下は通すけど、それ以上は0に減衰します(あれ、低周波の部分はそのまま通さず、値が変わっちゃってるけどあってるのかな?)。波形も高周波のとげとげがなくなって滑らかになっているのがわかります。音声を聴くと、高周波成分がなくなくなるため「しゃー」から「ごー」に変わっています。そーいえば、聴力検査するときによくこの音聴きますね。

簡単にまとめるつもりがかなり長くなってしまいました。次回は、残りのハイパスフィルタ、バンドパスフィルタ、バンドストップフィルタをまとめて実験します。