TensorFlow始めました その4。方形波をフーリエ級数展開して、その係数を学習する。

TensorFlow始めました。
TensorFlow始めました その2。placeholderを使う。
TensorFlow始めました その3。1次直線の傾きと切片を機械学習で求める。
の続きです。

今回は矩形波をフーリエ級数展開して、その係数を機械学習で求めてみます。

1. フーリエ級数とは

フーリエ級数(フーリエきゅうすう、Fourier series)とは、複雑な周期関数や周期信号を、単純な形の周期性をもつ関数の(無限の)和によって表す方法である。フーリエ級数は、フランスの数学者ジョゼフ・フーリエによって金属板の中での熱伝導に関する研究の中で導入された。

Wikipediaより引用

関数\(f(x)\)が周期関数であれば、次式のように、定数(右辺第1項)と三角関数の無限和(右辺第2項)で表現できてしまう…という法則です。

\[f(x) = \frac{a_0}{2} + \sum_{k=1}^\infty(a_k\cos kx + b_k\sin kx) \tag{1}\]

2. 方形波のフーリエ級数展開

今回のお題である方形波(周期\(2\pi\)、振幅\(\pm\pi/4\))を次のように定義します。

\[
f(x) = \begin{cases}
-\pi / 4 & (-\pi \leq x < 0) \\
\pi / 4 & (0 \leq x < \pi)
\end{cases}
\tag{2}
\]

これは奇関数ですのでフーリエ級数展開したとき、\(a_0=0\)、\(a_k=0\) (\(k = 1,2, … \infty\))になります。
また、\(k\)が偶数のとき\(b_k=0\)になり、結局、\(k\)が奇数のときのsinの項のみが残ります。そこで、(1)の余計な項を消去して、次のように書き直します。
\[
f(x) = \sum_{k=1}^{\infty}b_{2k-1}\sin{ \left (\left(2k-1 \right)x \right)}
\tag{3}
\]

無限和を計算できないので、\(k\)の上限を\(k_{max}\)とし、方形波を近似したのが次式。

\[
f(x) \approx \sum_{k=1}^{k_{max}}b_{2k-1}\sin{ \left (\left(2k-1 \right)x \right)}
\tag{4}
\]

\(b_{2k-1}\) \((k=1, 2, 3, …, k_{max})\)がフーリエ係数。それでは、方形波のフーリエ係数を機械学習で求めてみましょう。

3. フーリエ係数を機械学習するPythonコード

Pythonによるコーディング。
TensorFlow他、必要なパッケージをインポート。

 

TensorFlowで(4)式を記述します。
学習対象のフーリエ係数\(b\)の初期値は\(\pm1\)の範囲の乱数(一様分布)。

 

続けて教師データを用意します。
\(-\pi \leq x<0 \)で\(y=-\pi/4\)、\(0 \leq x<\pi \)で\(y=\pi/4\)です。

 

学習部分です。

 

グラフ描画。
学習前後のグラフを描画します。

 

最後に全コードを掲載。

 

4. 機械学習!

実行しました。繰り返し学習でlossが徐々に減少して12.53ぐらいに収束。

 

理屈では,(2)式の方形波のフーリエ係数は次表のとおり。
\(k=n\)のとき係数は\(b_{2n-1}=1/(2n-1)\)になります。

\(k\)\(2k-1\)フーリエ級数 \(b_{2k-1}\)
\(1\)\(1\)\(1/1\)
\(2\)\(3\)\(1/3\)
\(3\)\(5\)\(1/5\)
\(n\)\(2n-1\)\(1/(2n-1)\)
\(k_{max}\)\(2k_{max} -1\)\(1/(2k_{max}-1)\)

 

学習結果は,\(b\)がほぼ理屈どおりに、1/1, 1/3, 1/5, 1/7, ….になりました。

 

学習前後のグラフは次のとおり。
学習前は、ぐちゃぐちゃな波形。学習後はかなり方形波に近づいています。

 

\(k_{max}=100\)で実行してみました。結果は次のとおり。
学習後はほぼ方形波。ギブスの現象のツノもしっかり。

 

5. 余談:方形波のフーリエ級数展開から円周率\(\pi\)を求める。

方形波のフーリエ級数展開は(3)式。そのフーリエ係数は\(b_{2k-1}=1/(2k-1)\)。
なので、こうなります。
\[
f(x) = \sum_{k=1}^{\infty}b_{2k-1}\sin{\left (\left(2k-1 \right)x \right)} \\
= \frac{1}{1}\sin(1x) + \frac{1}{3}\sin(3x) + \frac{1}{5}\sin(5x) + …
\tag{5}
\]

\(x=\pi/2\)だったらどうなるでしょうか?
(2)式から\(f(\pi/2)=\pi/4\)。さらに(5)式より次式が得られます。
これはライプニッツの級数。美しい!
\[
f(\pi/2) = \frac{1}{1}\sin(\frac{1\pi}{2}) + \frac{1}{3}\sin(\frac{3\pi}{2}) +\frac{1}{5}\sin(\frac{5\pi}{2}) + …\\
= \frac{1}{1} – \frac{1}{3} + \frac{1}{5} – \frac{1}{7} + …
= \frac{\pi}{4}
\tag{6}
\]

なので\(\pi\)は次式で求められます。
\[
\pi = 4\left( \frac{1}{1} – \frac{1}{3} + \frac{1}{5} – \frac{1}{7} + \frac{1}{9} + … \right) = \sum_{k=1}^{\infty}\frac{(-1)^{k-1}}{2k-1}
\tag{7}
\]

ライプニッツの級数は,単純で美しい形で大好きなんですが,重大な欠点があります。

この公式は単純な形をしているが実際の円周率の計算に用いるには収束が非常に遅いために全く適していない。10進法での正確な値 (= 3.1415926535…) を10桁分計算するだけで100億回以上の計算を要するほどである。

Wikipediaより引用

収束おっそ!

 

 

 

Leave a Comment

メールアドレスが公開されることはありません。