# coding:utf-8
# spectrum.py --- WAVEファイルを読んで離散フーリエ変換する
#
# https://m-katsurada.sakura.ne.jp/lecture/fourier-2024/spectrum.py
# https://m-katsurada.sakura.ne.jp/lecture/fourier-2024/guitar-5-3.wav

# 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()