【Python】WAVファイルの波形データにFFTかけて周波数スペクトルを複数表示する【サウンドプログラミング】

どうもたっきーです。

最近の悩みは虫歯です；；

これは前回からの続きなので前回のヤツやってないと動作しません（）

参考
なにする？
環境
FFTとは
周波数スペクトルとは
窓関数とは
コード
おわりに・感想

参考

aidiary.hatenablog.com

窓関数を用いる理由 - ロジカルアーツ研究所

tacky0612.hatenablog.com

なにする？

WAVファイルのデータにFFTかけて周波数解析し、matplotlibで周波数スペクトルを複数表示する。

調べてみたところ、FFTかけるときに窓関数というものをかけるのが一般的らしいのでそれも実装。

環境

macOS High Sierra 10.13
Python 3.6.3
numpy 1.13.3
matplotlib 2.1.0

FFTとは

高速フーリエ変換（fast Fourier transform）とは、離散フーリエ変換（discrete Fourier transform, DFT）を計算機上で高速に計算するアルゴリズムのこと。

ちなみにDFTを直接計算するとめちゃめちゃ時間かかる。。。

そこで考えられたのがFFTってわけ。

めちょ早い。

数学が苦手なので中身の数式には触れない。

Numpyを使えばnp.fft.fft()で一瞬で求められる。

周波数スペクトルとは

光や音や電磁波信号は様々な周波数の成分から構成されている。そのようなものから周波数毎の強さを定量的に求める処理をスペクトル解析（spectrum analysis）と呼ぶ。

んで、そのスペクトル解析で求められたものを、x軸に周波数、y軸に強度をグラフにプロットしたものを周波数スペクトルと呼ぶ。

例↓↓↓（一応、強度については０〜１で正規化してある。）

f:id:tacky0612:20171212170734p:plain

窓関数とは

こんな関数↓

f:id:tacky0612:20171212171401p:plain

参考

僕の窓関数に対する理解は、『数値解析が容易になるべんりツール』って感じ。(簡単になるってニュアンスとはまたちょっと違う)

まぁ色んな窓関数があるわけだが、王道はハミング窓らしいので今回は無心でハミング窓を使用。

窓関数を用いる理由についてはここが分かりやすかった。

コード

まずは関数の定義。

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

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

    # 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                  # 振幅正規化
    data = data[::channels]
   
    return data


def fft_load(count,instrument,size,start,end):
    '''
    count = グラフに入れたい数
    instrument = 楽器のディレクトリ
    size = FFTのサンプル数（２＊＊ｎ）
    start　=乱数の開始位置
    end = 乱数の終点位置

    '''
    st = 10000   # サンプリングする開始位置
    hammingWindow = np.hamming(size)    # ハミング窓
    fs = 44100 #サンプリングレート
    d = 1.0 / fs #サンプリングレートの逆数
    freqList = np.fft.fftfreq(size, d)
    
    for i in range(count):
        n = random.randint(start,end)
        filename = "./" + instrument + "/output/" + str(n) +".wav"
        wave = wave_load(filename)
        windowedData = hammingWindow * wave[st:st+size]  # 切り出した波形データ（窓関数あり）
        data = np.fft.fft(windowedData)
        data = data / max(abs(data)) # 0~1正規化
        plt.plot(freqList,abs(data))
         
    plt.axis([0,fs/16,0,1]) #第二引数でグラフのy軸方向の範囲指定
    plt.title(instrument)
    plt.xlabel("Frequency[Hz]")
    plt.ylabel("amplitude spectrum")
    plt.show()

    return data