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

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

FFT(高速フーリエ変換)の実験、パラメータなどの解釈(2)

 続きです。グラフの縦軸の目盛りに関しては別の機会に説明するとして、今回は特に横軸の目盛りの読み方について。まずは基本的な事実を確認しておきます。みなさん分かりきっているからなのか、あまりこういう形にまとめられた資料を見ることがありません。自分で、分かりにくい点を含め書いておきます。

●サンプリング周波数をfsと書くことにする。1秒間にこれだけの回数、波形を相手にサンプリングする。前回の記事のコードでもこれを使っている。
フーリエ変換では、もとの波形を(定数と)3角関数の和で表す。(定数分)+(第1次高調波)+(第2次高調波)+(第3次高調波)+……となる。
●具体的な式は以下。

f:id:Inuosann:20191208131311p:plain

あるいは

f:id:Inuosann:20191208131345p:plain
ここでwは基本角周波数。Tをサンプリングされたデータの長さ(時間)としたとき、w=2π/Tである。Tを基本周期、f=1/Tを基本周波数と言う。
●第n次高調波では、基本周期をT、基本周波数をf、基本角周波数をw=2πf=2π/Tとしたとき

f:id:Inuosann:20191208131514p:plain

周期=T/n
周波数(=1/周期)=n/T=nf
角周波数=nw=2πn/T=2πnf
●各高調波で、sin(2π*□*t)の形に書かれた式の□の部分が周波数
●1秒で fs 回サンプリングするから1回のサンプリングに 1/fs 秒かかる。だからN回サンプリングするにはT=N/fs 秒かかる。Tは基本周期。T秒間、音が聞こえている。また基本周波数f=1/T=fs/N。以下、Nは偶数であるとする。
●サンプリング定理により、サンプリング周波数fsはデータに含まれる最大の周波数成分の周波数の2倍以上にしなければならない(そうしないと、サンプリングしたデータからもとのデータを復元できない)。例えば3000ヘルツの周波数成分を調べたいなら、サンプリング周波数は6000ヘルツ以上が必要。
FFTで得られる周波数成分は基本周波数f(=1/T)ヘルツの整数倍。これを「周波数分解能はfヘルツ」と表現する。より正確な(細かな)周波数成分が欲しければfを小さくするしかない。f=fs/Nなので、fsが一定ならNを大きくすればよい(Tを大きくする、と言っても同じ)。
●サンプル数Nのとき、FFTの結果の周波数成分はN個求まる。そのうち、意味のある周波数成分は第N/2-1次高調波まで。例えばN=600(600回サンプリング)なら第299次高調波まで。
FFTでは第n次高調波の周波数の他、離散フーリエ係数Xn(一般には複素数)が得られる。

f:id:Inuosann:20191208131636p:plain

従ってn=N/2 に関して(左右で)離散フーリエ係数の絶対値は対称になる。

f:id:Inuosann:20191228205155p:plain


以上によると今回 N=300なので(コードの yf = sp.fftpack.fft(y[:300]) という行に注目)、基本周波数は1/T=fs/N=600/300=2ヘルツ。従って第k次高調波の周波数は2kヘルツです。使った式は
y = sin(2π*50t)+2*np.sin(2*π*80t)+3*np.sin(2π*100t)
ですから、FFTの結果、50ヘルツ、80ヘルツ、100ヘルツの高調波が得られるはずです。つまり、第25次、40次、50次の高調波が求められるはずです。グラフでは確かに横軸の25、40、50のところにピークがあることが分かります。今回N=600なので、グラフはn=N/2=600/2=300に関して左右で対称です。

 また今回 N=300なので、グラフはn=N/2=300/2=150に関して左右で対称です。

(2020.1.5一部訂正しました)