7.3 DFT してスペクトルを調べる

Python でのFFTライブラリィは選り取り見取り。

“FFT scipy vs numpy” で検索したら

python - What is the difference between numpy.fft and scipy.fftpack
https://stackoverflow.com/questions/6363154/what-is-the-difference-between-numpy-fft-and-scipy-fftpack

とか出て来た。 scipy と numpy だと、scipy がずっと速いという人がいたり (僕には「ちょっと」としか思えなかった)、 多次元だと numpy の方が速いとか、 FFTW を使う方が速いとか。

スピード競争は好きな人に任せる。 某授業で Mathematica でやったことをほぼそのまま後追いしてみる。

spectrum.py

# coding:utf-8
# spectrum.py --- WAVEファイルを読んで離散フーリエ変換する

# https://docs.python.jp/3/library/wave.html
# http://aidiary.hatenablog.com/entry/20110618/1308367728

import wave
import numpy as np
import scipy.fftpack # サブパッケージは一々 import する
import matplotlib
import matplotlib.pyplot as plt #

# WAVE ファイルをオープンする
wf = wave.open("guitar-5-3.wav", "r")

# サンプリング・レートを得る (Python 的には frame rate という)
sr = wf.getframerate()

# チャンネル数を得る (モノラルか
nchannels = wf.getnchannels()

# サンプルあたりのバイト数
samplewidth = wf.getsampwidth()

# フレームの総数
nframes = wf.getnframes()

# すべてのフレームを読み込む
snd = wf.readframes(nframes) # snd=wf.readframes(wf.getnframes()) とも出来る

# もう閉じていいでしょう
wf.close()

# 
if samplewidth == 2:
    tbl=np.frombuffer(snd, dtype="int16")/32768.0
elif samplewidth == 4:
    tbl=np.frombuffer(snd, dtype="int32")/2147483648.0

# ステレオの場合は左チャンネルだけ
if nchannels == 2:
    tbl = tbl[::2]  # tbl=tbl[1::2] とすると右チャンネル
elif nchannerls == 1:
    print(u"モノラル")
else:
    print(u"あれ？チャンネル数=", nchannels)

# とりあえず先頭から3秒分プロットしてみる
plt.plot(tbl[:3*sr])
plt.show()

# 62800番目から1秒分取り出してプロットする (そうすると周波数簡単)
start=62800
N = sr
plt.plot(tbl[start:start+N])
plt.show()

#
c=scipy.fftpack.fft(tbl[start:start+N])
# c = np.fft.fft(tbl[start:start+N])
plt.plot(abs(c))
plt.show()

# n1からn2までの範囲で図示
def graph(c,n1,n2):
    plt.plot(np.array(range(n1,n2)),abs(c[n1:n2]))
    plt.xlabel('Frequency')
    plt.ylabel('Amplitude Spectrum')
    return

graph(c,1,1600)
plt.show()

graph(c,120,140)
plt.show()

図 5: 3秒分プロット

図 6: 鳴り始めから1秒分プロット

図 7: スペクトル (離散Fourier係数の絶対値)

図 8: 周波数のピークを探す (1600Hzまで)

図 9: