続きです。グラフの縦軸の目盛りに関しては別の機会に説明するとして、今回は特に横軸の目盛りの読み方について。まずは基本的な事実を確認しておきます。みなさん分かりきっているからなのか、あまりこういう形にまとめられた資料を見ることがありません。自分で、分かりにくい点を含め書いておきます。
●サンプリング周波数をfsと書くことにする。1秒間にこれだけの回数、波形を相手にサンプリングする。前回の記事のコードでもこれを使っている。
●フーリエ変換では、もとの波形を(定数と)3角関数の和で表す。(定数分)+(第1次高調波)+(第2次高調波)+(第3次高調波)+……となる。
●具体的な式は以下。
あるいは
ここでwは基本角周波数。Tをサンプリングされたデータの長さ(時間)としたとき、w=2π/Tである。Tを基本周期、f=1/Tを基本周波数と言う。
●第n次高調波では、基本周期をT、基本周波数をf、基本角周波数をw=2πf=2π/Tとしたとき
周期=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/2次高調波まで。例えばN=600(600回サンプリング)なら第299次高調波まで。(2022.1.25修正)
●FFTでは第n次高調波の周波数の他、離散フーリエ係数Xn(一般には複素数)が得られる。
従ってn=N/2 に関して(左右で)離散フーリエ係数の絶対値は対称になる。
以上によると今回 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一部訂正しました)