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

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 でやったことをほぼそのまま後追いしてみる。

sound1.py

# sound1.py
# coding:utf-8
# 原始的なやり方で WAVE ファイルを扱う

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

# sudo port install py27-scipy
# sudo port install py27-matplotlib これに +gtk2 つける?
#   ~/.matplotlib/matplotlibrc に何か書く。例えば
#   backend : TkAgg
# sudo port install py27-pygtk

# import モジュール名 とすると、モジュール名のみがシンボルテーブルに入る
#   モジュール名.関数名 でアクセス
#   好きな名前 = モジュール名.関数名 という小技もある
#
# from モジュール名 import * はモジュールで定義されている名前をすべて入れる
#

import wave
import numpy as np
import scipy.fftpack # サブパッケージは一々 import する
import matplotlib
#matplotlib.use('GTKAgg')# GTK GTKAgg tkagg とかあるけど、GTKAgg しか動かん
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()

桂田 祐史
2018-01-08