Pythonでフーリエ級数展開、その処理の背景

2019/03/31 2019/11/04

大学の卒業研究で，フーリエ級数展開のプログラムを実装することがありました。

フーリエ変換であれば，numpy.fft.fft(f)で簡単に高速フーリエ変換ができるのですが，フーリエ級数展開の場合はモジュールが存在していないっぽかったので，自分で作る必要がありました。

特別なモジュールや記法はないと思うので、少しいじればPythonに限らず使えます。

数式

フーリエ級数展開は，

$\displaystyle f(t)=\sum_{n=1}^\infty (a_n \cos n \omega t+b_n \sin n \omega t)$

で表される式で，nは展開の次数，ωは角速度(=2πf=2π/T)です。

フーリエ係数an, bn は

$\displaystyle a_n=\frac{2}{T} \int_{-T/2}^{T/2} f(t) \cos n \omega t dt$

$\displaystyle b_n=\frac{2}{T} \int_{-T/2}^{T/2} f(t) \sin n \omega t dt$

で計算されます。

コード

	import numpy as np

	def fourier(i_data, t_data):
	kaisu = 160 # フーリエ級数展開の回数

	def cos(t, n):
	return np.cos(n * 2 * np.pi / T * t) # cos(n * ωt)
	def sin(t, n):
	return np.sin(n * 2 * np.pi / T * t)

	i_gu = i_ki = 0
	for n in range(1, kaisu + 1):
	i_cos = i_data * cos(t_data, n)
	an = (2 / T) * i_cos.sum() * dt
	i_gu += an * cos(t, n)

	i_sin = i_data * sin(t_data, n)
	bn = (2 / T) * i_sin.sum() * dt
	i_ki += bn * sin(t, n)
	return [i_gu, i_ki]

	# 本体----------------------------------------

	Current = フーリエ級数展開したいデータを配列で用意
	T = フーリエ級数展開したいデータの周期を用意

	point = len(Current)
	dt = 任意
	t = np.arange(0, dt * point, dt)

	# フーリエ級数展開の実行
	i = fourier(Current, Time)[0] + fourier(Current, Time)[1]

view raw FourierExpansion.py hosted with ❤ by GitHub

説明

cos(nωt), sin(nωt) は係数の計算と展開で2回使うので，関数内関数として定義しました。

あとは大体数式通りの処理をしているだけですが，「フーリエ係数を求めるときのt」と「f(t)を求めるときのt（プログラムではTime）」は使っているものが異なります。
これで測定したデータ間の値を予測できるようになるためです。詳しくは後述。

積分について

積分についてはsum()で単純に合計をとっていますが，これは元データがオシロスコープ（電流を測定できる装置）で測った離散的なデータだったからです。

下図の左のように，元のデータが離散的であれば単純に合計をとったほうが正確な値が出ますが，連続的な数式やモデルをプログラムに落とし込むにあたって仕方なく離散的なデータになってしまっている場合は下図の右のように，scipyのようなモジュールでcumtrapz（台形則）なり，simps（シンプソン則）なり，quad（ガウス求積）なりの積分をしたほうが正確な値が出ます。

しかし，連続的な数式やモデルが既に分かっている場合にフーリエ級数展開をしたくなるような状況って，あまりないと思うので，基本的には単純に合計とるのがよいです。

フーリエ級数展開をして何がうれしいのか

プログラムでデータを扱うにあたって，フーリエ級数展開をする理由は以下だと思っています。

偶関数成分と奇関数成分の分解が行える
測定したデータ間の値を（それっぽく）予測できるようになる
高周波のノイズ除去ができる

偶関数成分と奇関数成分の分解が行える

以下のように，よく分からない変化をしている波形を偶関数と奇関数に分解することで，変化の状況を把握しやすくなったりします。

私の専攻だった電気系では，有効電力とか何とかの概念があって，一周期でみると電力を消費するのは偶関数だけということになっていました。

だから偶奇の分解をすることは重要だったのかなと思います（小並）

測定したデータ間の値を（それっぽく）予測できるようになる

離散的なデータって，測定した微小時間以上の近くまで拡大すると，線形補完みたいな感じでカクカクしたグラフになってしまいます。

しかし，フーリエ級数展開をすることによってこういったデータも連続的なものにできるので，拡大しても滑らかにそれっぽく変化したグラフにすることができます。

あとは，単純に物理的な制約（※）で測定データ間のポイント数が合わなくなることってあるので，その対策にもなります。
※ 私の場合，周波数をパラメータとしてオシロスコープの分解能の限界をまたぐような範囲で電流のデータを測定していた

そのため，本プログラムではdtを任意の値にすれば時間軸も自由に変化させられるようにしています。

高周波のノイズ除去ができる

これはある種のデメリットかもしれませんが，実用的なフーリエ級数展開には次数の限界があるので，元のデータと完全に同じ数式を作るということは不可能になります。

なので，高周波のノイズ（かもしれない）データが欠落したグラフを作れるということです。

フーリエ級数展開の中身

$a_n \cos n \omega t+b_n \sin n \omega t$

は三角関数の和なので， $\tan \theta_n=a_n/b_n$ として

$\sqrt{a_n^2+b_n^2} \sin (\theta_n+n \omega t)$

のように表すことができるので，下図のようなスペクトラムを作れます。

この図だと，n=50くらい以上の展開はほとんどノイズのようだから展開の次数を減らして綺麗なグラフを描こう，みたいな判断ができます。

こういうノイズの除去はフーリエ変換と逆変換を使って振幅と周波数の両面から行う場合が多いと思うのですが，フーリエ級数展開でも同じことができそうですね。

私が深く理解してないので，数学的には同じことをしてるのかもしれないですが，ちょっとよく分からないです（投げやり）

-IT, 電気電子
-Python, ツール