第5回
フーリエ級数
（2021/03/02 最終修正）
内容
フーリエ級数
フーリエ級数
今回は，フーリエ級数の形で表現されたものをグラフとして表示しよう．（最後の課題で数値的にフーリエ級数を求める事もやる。）
細かい事は抜きにして天下り的に，周期2πを持つ周期関数 u(t) は次の級数の形でかくことができる(詳しくはフーリエ解析に関する参考書を見る事)．
これを u(t) のフーリエ級数展開と呼ぶ．ここに，
である．具体的に u(t) が与えられた場合に（２）の積分を計算し，（１）の形の級数を求める事をフーリエ級数展開するいう．フーリエ級数展開するには，（２）の積分計算が必要となる（本ページ最後の課題4）．
ここでは具体例としてフーリエ級数展開された結果（１）を例示し，元の関数 u(t) とともに matplotlibでグラフとして表示するという内容について演習を行う．（１）の級数は無限級数だが実際には無限に足し合わせる事等できない．ある有限個数を足し合わせる事で u(t) を近似する（有限の足し合わせでも結構うまく元の関数を再現する事ができる。このことは音声圧縮MP3等の原理の一部である）．足し合わせ個数が増えるごとに近似が良くなる事を確認する．具体的には，以下の例に取り組む．
ここに示した，２枚目のスライドにある絵を実際に作成する．つまり，
を例えば，N=5 について画面に表示するプログラムを作成する．過去に，指数関数など描いたが，それと基本的に同じ事であることを注意しておく．
課題1（演習）
次の関数を N=5 についてグラフを表示するプログラムを作成せよ．
なお，πは，math モジュールを import した上で、
math.pi
で利用できる．グラフは，-π < t < πについて （上記スライドでは、-2π<t<2π について描いてある） uN(t)を青い実線描き，元の関数を赤い細い線で描くようにせよ．
ヒント： t に関する繰り返し文（過去に sin(t) 等のグラフを描いたのでわかるであろう．t に関して離散化が必要となる．曲線を折れ線で近似的に描くのであった．）と，Σ部分の繰り返し部が入れ子になった２重の繰り返し文が必要となる．これは、Σと繰り返し文の関係を知る良い例になるが、困惑する人も多いので注意して学ぶこと。
念のため書くが，（３）は N=5 の場合，次のようになる．
作成したプログラムを実行すると大体次のようになる．これらの図は，-2π<t<2π を200等分して描いたものである．N=15 程度でほぼ元の関数と重なっている．作成したプログラムについても、-2π<t<2π で描く事もチャレンジして欲しい。これができないとなると、ここまでの内容がわかっていないということだろう。（uNは簡単に -2π<t<2π で描けるが、 u(t) はそのままではだめなので工夫が必要である。）
Code5-1: フーリエ級数
# Code5-1.py
import math
import matplotlib.pyplot as plt
T = 2*math.pi
M = 100
dt = T/M
N = 100
def fu(t):
return abs(t)
t, u, uN = [0.0]*(M + 1), [0.0]*(M + 1), [0.0]*(M + 1)
for i in range(0, M + 1):
t[i] = i*dt - math.pi
u[i] = fu(t[i])
for i in range(0, M + 1):
uN[i] = math.pi/2
for i in range(0, M + 1):
for n in range(1, N + 1, 2):
uN[i] = uN[i] - 4/math.pi*(math.cos(n*t[i])/(n*n))
plt.plot(t, u, color="blue", linewidth=1, label="$u(t)$")
plt.plot(t, uN, color="red", linewidth=2, label="$u_N(t)$")
plt.title("5-1 (N=" + str(N) + ")")
plt.xlabel("t")
plt.ylabel("u(t)")
plt.legend()
plt.show()
課題2
課題1で作成したプログラムで N を大きくした場合に，元の関数をより良く近似することを視覚的に確認せよ．
課題3
他の周期関数についてもフーリエ級数展開を行い，それをグラフとして表示するプログラムを作成せよ．
周期2Lを持つ周期関数 u(t) は次の級数の形でかくことができる（なぜこうなるかは各人参考書等で確認してください)．
ここに，
であります．具体的に u(t) が与えられた場合に（２L）の積分を計算し，（１L）の形の級数を求める事をフーリエ級数展開するという．フーリエ級数展開するには，（２L）の積分計算が必要だが，それは皆さん得意であろう。ただし、u(t)が例えばデータとして与えられた場合には、通常その積分は難しい。数値的な積分であれば可能であり、フーリエ級数を求める事ができる。
（２L）の積分計算を数値的に行う方法は、離散フーリエ変換、FFT等が調べる場合のキーワードで調べると良い。例えば音声データなどをフーリエ級数展開すると、どういう波数の波が多く含まれているか等の情報が得られる上に、MP3等の音楽の圧縮アルゴリズムもフーリエ級数展開の応用である。
課題4（演習）
課題1の問題について（２）を数値的に計算するデモを行う。本来はFFT等のアルゴリズムを使うのだろうが、ここでは素朴に数値積分を行うデモだ。直ぐにFFTだ！と例題コードをコピーしにいくのではなく、まずは素朴にやってみるのはプログラミングの練習としては良いだろうと思う。
実用上はこのような単純な方法では無く、より計算量が少なくなるよう工夫された方法FFT等が使われる。調べると良い。なるほど、うまいやり方があるもんだなと思うだろう。
Code5-4 : 数値積分でフーリエ級数を求めてみる
# Code5-4.py
import math
import matplotlib.pyplot as plt
T = 2*math.pi
M = 1000
dt = T/M
N = 50
#def fu(t):
# return abs(t)
def fu(t):
if t < 0:
return -1.0
else:
return 1.0
a, b = [0.0]*(N + 1), [0.0]*(N + 1)
t, u, uN = [0.0]*(M + 1), [0.0]*(M + 1), [0.0]*(M + 1)
for i in range(0, M + 1):
t[i] = i*dt - math.pi
u[i] = fu(t[i])
for i in range(0, M):
a[0] = a[0] + u[i]*dt
a[0] = a[0]/(2*math.pi)
for n in range(1, N + 1):
for i in range(0, M):
a[n] = a[n] + u[i]*math.cos(n*t[i])*dt
b[n] = b[n] + u[i]*math.sin(n*t[i])*dt
a[n] = a[n]/math.pi
b[n] = b[n]/math.pi
for i in range(0, M + 1):
uN[i] = a[0]
for i in range(0, M + 1):
for n in range(1, N + 1):
uN[i] = uN[i] + a[n]*math.cos(n*t[i]) + b[n]*math.sin(n*t[i])
for n in range(0, N + 1):
print(n, a[n], b[n])
plt.plot(t, u, color="blue", linewidth=1, label="$u(t)$")
plt.plot(t, uN, color="red", linewidth=2, label="$u_N(t)$")
plt.title("5-4 (N=" + str(N) + ")")
plt.xlabel("t")
plt.ylabel("u(t)")
plt.legend()
plt.show()
上記コード中に、
u[i] = fu(t[i])
というところがあるが、ここで例えば、u[i] に対して音声の波形データを入れれば、音声データのフーリエ解析ができる事になる。求まった a[i], b[i] からパワースペクトルを得たり、一部を 0 にした上で、元に戻すこと（Codeの後半）を行えばノイズを減らすとか音質を変えることができるのだ。
課題5
課題4で作成したコードを一般の区間幅に対応できるように変更せよ。