![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
2013/12/23からのアクセス回数 4364
ここで紹介したSageワークシートは、以下のURLからダウンロードできます。
http://www15191ue.sakura.ne.jp:8000/home/pub/31/
また、Sageのサーバを公開しているサイト（http://www.sagenb.org/, http://www15191ue.sakura.ne.jp:8000/）にユーザIDを作成することで、ダウンロードしたワークシートを アップロードし、実行したり、変更していろいろ動きを試すことができます。
Interface 2013/11から連載が始まった「実験で入門！音声合成のメカニズム」をSageを使ってお復習いしていきます。
この説明は、私のブログの lbed/LM4F120 Launchpadで音声合成 と連携しています。 こちらも参考にしてください。
フーリエ級数は、周期的に繰り返される信号をSinとCosで表現したものです。
基本周期\(T_0\)として、周期信号を式で表すと次のようになります。 $$ f(t) = a_0 + \sum_{k=1}^\infty \left( a_k cos \frac{2 \pi k}{T_0} t + b_k sin \frac{2 \pi k}{T_0} t \right) $$
Interface 2013/11号では、以下の様な三角波を例題としています。この波形はY軸対称（偶関数）なので、求めるフーリエ係数の\(b_k\)が0となります。
Sageで三角波を表現する場合には、区分関数（Piecewse）を使用します。 三角波では、区間を-0.5から0と0から0.5の2つに分け、それぞれf1, f2で関数を定義します。
sageへの入力:
# [-0.5, 0.5]の区間で以下の関数を定義
n, t = var("n t")
f1 = lambda t: 4*t + 1
f2 = lambda t: -4*t + 1
f = Piecewise([[(-1/2, 0), f1], [(0, 1/2), f2]])
plot(f, figsize=4)
基本周期\(T_0\)とするとフーリエ係数は、以下の様に求めることができます。 $$ \left\{ \begin{eqnarray} a_0 & = & \frac{1}{T_0} \int_{-T_0/2}^{T_0/2} f(t) dt \\ a_k & = & \frac{1}{T_0} \int_{-T_0/2}^{T_0/2} f(t) cos \frac{2 \pi k}{T_0} t dt, k = 1, 2, ... \\ b_k & = & \frac{1}{T_0} \int_{-T_0/2}^{T_0/2} f(t) sin \frac{2 \pi k}{T_0} t dt, k = 1, 2, ... \end{eqnarray} \right. $$
さっそく、数式処理システムのSageを使って三角波のフーリエ級数を求めてみましょう。 残念ながらPiecewiseの関数に直接sin, cosを掛けることができないので、 sin, cosを掛けた区間関数を定義して、それを積分することにします。
以下がSageでの処理です。最後に$a_0, a_1, a_2, a_3, b_1$を出力しています。
この結果は、Interface 2013/11号の式(4)と同じ結果になっているのが確認できました。 $$ \left\{ \begin{eqnarray} a_k & = & \left\{ \begin{eqnarray} 0, & k: 偶数 \\ \frac{8}{(k \pi)^2}, & k: 奇数 \\ \end{eqnarray} \right. \\ b_k & = & 0, k: すべての整数 \end{eqnarray} \right. $$
sageへの入力:
# Piecewiseの関数にsin, cosを掛けることができないので(たぶんバグ)、
# フーリエ級数を求める積分は、sin, cosを掛けた区分関数を使って積分する
T0 = 1
k = var('k')
f1s = lambda k: f1(t)*sin((2*pi*k)/T0*t)
f2s = lambda k: f2(t)*sin((2*pi*k)/T0*t)
f1c = lambda k: f1(t)*cos((2*pi*k)/T0*t)
f2c = lambda k: f2(t)*cos((2*pi*k)/T0*t)
b(k) = (2/T0)*Piecewise([[(-1/2, 0), f1s(k)], [(0, 1/2), f2s(k)]]).integral(t, -T0/2, T0/2)
a(k) = (2/T0)*Piecewise([[(-1/2, 0), f1c(k)], [(0, 1/2), f2c(k)]]).integral(t, -T0/2, T0/2)
a0 = (2/T0)*Piecewise([[(-1/2, 0), f1c(0)], [(0, 1/2), f2c(0)]]).integral(t, -T0/2, T0/2)
# a0 = Piecewise([[(-1/2, 0), f1], [(0, 1/2), f2]]).integral(t, -T0/2, T0/2)
# とするとTypeErrorとなります（これもバグ？）
# 値の確認
print a0, a(1), a(2), a(3), b(1)
0 8/pi^2 0 8/9/pi^2 0
求めた\(a_k\)を使ってkの値を1, 5, 19と変えた結果をプロットしてみます。
sageへの入力:
def fk(t, k):
return sum(a(n)*cos((2*pi*n/T0)*t) for n in range(1, k+1))
show(fk(t, 3))
# n=1,5,19の合成波をプロット
p = plot(fk(t, 2), -1,1, color='blue',legend_label='n=1')
p += plot(fk(t, 5), -1,1, color='red',legend_label='n=5')
p += plot(fk(t, 10), -1,1, color='green',legend_label='n=19')
p.show(figsize=5)
Piecewiseには、フーリエ級数を計算する機能が備わっています。 先のfkに相当する関数がfourier_series_partial_sumです。
fourier_series_partial_sumの使い方は、以下の様に使います。
g.fourier_series_partial_sum(N, L)
ここで、Nは計算する次数、Lは区間の長さL（区間は、-LからLです）で、fourier_series_partial_sumの返す値は、以下の式で表されます。 $$ f(x) \sim \frac{a_0}{2} + \sum_{n=1}^N [a_n\cos(\frac{n\pi x}{L}) + b_n\sin(\frac{n\pi x}{L})] $$
それでは、同じ問題をfourier_series_partial_sumで計算してみましょう。
sageへの入力:
# 同じ問題をPiecewiseのフーリエ級数関数を使って解いてみる g = Piecewise([[(-1/2, 0), f1], [(0, 1/2), f2]]) plot(g, figsize=4)
sageへの入力:
show(g.fourier_series_partial_sum(5, 1/2))
計算結果をプロットしてみます。N=100とするとほとんど三角波になっています。
sageへの入力:
p = plot(g.fourier_series_partial_sum(5, 1/2), -1,1, color='blue',legend_label='n=5') p += plot(g.fourier_series_partial_sum(10, 1/2), -1,1, color='red',legend_label='n=10') p += plot(g.fourier_series_partial_sum(50, 1/2), -1,1, color='green',legend_label='n=50') p += plot(g.fourier_series_partial_sum(100, 1/2), -1,1, color='purple',legend_label='n=100') p.show(figsize=5)
皆様のご意見、ご希望をお待ちしております。
