ぽきたに 〜ありきたりな非凡〜

現役F欄大学生が送るゴミ溜めと独り言

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

どうもたっきーです。

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

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

 

 

参考

aidiary.hatenablog.com

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軸に強度をグラフにプロットしたものを周波数スペクトルと呼ぶ。

 

例↓↓↓(一応、強度については0〜1で正規化してある。)

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のサンプル数(2**n)
    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

 んで、各引数は以下の通り。

  1. グラフに入れたいWAVファイルの数
  2. 楽器のディレクトリ名(例. "piano")
  3. FFTのサンプル数 (2の乗数)
  4. 乱数の開始位置
  5. 乱数の終点位置

それでは、スペクトルを一個だけ表示してみる。

IN

fft_load(1,"piano",1024,0,1000)

 OUT

f:id:tacky0612:20171212173325p:plain

複数表示してみる。

IN

fft_load(7,"Eguitar",1024,0,1000)

OUT

f:id:tacky0612:20171212173500p:plain

うまく表示できました✋

 

サウンドプログラミング入門――音響合成の基本とC言語による実装 (Software Design plus)

サウンドプログラミング入門――音響合成の基本とC言語による実装 (Software Design plus)

 
入門 Python 3

入門 Python 3

 
退屈なことはPythonにやらせよう ―ノンプログラマーにもできる自動化処理プログラミング

退屈なことはPythonにやらせよう ―ノンプログラマーにもできる自動化処理プログラミング

 

 

 

 

 

おわりに・感想

この後、色んな楽器の周波数スペクトルを表示してみたが、目視では楽器ごとの周波数スペクトルの規則性みたいなのはあまり感じられなかった。(あくまで主観。)

しかし、低音楽器であるベースは周波数の低い位置に周波数スペクトルが固まっていた。↓

f:id:tacky0612:20171212173923p:plain

って、まぁ当たり前なんだが。

 

 あ、この前に鹿と戯れてきましたカワイイ。

 

 

スポンサードリンク