自由研究の準備（その１１）Pythonで音をフーリエ変換(FFT)して分析します

はじめに

さて、自由研究の素材を「画像:MNIST」から「音:Sound」に変更すると宣言してから幾年月。
ようやく大詰めを迎えました。いよいよ大御所にしてラスボス「フーリエ変換」FFTの登場です。前回のlibrosaでも内部的には大活躍でしたが、直接対決です。

おおまかな内容

音を解析し、再構築するには「フーリエ変換」（もしくは高速フーリエ変換：FFT)を使うことが多いです。

フーリエ変換は $F(j\omega)=\int_{-\infty}^{\infty}f(t)e^{-j\omega t}dt$

離散フーリエ変換は $F(j\omega)=\sum_{n=-\infty}^{\infty}f(nT_{0})e^{jnT_{0}\omega}$

...
みなかったことにしましょう。
この数式の意味は分からずとも、ライブラリに準備された高速フーリエ変換を実際に使い、音を周波数の数値に変換してみます。また、元の音に戻してみます。

いきなりやってみました

元にする波形は準備（その５）でお世話になったプログラムで作成した９倍音まで含むノコギリ波です。

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

if __name__ == "__main__":
    # あらかじめ準備した ノコギリ波 C音（基底 261.626Hz:9倍音まで)を読み込み
    wf = wave.open("Music/saw-sample-c.wav", "r")
    fs = wf.getframerate()  # サンプリング周波数 上のサンプルでは48K
    x = wf.readframes(wf.getnframes())
    x = frombuffer(x, dtype="int16") / 32768.0  # -1 - +1に正規化
    wf.close()

    start = 0  # サンプリングする開始位置
    N = 2 ** 12  # FFTのサンプル数  2**12は4096です 

    # サンプリング数にあったハニング窓を準備します。
    win = signal.hann(N)

    # FFTを実行してスペクトラム情報を獲得します。
    spectrum = scipy.fftpack.fft(x[start:start + N] * win)
    freqList = scipy.fftpack.fftfreq(N, d=1.0 / fs)

    # 振幅と位相をグラフ表示するために抽出します。
    amplitudeSpectrum = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in spectrum]  # 振幅スペクトル
    phaseSpectrum = [np.arctan2(int(c.imag), int(c.real)) for c in spectrum]  # 位相スペクトル

    # スペクトラム情報から元の波形に復元します。
    resyn_sig = ifft(spectrum)

    # 1.元の波形を描画
    figure(figsize=(8, 8))
    subplot(511)  # 5行1列のグラフの1番目の位置にプロット
    subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.2, hspace=0.5)
    plot(range(start, start + N), x[start:start + N])
    axis([start, start + N, -1.0, 1.0])
    xlabel("Fig 1. Original Wave / time [sample]")
    ylabel("amp")

    # 2.窓関数を通した波形を描画
    subplot(512)  # 5行1列のグラフの2番目の位置にプロット
    plot(range(start, start + N), x[start:start + N] * win)
    axis([start, start + N, -1.0, 1.0])
    xlabel("Fig 2. through window Wave / time [sample]")
    ylabel("amp")

    # 3.振幅スペクトルを描画
    subplot(513)
    plot(freqList, amplitudeSpectrum, linestyle='-')
    axis([0, 3000, 0, 600])
    xlabel("Fig 3.  frequency [Hz]")
    ylabel("amp spect.")

    # 4.位相スペクトルを描画
    subplot(514)
    plot(freqList, phaseSpectrum, linestyle='-')
    axis([0, 3000, -np.pi, np.pi])
    xlabel("Fig 4.  frequency [Hz]")
    ylabel("phase spect.")

    # 5.元の波形を描画
    subplot(515)  # 5行1列のグラフの5番目の位置にプロット
    plot(range(start, start + N), resyn_sig[start:start + N])
    axis([start, start + N, -1.0, 1.0])
    xlabel("Fig 5. resynth Wave / time [sample]")
    ylabel("amp")

    show()

ちなにみ、ちょっとだけscipyのFFTのほうが、numpyのFFTよりも早いという噂がありscipyのFFTを使っています。

図１.元にする波形ののこぎり波です

あらかじめ別のpythonプログラムで生成したノコギリ波（に近い）の波形です。

少し低いド(C) で基底の周波数は261.626Hzです。９倍音までを重ねています。

f:id:np2LKoo:20160924172808p:plain

音も貼ってみます。4096サンプルなのでとても短いです。

図２.窓関数を通した波形

そのまま使うと両端でクリッピングノイズが発生するので、窓関数で両端をゼロにしてやります。(正確にはサンプリングの「窓」（上記プログラムでは N ）の範囲で周期関数であると仮定して計算するため、終端がゼロでないと元の信号にない成分が導出されるのでこれを防ぐためためです）

f:id:np2LKoo:20160924172857p:plain

これも音を貼ってみます。立ち上がり、立下りがクリップノイズなどがなく、スムーズに聞こえると思います。一層短く聞こえます。

図３.スペクトラムの振幅情報

フーリエ変換(FFT)して得たスペクトラムの振幅情報です。グラフで９倍音までちゃんと検出されていることがわかります。

f:id:np2LKoo:20160924172934p:plain

図４.スペクトラムの位相情報

フーリエ変換(FFT)して得たスペクトラムの位相情報です。おもいのほか重要な情報です。

f:id:np2LKoo:20160924172956p:plain

図５.スペクトラムから元に戻した波形

スペクトラム情報からFFT逆変換(scipy.ifft())を使い元の波形にしています。

f:id:np2LKoo:20160924173046p:plain

ほぼ完ぺきなくらいに復元できてます。ここまで復元できるのは技術的には当たり前なのでしょうが、感覚的にはまるで魔法のようです。元の音と違いはわからないですが、いちおう貼っておきます。

ちなみに、窓関数をつかわないとこうなります。

f:id:np2LKoo:20160924173118p:plain

フーリエ変換のまとめ

フーリエ変換は元になる波形情報からある範囲(窓)を計算し、以下を導出します

周波数
振幅
位相

窓関数を入れて計算することが多いです。

窓関数には用途に応じいろいろな種類がありますが、正弦波、または正弦波の組み合わせの解析にはハニング窓関数がよいとのことでした。

以下がハニング窓関数です。

$w(n) = 0.5−0.5 \cos ( {\frac{2 \pi n}{M−1}} )$ $\ \ \ 0\ \leq \ n \ \leq \ M−1$

...
みなかったことにしましょう。

ちなみにグラフはこちらです。

f:id:np2LKoo:20160924193410p:plain

ハニング関数はscipy.signal.hanning()を使えます。scipy.signal.hann()も同じです。似た名前のハミング関数scipy.signal.hamming()は似ているけど別のものなので注意しましょう。（ハミングは両端がゼロでない関数です）

ハニング（ハン窓）関数は気象学者のJulius von Hannにより名づけられました。

ハミング関数はR. W. Hammingにより名づけられました。

低域の正確な検出には多くのサンプルが必要です

いろいろな周波数で確認していると問題点がでてきました。
高域ではまったく問題ない解像度なのですが、低域から中域にかけてサンプル数が少ないため分解能がひくいのです。
たとえば、サンプリング周波数が48kでサンプル数が1024の場合はfreqList=[0.0, 46.875, 93.75, ..]と約46.8Hz刻みです。これでは低域の半音検出もできません。

低音の「ド」は、C=32.703, C#=34.648 で、差は1.945Hz です。sin波の音をCかC#を特定するにはこの周波数付近では1Hzの解像度はほしいところです。
とするとサンプル数（＝窓の大きさ）Nは上限（＝サンプリング周波数）の48000までいってしまいました。これで低音の「ド」を解析させると、1Hz刻みで
amplitudeSpectrum[30] = 90.05 (30.0Hz)
amplitudeSpectrum[31] = 474.23 (31.0Hz)
amplitudeSpectrum[32] = 4315.14 (32.0Hz)
amplitudeSpectrum[33] = 5665.88 (33.0Hz)
amplitudeSpectrum[34] = 1734.12 (34.0Hz)
なので、原音C=32.703Hzは12半音識別できます。(ちなみにB=30.868, C#=34.648)

中域はA=440.000Hz, G#=415.305Hz, A#=466.164 で14Hz程度の識別は最低必要です。例えば「窓」サイズN = 4096 (2**12) で11.72Hz刻みになります。

高域A=3520.000Hz, G#=3322.438Hz, A#=3729.310Hz で209Hz程度の識別は最低必要です。例えば「窓」サイズN = 256 (2**8) で187.5Hz刻みになります。

低音の検出の「窓」は思いっきり広くとり、中音では程々に、高音では結構狭くてもかまわない感じです。（ただ、ビブラートとかチョーキングとか、フレットレスベース感を出したいので精度をどれくらいにするか、各種資源や処理速度とのせめぎあいです）どのように実装に反映するかは、実験しながら決めたいです。