今回も基本シリーズ。
以前ウワンが記述した記事は、FFTに主眼がなく、実際に応用するとき参考にしずらいので、もっと簡単に具体的な例で記述しようと思う。
とにかく、正しい周波数を返すように注意したいと思います。
【参考】
①【Scipy】FFT、STFTとwavelet変換で遊んでみた♬～②不確定原理について～
この記事もFFTの基本は理解しているものとして、コード解説を中心に話を進めます。
コード解説
from scipy import signal
import matplotlib.pyplot as plt
import numpy as np
import wave
from scipy.fftpack import fft, ifft
import matplotlib
今回は、昨夜の記事で作成したwavfileを利用する。
wavfile = 'output400_50'
以下がwavfile読込の関数です。
def wave_input(wavfile):
wr = wave.open(wavfile+'.wav', "rb")
ch = wr.getnchannels() #monoral ch=1
width = wr.getsampwidth() #width=2
fr = wr.getframerate() #sampling freq ; RATE
fn = wr.getnframes() #sampling No. of frames
fs = fn / fr #sampling time
data = wr.readframes(fn)
sig = np.frombuffer(data, dtype="int16") /32768.0
t = np.linspace(0,fs, fn) #, endpoint=False)
return t,sig,fn,fr
t,sig,fn,fr=wave_input(wavfile)
描画するために座標等を定義します。
matplotlib.rcParams.update({'font.size': 18, 'font.family': 'sans', 'text.usetex': False})
fig = plt.figure(figsize=(12,12)) #(width,height)
axes1 = fig.add_axes([0.1, 0.55, 0.8, 0.4]) # main axes
axes2 = fig.add_axes([0.15, 0.7, 0.2, 0.2]) # insert axes
axes3 = fig.add_axes([0.1, 0.1, 0.8, 0.35]) # FFT axes
信号を描画します。描画方法は昨夜と同様な描画です。
axes1.plot(t, sig)
axes1.grid(True)
axes1.set_xlim([0, 0.1])
#plt.pause(1)
axes2.set_ylim(-1.2,1.2)
axes2.set_xlim(0,5)
axes2.plot(t,sig)
#plt.pause(1)
#plt.savefig("./fft/sig_fig_"+wavfile+".jpg", dpi=200)
以下がFFT関数です。
def FFT(sig,fn,fr):
freq =fft(sig,fn)
Pyy = np.sqrt(freq*freq.conj())/fn
f = np.arange(0,fr,fr/fn)
以下でピーク振動数と振幅を求めます。
order=10は（両サイドの）連続する10個のデータ中相対値が一番大きな値を求めます；10回で勝ち抜け
order : How many points on each side to use for the comparison to consider comparator(n, n+x) to be True.
ld = signal.argrelmax(Pyy, order=1) #相対強度の最大な番号をorder=10で求める
ssk=0
fsk=[]
Psk=[]
maxPyy=max(np.abs(Pyy))
for i in range(len(ld[0])): #ピークの中で以下の条件に合うピークの周波数fと強度Pyyを求める
if np.abs(Pyy[ld[0][i]])>0.25*maxPyy and f[ld[0][i]]<20000: # and f[ld[0][i]]>20:
fssk=f[ld[0][i]]
Pssk=np.abs(Pyy[ld[0][i]])
fsk.append(fssk)
Psk.append(Pssk)
ssk += 1
print('{}'.format(np.round(fsk[:len(fsk)],decimals=3))) #標準出力にピーク周波数fskを小数点以下二桁まで出力する
print('{}'.format(np.round(Psk[:len(fsk)],decimals=4))) #標準出力にピーク強度Pskを小数点以下6桁まで出力する
return freq,Pyy,fsk,Psk,f
freq,Pyy,fsk,Psk,f=FFT(sig,fn,fr)
いよいよFFTの結果描画です。
※抽出したピークの近傍を描画させています
Pyy_abs=np.abs(Pyy)
axes3.plot(f,Pyy_abs)
axes3.axis([min(fsk)*0.5, max(fsk)*2, 0,max(Pyy_abs)*1.5])
axes3.grid(True)
axes3.set_xscale('log')
axes3.set_ylim(0,max(Pyy_abs)*1.5)
以下は必須ではありませんが、求めたピークをグラフに記述しているところです。
タイトルと図中のピークにそれぞれ記載するようにしました。
axes3.set_title('{}'.format(np.round(fsk[:len(fsk)],decimals=1))+'\n'+'{}'.format(np.round(Psk[:len(fsk)],decimals=4)),size=10) #グラフのタイトルにピーク周波数とピーク強度を出力する
axes3.plot(fsk[:len(fsk)],Psk[:len(fsk)],'ro') #ピーク周波数、ピーク強度の位置に〇をつける
# グラフにピークの周波数をテキストで表示
for i in range(len(fsk)):
axes3.annotate('{0:.1f}'.format(fsk[i]), #np.round(fsk[i],decimals=2) でも可 '{0:.0f}(Hz)'.format(fsk[i])
xy=(fsk[i], Psk[i]),
xytext=(10, 20),
textcoords='offset points',
arrowprops=dict(arrowstyle="->",connectionstyle="arc3,rad=.2")
)
plt.pause(1)
plt.savefig('./fft/figure_'+wavfile+'.jpg')
【参考】
②Pythonで高速フーリエ変換（FFT）の練習-5 周波数ピークを自動で検出
結果
まず、基準振動の400Hz、少数第一位まで正しく表示できました。
そして、昨夜と同じようにGifアニメーションで示すと以下のようになりました。
※近接した周波数の重畳を分離したくて、途中描画の最大値と最小値の決め方を変更したので、少し変な動きになっています。
結果の評価
１．二つの周波数を入力としたとき、１Hzの差は分離できている
２．20Hz以下の低周波も1Hzまではきちんと取得できているが、0.5Hz以下は取得できていない
３．周波数は少数第一位まで正しく表示できている
コードは以下に置いた
まとめ
・比較的簡単なコードで正しくFFTできるようになった
・一連の話として可聴音以下の低周波でも重畳する音の周波数と強度が適切ならば、人が感知できることになる
⇒これが実際のマイク（やスピーカー）で実現できるのかを検証する
・フォルマント合成を本格的に実施する
・stftとwaveletで分析する