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

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

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

 少し前、このブログでフーリエ変換について書きました。

www.omoshiro-suugaku.com

 量が多くなりますので何回かに分けますが、実際にFFTのコードを書くにはどうしたらいいのか、出てきた結果をどう解釈したらいいのか、など具体的に書いてみます。備忘録を兼ねていますが、ぼくが引っかかってあれこれ勉強して理解したことをまとめたものです。みなさんの勉強の参考になればうれしいです。まず少しお話を。

 あっちでもこっちでも「(フーリエ変換などの)コードを自分で書くのはダメ。アマチュアはそうしたがる」みたいに言われているようです。ぼくのようなアマチュアには、特に「いつまでにこの仕事を仕上げる」という制限はありません。だから時間をかけて自分で納得のいくまでコードを書けます。ぼくはとにかく理屈を理解したいわけで、それには自分で書くに限ります。ただ、実験室レベルではそれでよくても実用にする段階でまずいことが起こる(桁落ちとか。計算過程で大きな値同士の差などが出てきたときに有効数字が失われる)のは確かでしょう。そこで、ぼくは1回は自分で書いてみて理屈が分かったらそれ以降はライブラリなどを使う、ということでよしとしています。もちろん場合によりますが。

 さて、ライブラリだから使うのも簡単、というわけにはきません。やはり勉強しないとパラメータの解釈などもできません。Pythonも勉強中なので、試してみました。参考文献は、Pythonの利用に関しては『Pythonで学ぶ実践画像処理・音声処理入門』(伊藤克亘ほか2018コロナ社)、フーリエ変換については『これなら分かる応用数学教室―最小二乗法からウェーブレットまで』(金谷健一2003共立出版)です。 

Pythonで学ぶ実践画像・音声処理入門

Pythonで学ぶ実践画像・音声処理入門

 
これなら分かる応用数学教室―最小二乗法からウェーブレットまで

これなら分かる応用数学教室―最小二乗法からウェーブレットまで

 

 Pythonのコードと実行結果は以下の通り。 ほとんど最小限のコードなので分かりやすいと思います。フーリエ変換Pythonのモジュール、scipyのfftpackというのを使っています。便利ですね。fftというのはFFT、「高速フーリエ変換」(fast Fourier transform)の略です。

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

fs = 600
t = np.arange(0, 6, 1/fs)
y = np.sin(2*np.pi*50*t)+2*np.sin(2*np.pi*80*t)+3*np.sin(2*np.pi*100*t)
#y = np.sin(2*np.pi*50*t)
yf = sp.fftpack.fft(y[:300])
plt.figure(1)
plt.plot(np.abs(yf))
plt.show()

f:id:Inuosann:20191208130210p:plain

(周期)関数yを周波数の異なる3角関数の和に分解するのがフーリエ変換です。そして関数そのものではなくサンプリングされたデータをもとにするのが離散フーリエ変換(DFT , discrete Fourier transformです。PC向きです。DFTの高速版がFFTです。例えば y = sin(2πt)+5sin(2π*20t)+7sin(2π*40t) のように分解します。もちろん最初は右辺は分かっておらず、分かっているのはサンプリングされた時点での y の値だけ。サンプリングされたデータから右辺を求めるのです。

 とりあえず今回はここまで。次回はFFTを使う上で理解しておかなければならない基本的な事実を説明します。