人工知能に関する断創録

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

離散フーリエ変換

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

今回は、信号処理の肝とも言える離散フーリエ変換(Discrete Fourier Transform: DFT)を試してみようと思います。ときどき感動するアルゴリズムに出会うけれど、フーリエ変換はその一つです。最初に考え出したフーリエさんはすごい!フーリエ変換を扱った本は参考文献に挙げている何冊かを読んだのですが、理解するのにけっこう苦労しました。ここでも間違ったこと書いていたらコメントもらえると助かります。

前回の正弦波の合成(2011/06/07)で試したように、任意の周期波形はさまざまな周波数を持つ正弦波の合成で表せます。フーリエ変換は各周波数の正弦波がどれくらいの割合で含まれているかを求める技術。ここら辺の定性的な理解は、

の説明が大変わかりやすかったです。まず、フーリエ変換の式と逆フーリエ変換の式、


X(f) = \displaystyle \int_{-\infty}^\infty x(t) \exp(-j 2 \pi f t) dt \hspace{100pt} (-\infty \le f \le \infty)
x(t) = \displaystyle \int_{-\infty}^\infty X(f) \exp(j 2 \pi f t) df \hspace{100pt} (-\infty \le t \le \infty)

暗記したいのできれいに書いてみました(笑)時間tも周波数fも両方とも連続値なので積分を使ってます。数学の本だとこっちの式を使って解析的にがしがし解くような例題がよく出てますね。これだとコンピュータでどう実装するかいまいちわからなかった覚えがあるんですが、コンピュータで扱うときは時間tも周波数fも離散化して離散フーリエ変換という形式に変換します。


t = n t_s \hspace{100pt} (0 \le n \le N-1)
f = k f_s / N \hspace{100pt} (0 \le k \le N-1)

ここで、tsは標本化周期、fsは標本化周波数(8000Hz、16000Hz、22050Hz、44100Hzなど)。ts = 1/fsの関係が成り立つ。Nはサンプル数。離散化すると新しい変数nとkが出てくる。nは離散化した時間のインデックス(n番目の時間)、kは離散化した周波数のインデックス(k番目の周波数)です。この記号が本によって違ってややこしい。時間の離散化は比較的わかりやすかったけど周波数も離散化されるので注意。周波数の場合、fs/Nが1単位となってその整数倍(k次高調波)の値しか取れない。サンプル数Nが大きいほど1単位であるfs/Nが小さくなるので周波数分解能が高くなる。ここら辺は後で実験で確かめてみます。で、この離散化を導入すると離散フーリエ変換(DFT)と逆離散フーリエ変換(IDFT)の式が得られる。


X(k) = \displaystyle \sum_{n=0}^{N-1} x(n) \exp \Bigl( \frac{-j 2 \pi k n}{N} \Bigr) \hspace{100pt} (0 \le k \le N-1)
x(n) = \frac{1}{N} \displaystyle \sum_{k=0}^{N-1} X(k) \exp \Bigl( \frac{j 2 \pi k n}{N} \Bigr) \hspace{100pt} (0 \le n \le N-1)

こっちも暗記したいのできれいに書いてみました(笑)一般的にX(k)は複素数になる。一方、x(n)は波形の値なので実数。というわけでさっそくPythonでDFTを実装してみます(numpymatplotlibのインストールが必要)。上のDFTの式をそのまま実装しただけです。Pythonは複素数がそのまま使えるのですごい便利ですね。

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

def dft (start, x, N):
    """xのstartサンプル目からのNサンプルを周期波形とみなしたDFT値を返す"""
    X = [0.0] * N
    for k in range(N):
        for n in range(N):
            real = np.cos(2 * np.pi * k * n / N)
            imag = - np.sin(2 * np.pi * k * n / N)
            X[k] += x[start + n] * complex(real, imag)
    return X

if __name__ == "__main__" :
    wf = wave.open("data/sine.wav" , "r" )
    fs = wf.getframerate()  # サンプリング周波数
    x = wf.readframes(wf.getnframes())
    x = frombuffer(x, dtype= "int16") / 32768.0  # -1 - +1に正規化した波形
    print len(x)
    wf.close()

    start = 0        # サンプリングする開始位置
    N = 256          # サンプル数
    X = dft(start, x, N)    # 離散フーリエ変換
    freqList = [k * fs / N for k in range(N)]    # 周波数のリスト
    amplitudeSpectrum = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in X]    # 振幅スペクトル np.abs()
    phaseSpectrum = [np.arctan2(int(c.imag), int(c.real)) for c in X]      # 位相スペクトル np.angle()

    # 波形サンプルを描画
    subplot(311)  # 3行1列のグラフの1番目の位置にプロット
    plot(range(start, start+N), x[start:start+N])
    axis([start, start+N, -1.0, 1.0])
    xlabel("time [sample]")
    ylabel("amplitude")

    # 振幅スペクトルを描画
    subplot(312)
    plot(freqList, amplitudeSpectrum, marker='o', linestyle='-')
    axis([0, fs/2, 0, 15])    # ナイキスト周波数まで表示すれば十分
    xlabel("frequency [Hz]")
    ylabel("amplitude spectrum")

    # 位相スペクトルを描画
    subplot(313)
    plot(freqList, phaseSpectrum, marker='o', linestyle='-')
    axis([0, fs/2, -np.pi, np.pi])
    xlabel("frequency [Hz]")
    ylabel("phase spectrum")

    show()

startサンプル目から始まるNサンプルを対象にDFTを求めています。DFTはこの抽出したNサンプルが周期波形であることを仮定しています。なので、取り出した区間が周期波形かを確認するため波形も表示してみました。サンプルの取り出し方によっては端っこが滑らかにつながらないこともあるためハニング窓を使うことが多いようです。これは後で試してみます。使った正弦波のデータ(sine.wav)は、振幅0.25、基本周波数250Hz、サンプリング周波数8000Hz、量子化ビット数16です。この設定の正弦波WAVEファイルの作り方は前回の正弦波の合成(2011/06/07)を参照してください。実行すると、

f:id:aidiary:20110611073340p:plain

となります。250Hzの成分だけがもとの正弦波に含まれていることがわかりました。大成功!DFT値は、k=N/2にあたるナイキスト周波数fs/2を中心として対称構造を取るため横軸の周波数はfs/2まで取れば十分です。今回は、サンプリング周波数8000Hzなのでナイキスト周波数は4000Hzまでです。ちなみにナイキスト周波数を超えてfsまで全部描画すると

f:id:aidiary:20110611073339p:plain

のようになります。振幅スペクトルは線対称、位相スペクトルは点対称になってます。説明が飛びましたが、dft()が返すXには複素数が入っています。複素数はそのままでは上のようにプロットできないため振幅スペクトル位相スペクトルというのを求めます。普通、スペクトルといったら振幅スペクトルの方ですが、位相スペクトルも信号を復元する上で必要な情報です(参考文献によると人間は位相スペクトルの違いに鈍感だそうであまり重視されてこなかったとのこと)。振幅スペクトルは複素数のDFT値の絶対値(np.absで計算可)、位相スペクトルは複素数のDFT値の偏角(np.angleで計算可)になります。振幅スペクトルを2乗したパワースペクトルを使うこともあります。


A[k] = \sqrt{real(X[k])^2 + imag(X[k])^2}
\Theta[k] = \tan^{-1} \frac{imag(X[k])}{real(X[k])}

Pythonで逆正接を求めるにはnumpyのarctan2(y, x)を使います(mathのatan2()でもOK)。虚数軸の方が第1引数にくるので注意!もう1つ非常に悩んだのだけど上のプログラムで式の通りに

phaseSpectrum = [np.arctan2(c.imag, c.real) for c in X]

とやってしまうと無茶苦茶な位相スペクトルが出てきてはまる。その理由は、浮動小数点で正確に0を表せないためでした。逆正接は角度を返すけれど、微妙にxとyが変わるだけで角度は大きく変わってしまうわけです。下のような感じ。

>>> np.arctan2(0, 0)
0.0
>>> np.arctan2(1, 0)
1.5707963267948966   # pi/2
>>> np.arctan2(-1, 0)
-1.5707963267948966  # -pi/2
>>> np.arctan2(0.000001, 0.000002)
0.46364760900080609
>>> np.arctan2(0.000002, 0.000002)
0.78539816339744828

なので仕方なくint型に変換してからarctan2()に渡しました。int型にしてしまっていいのか微妙だけど・・・もしいい方法あったら教えてください。

いろんな波をフーリエ変換

正弦波だけでは面白くないので他の音声信号もDFTでスペクトルを求めてみましょう。前に作ったノイズ入り正弦波(f0=250Hz)、三角波(f0=250Hz)、矩形波(f0=250Hz)、ノコギリ波(f0=250Hz)とギター(f0=440Hz)のスペクトルを見てみます。

f:id:aidiary:20110611073341p:plain

ノイズ入り正弦波。基本周波数の250Hzの成分はいっぱい、高周波ノイズの成分が微妙に混じるって感じですね。ローパスフィルタを使うと高周波ノイズが除去できます。これは後で試してみたいです。

f:id:aidiary:20110611073343p:plain

三角波。基本周波数の250Hzの成分が一番多く、基本周波数の3倍の750Hz、5倍の1250Hz、7倍の1750Hzの成分が徐々に減衰しながら含まれていることがわかります。これは三角波の正弦波合成式を見ると明らかです。


x(t) = \frac{8}{\pi^2} \sum_{k=0}^\infty (-1)^k \frac{\sin((2k+1) \omega t)}{(2k+1)^2}

基本周波数 w =2πfの2k+1倍の波が重ね合わさってます。この2k+1は、1、3、5、7といった奇数になります。また、分母に(2k+1)^2があるようにkが大きい高周波成分ほど減衰されます。詳しくは、正弦波の合成(2011/06/07)を参照してください。

f:id:aidiary:20110611073342p:plain

矩形波。基本周波数の250Hzの成分が一番多く、基本周波数の3倍の750Hz、5倍の1250Hz、7倍の1750Hzの成分が徐々に減衰しながら含まれていることがわかります。三角波と同じですが、減衰が三角波より緩やかです。これも矩形波の正弦波合成式を見ると明らかです。


x(t) = \frac{4}{\pi} \sum_{k=1}^\infty \frac{\sin((2k-1) \pi f t)}{(2k-1)}

三角波は分母がべき乗だったので減衰が激しいけど矩形波は分母がべき乗じゃないですね。

f:id:aidiary:20110611073338p:plain

ノコギリ波。基本周波数の250Hzの成分が一番多く、基本周波数の2倍の500Hz、3倍の750Hz、4倍の1000Hz、5倍の1250hzの成分が徐々に減衰しながら含まれていることがわかります。矩形波、三角波と違って偶数倍の周波数も含まれることがわかります。これもノコギリ波の正弦波合成式を見ると明らかです。


x(t) = \frac{2}{\pi} \sum_{k=1}^\infty \frac{\sin(2 \pi kft)}{k}

f:id:aidiary:20110611073344p:plain

ギターの音。普通の楽器の音だけど人工的に作った三角波、矩形波、ノコギリ波のようにちゃんと基本周波数440Hzの整数倍の周波数を含む倍音構造があるのがわかります。この倍音がどれくらい含むかによって音色が決まるとのこと。

最後に、サンプル数Nを変えてみましょう。今まではN=256でやってましたが、N=64に落としてみます。すると、

f:id:aidiary:20110611084514p:plain

となって周波数の離散化が粗くなることがわかります。周波数分解能が落ちてしまったわけです。

  • サンプル数Nが大きい ⇔ 周波数分解能が高くなる ⇔ 長い区間のスペクトルなので時間分解能は低い
  • サンプル数Nが小さい ⇔ 周波数分解能が低くなる ⇔ 短い区間のスペクトルなので時間分解能は高い

が成り立つとのこと。周波数分解能と時間分解能はトレードオフ。一般的に、スペクトルの時間的変化を知りたい場合は、サンプル数Nを小さくしたほうがよいそうです。実際はどれくらいのNがよく使われるんでしょう?

まとめ

今回はDFTの基本式に忠実に計算しましたがこの方法は遅いです。よく音楽プレイヤーで波打っているようなスペクトルアナライザは、DFTを高速化した高速フーリエ変換(Fast Fourier Transform: FFT)というアルゴリズムが使われるそうです。これは、numpyやscipyにも実装されているので次回試してみたいと思います。

f:id:aidiary:20110611084515p:plain

参考文献

フーリエの冒険

フーリエの冒険

  • 作者: トランスナショナル・カレッジ・オブ・レッ
  • 出版社/メーカー: ヒッポファミリークラブ
  • 発売日: 1988/07
  • メディア: 単行本
  • 購入: 7人 クリック: 126回
  • この商品を含むブログ (23件) を見る
やり直しのための信号数学―DFT、FFT、DCTの基礎と信号処理応用 (ディジタル信号処理シリーズ)

やり直しのための信号数学―DFT、FFT、DCTの基礎と信号処理応用 (ディジタル信号処理シリーズ)

今日から使えるフーリエ変換 (今日から使えるシリーズ)

今日から使えるフーリエ変換 (今日から使えるシリーズ)

C言語ではじめる音のプログラミング―サウンドエフェクトの信号処理

C言語ではじめる音のプログラミング―サウンドエフェクトの信号処理