読者です 読者をやめる 読者になる 読者になる

すこしふしぎ.

VR/HI系院生による技術ブログ.まったりいきましょ.(友人ズとブログリレー中.さぼったら焼肉おごらなきゃいけない)

【python】DFT(離散フーリエ変換)してみる【サウンドプログラミング】

こんばんは,1000chです. 前回波形を弄れるようになったので,今回は早速DFTしてみます. 引き続き参考はこちらの書籍.

ではいきます.

離散フーリエ変換の定義

ざっくりメモするとこんな感じ.

音声信号は,1次元の波として表現されています. 前回見た波形は,フレーム毎の振動を表しています. これは"時間軸上の情報"であり,変数tを用いてx(t)と書けます.

このx(t)を,周波数軸で表現する手法がフーリエ変換(Fourier Transform, FT)です. 音声信号を様々な周波数の正弦波の重ね合わせで表現するということです.

{ \displaystyle
X(f) = \int_{-\infty}^{\infty} x(t)\exp(-j2\pi ft)dt \hspace{3.5em} (-\infty \leq f \leq \infty)
}

周波数軸で表現された信号を時間軸に戻すには逆フーリエ変換を用います.

{ \displaystyle
x(t) = \int_{-\infty}^{\infty} X(f) \exp(j2\pi ft)df \hspace{3.5em} (-\infty \leq t \leq \infty)
}

なお,x(t)は連続値の複素数であるとします. しかし,コンピュータ中では音声信号はAD変換され,離散値で扱われています. このままではFTを利用することができません.

そこで離散値を対象として周波数分解できるようにしたのが離散フーリエ変換(Discrete Fourier Transform, DFT)です. 式の導出に関してはここでは触れません. 本を読むなり,ここを参考にしてみて,実際に手を動かしてみるとわかりやすいかと思います. 先に上げた本はなかなかわかりやすかったです(ちょっと変数変換おかしいところがあるような気がしないでも無いけれど).

DFTの定義式はこのようになります.

{ \displaystyle
X(k) = \sum_{0}^{N-1} x(t) \exp( \frac{-j2\pi kn}{N}) \hspace{3.5em} (0 \leq k \leq N-1)
}

その逆変換である逆離散フーリエ変換(Inverse DFT, IDFT)の定義式はこちら.

{ \displaystyle
x(n) = \frac{1}{N} \sum_{0}^{N-1} X(k) \exp( \frac{j2\pi kn}{N}) \hspace{3.5em} (0 \leq n \leq N-1)
}

FTにおいて積分出会った処理がDFTでは総和に書き換えられています. コンピュータで処理しやすくなっていますね.

ちなみに,周波数fは

{ \displaystyle
f = \frac{k f_s}{N}
}

で計算できます.

振幅スペクトルと位相スペクトル

では,入力値がどのように変換されるのか少し考えてみます.

入力値の音声信号は実数であるとします. これをDFTすると,複素数としてX(k)が計算されます.

{
X(k) = A(k) \exp(j\Theta(k))
}

このX(k)をIDFTします.先ほどの式に代入すると,

{ \displaystyle
x(n) = \frac{1}{N} \sum_{0}^{N-1}  A(k) \exp( \frac{j2\pi kn}{N} + j\Theta(k) ) \hspace{3.5em} (0 \leq n \leq N-1)
}

となります.

実際に信号として利用されるのはこのx(n)の実部なので,オイラーの公式により,出力信号x(n)は

{ \displaystyle
x(n) = \frac{1}{N} \sum_{0}^{N-1}  A(k) \cos( \frac{2\pi kn}{N} + \Theta(k)) \hspace{3.5em} (0 \leq n \leq N-1)
}

と表現できるという事になります.

これはすなわち,k番目の周波数の振幅がA(k)ということを表しています.

振幅スペクトルをグラフ化する

以上説明したように,信号の周波数特性を知る為には振幅スペクトルA(f)を求めることが重要です. というわけで,定義に基づいて信号をDFTし,その振幅スペクトルをプロットしてみます.

書籍の参考サイトに置いてあったギター,リコーダーの音を解析してみました.

#coding:utf-8
import wave
import numpy as np
import cmath as cm
import matplotlib.pyplot as plt

def wave_load(filename):
    # open wave file
    wf = wave.open(filename,'r')
    channels = wf.getnchannels() # 追記
    print(wf.getparams())

    # load wave data
    chunk_size = wf.getnframes()
    amp  = (2**8) ** wf.getsampwidth() / 2
    data = wf.readframes(chunk_size)   # バイナリ読み込み
    data = np.frombuffer(data,'int16') # intに変換
    data = data / amp                  # 振幅正規化

    return data


def DFT(data):
    res = []
    N = len(data)
    for k in range(N): #各周波数に関して
        w = cm.exp(-1j * 2 * cm.pi * k / float(N))
        X_k = 0
        for n in range(N): #信号*重みの総和をとる
            X_k += data[n] * (w ** n)
        res.append(abs(X_k))
    return res

if __name__ == '__main__':
    # 波形データ読み込み
    fs = 8000
    wave1 = wave_load('guitar_A4.wav')
    wave2 = wave_load('recorder_A4.wav')

    # DFT.時間かかるので一部のみを利用
    dt1 = DFT(wave1[10000:11024])
    dt2 = DFT(wave2[10000:11024])
    
    # 周波数リストを作成
    # 標本化周波数をデータ数で分割 
    frq = np.arange(1024) * fs / 1024

    # グラフ表示
    # ギターA4
    plt.subplot(2,1,1)
    plt.title('guitar_A4')
    plt.plot(frq,dt1)

    # リコーダーA4
    plt.subplot(2,1,2)
    plt.title('recorder_A4')
    plt.plot(frq,dt2)

    # グラフ表示
    plt.show()

実行結果

f:id:ism1000ch:20140520170126p:plain

こんな感じになりました. DFTではfs/2を境にデータが対象に格納されるという性質があり,解析できるのはfs/2までとなります.

A4は440Hzなので,だいたいその倍音が入っていることが確認できます. また,ギター(上)のほうが倍音成分が多く,リコーダー(下)は基音が強いようです.440Hzとその3倍音である1320Hzあたりにすこし波が存在する程度ですね. 同じA4でも,楽器が違うと含まれる成分が違うことがわかります.

実際の波形を見てみると,

f:id:ism1000ch:20140520171040p:plain

こんな感じになっています.

確かにギターの方がデータのぶれが激しく,高周波成分を多く含むことがわかります.

まとめ

とまぁこんな感じでDFTを実装しました. 楽器によって含まれる周波数成分が異なる,というのを可視化するのはなかなか面白いです.

しかし,動かしてみればわかるのですがDFTはなかなか計算が重く,処理に時間がかかってしまいます. そこで次回は,DFTの計算を工夫し高速化した,FFT高速フーリエ変換)について見ていこうと思います.