Python - フーリエ級数展開！
Updated:
Python3 で、フーリエ級数展開を実装する方法についてです。
0. 前提条件
- LMDE 2 (Linux Mint Debian Edition 2; 64bit) での作業を想定。
- Python 3.6.4 での作業を想定。
- 当方は他のバージョンとの共存環境であり、
python3.6,pip3.6で 3.6 系を使用するようにしている。（適宜、置き換えて考えること）
1. アルゴリズムについて
当ブログ過去記事を参照。
2. Python スクリプトの作成
- 敢えてオブジェクト指向で作成している。
- Shebang ストリング（1行目）では、フルパスでコマンド指定している。（当方の慣習）
- 数値計算ライブラリ NumPy を使用しない。（この程度の計算では、逆に2倍程度時間がかかってしまうため）
- 必要であれば、スクリプト内の定数を変更する。（
Nは計算項数）
File: fourier_series_expansion.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
#! /usr/local/bin/python3.6
"""
Fourier's series expansion
f(t) = -1 (-pi < t <= 0 )
1 ( 0 < t <= pi)
"""
import math
import sys
import traceback
class FourierSeriesExpansion:
N = 3 # Number of items
CSV_FILE = "FourierSeriesExpansion.csv"
def expand_fourier_series(self):
""" Fourier's series expansion """
try:
with open(self.CSV_FILE, "w") as f:
f.write("t,f(t)\n")
y = 0
for t in range(int(math.pi * -1000), int(math.pi * 1000) + 1):
for i in range(1, self.N + 1):
y += self.__calc_term(i, t / 1000.0)
ft = 4 / math.pi * y
f.write("{:.3f},{:.6f}\n".format(t / 1000.0, ft))
y = 0
except Exception as e:
raise
def __calc_term(self, n, t):
""" Computation for each item
:param int n
:param int t
:return float
"""
try:
return math.sin((2 * n - 1) * t) / (2 * n - 1)
except Exception as e:
raise
if __name__ == '__main__':
try:
obj = FourierSeriesExpansion()
obj.expand_fourier_series()
except Exception as e:
traceback.print_exc()
sys.exit(1)
3. Python スクリプトの実行
まず、実行権限を付与。
$ chmod +x fourier_series_expansion.py
そして、実行。
$ ./fourier_series_expansion.py
同じディレクトリ内に “FourierSeriesExpansion.csv” という CSV ファイルが作成される。
File: FourierSeriesExpansion.csv
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
t,f(t)
-3.141,-0.002264
-3.140,-0.006083
-3.139,-0.009903
-3.138,-0.013723
-3.137,-0.017542
-3.136,-0.021361
-3.135,-0.025180
-3.134,-0.028999
-3.133,-0.032817
-3.132,-0.036635
-3.131,-0.040452
-3.130,-0.044269
-3.129,-0.048086
-3.128,-0.051901
-3.127,-0.055717
-3.126,-0.059531
-3.125,-0.063345
-3.124,-0.067159
-3.123,-0.070971
:
====< 途中省略 >====
:
3.122,0.074783
3.123,0.070971
3.124,0.067159
3.125,0.063345
3.126,0.059531
3.127,0.055717
3.128,0.051901
3.129,0.048086
3.130,0.044269
3.131,0.040452
3.132,0.036635
3.133,0.032817
3.134,0.028999
3.135,0.025180
3.136,0.021361
3.137,0.017542
3.138,0.013723
3.139,0.009903
3.140,0.006083
3.141,0.002264
4. グラフ化
数字だけを眺めてもよく分からないので、R でグラフ化（プロット）してみるとよいだろう。
グラフは C++, Ruby での実装を紹介した過去記事に掲載のものを参照。
電気工学、音響学、振動、光学等でよく使用する重要な概念です。応用範囲は広いので他にも利用できるかと思います。
以上
Comments