いぬおさんのおもしろ数学実験室

おいしい紅茶でも飲みながら数学、物理、工学、プログラミング、そして読書を楽しみましょう

PythonでscipyのFFT、IFFTを使った分かりやすい例(2)

 前回、FFTの例をお見せしました。さらに理解を確実なものにするため、もう少し別のパラメータで確認してみましょう。タイトルは「FFTとIFFT」となっていますが、今回はFFTのみで説明します。コードは以下の通りとします。

import numpy as np
import scipy as sp
from scipy.fftpack import fft
import matplotlib.pyplot as plt
from matplotlib import pyplot

pi = np.pi

#基本周期T = 3秒。基本周波数=1/3Hz。
fs = 40 #サンプリング周波数。データ1個で1/40秒。
N = 120 #データ数は3/(1/40)=120
#tは1/40秒刻みの時刻の配列。1/40秒ごとに3秒間、サンプリング

t = np.arange(0, 3, 1/40)   

y = np.sin(2*pi*4*t)+2*np.sin(2*pi*10*t)    #周波数は4Hz,10Hz
pyplot.plot(t, y)     #元の波形を描画
pyplot.show()

#yf[k]には第k次高調波のフーリエ係数が入る
yf = sp.fftpack.fft(y[:N])

pyplot.plot(t, np.abs(yf)) #含まれる周波数成分を描画
pyplot.show()

  データ数は120です。するとFFTでは第120次高調波まで得られます(が、そのうち意味があるのは第59次高調波までなのでした)。この音をT=3秒、流しています。このTが基本周期です。すると基本周波数はf=1/T=1/3Hzです。第1次高調波、第2次高調波、第3次、第4次、……、第120次高調波の周波数は順に1/3Hz、2/3Hz、3/3Hz、4/3Hz、……、120/3Hz=40Hz(=fs)です。だから、得られたFFTの結果の横軸にはこれらが目盛られていると考えます。横軸の左端は0Hz、右端は40Hz(=fs)です。

元の波形はこれです。

f:id:Inuosann:20191220211125p:plain

これに対してFFTを実行すると次のようになります。

f:id:Inuosann:20191220211321p:plain

右端の「3.0」のところには120/3=40Hzの成分(=0)が現れています(第120次高調波の周波数です)。例えば「2.5」のところには単純に比例配分して40×(2.5/3.0)=33.3Hzの成分が現れます(これも0ですが)。

左から2つ目のピークは、拡大すると以下の通りです。

f:id:Inuosann:20191220211837p:plain

「0.75」のところのピークの周波数はやはり単純に比例配分して40×(0.75/3)=10Hzです。これがコード通りの10Hzのピークなのです。

まとめましょう。

データ数をN、それだけの音が流れるに要する時間をT秒(従って基本周波数は1/Tヘルツ)とするときのFFTの結果の右端は第N次高調波の周波数成分でfsヘルツ(=N/Tヘルツ)。そこより左側は、単純に比例配分して周波数を求められる。

 なお、横軸に周波数を表示させることはすぐできます。右端を40Hz(=fs)にすればよいのですから、コードの下から2行目を

pyplot.plot(t * 40 / 3, np.abs(yf)) #含まれる周波数成分を描画

と書き換えればOKです。以下のように表示されます。

f:id:Inuosann:20191220214253p:plain

 一応、これでFFTの基本的な部分については迷わないはずです。おわかりかとは思いますが、うまくいった例はスニペットとして残しておき、再利用しましょう。