気象情報解析特論　第1回　フーリエ級数

担当教員: 神山 翼 (tsubasa/at/is.ocha.ac.jp, @t_kohyama, 理学部3号館703号室)

In [1]:
import numpy as np
from scipy.special import factorial
import matplotlib.pyplot as plt

マクローリン級数

まず，フーリエ級数に類似する例として，マクローリン級数について復習しましょう。

マクローリン級数を用いると，（無限回微分可能で，かつ剰余項が0に収束するような）関数$f(x)$を$x=0$のまわりで多項式でフィッティングすることができます。

$$f(x) = f(0) + f'(0)x + \frac{f''(0)}{2!}x^2 + \frac{f'''(0)}{3!}x^3 + ...$$

たとえば，指数関数$e^x$を多項式でフィッティングすることを考えましょう。

$$e^x = 1 + x + \frac{1}{2!}x^2 + \frac{1}{3!}x^3 + ...$$
In [2]:
# xを定義
x = np.arange(-3, 3, 0.01)

# numpyが計算してくれるe^xをいまは厳密解と思うことにする
exp_analytical = np.exp(x)

# e^xの近似解（マクローリン級数5項目まで）
exp_approx = np.zeros_like(x)
for ii in range(0, 5):
    exp_approx = exp_approx + (1/factorial(ii))*x**(np.real(ii))

厳密解と近似解を描画してみます。

In [3]:
plt.plot(x, exp_analytical, 'k') # 厳密解を黒
plt.plot(x, exp_approx, 'b') # 近似解を青
plt.show()

$x=0$のまわりで，指数関数がうまく近似できていることがわかります。

一般に，マクローリン級数を

$$ \begin{eqnarray} f(x) &=& c_0 + c_1 x + c_2 x^2 + c_3 x^3 + ...\\ &=& \sum_{n=0}^\infty c_n x^n \end{eqnarray} $$

と書くことにすると，マクローリン級数の係数は

$$c_n = \frac{f^{(n)}(0)}{n!}$$

となります。$e^x$については$c_n = 1/n!$ となり，$c_n$のグラフを書くと次のようになります。

In [4]:
n = np.arange(0, 6)
plt.plot(n, 1/factorial(n), 'r')
plt.scatter(n, 1/factorial(n))
plt.show()

このグラフを見ると，多項式の次数が上がるほど，その項の重みが小さくなることがわかります。

これを少し難しい言葉で言い直すと，指数関数をべき関数$x^n$を重み付きで重ね合わせたもの だと思ったとき，低次の項に由来する成分が多く含まれていることを意味します。

フーリエ級数

指数関数では，多項式フィッティングがうまく機能しました。実は，繰り返し周期をもつ関数の場合，別の面白い方法があります。

一般に，周期$2L$を持つ（区分的に滑らかな）周期関数$f(x)$をフィッティングしようと思った場合，べき関数$x^n$の重ね合わせではなく，色々な周期を持った三角関数，つまり$\cos(\frac{n\pi x}{L})$および$\sin(\frac{n\pi x}{L})$を重み付きで重ね合わせたものだと思う方法があります。

$$ \begin{eqnarray} f(x) &=& \frac{a_0}{2} + \left\{a_1 \cos(\frac{\pi x}{L}) + b_1 \sin(\frac{\pi x}{L})\right \} + \left\{a_2 \cos(\frac{2 \pi x}{L}) + b_2 \sin(\frac{2 \pi x}{L})\right \} + ... \\ &=& \frac{a_0}{2} + \sum_{n=1}^\infty \left\{a_n \cos(\frac{n \pi x}{L}) + b_n \sin(\frac{n \pi x}{L})\right \} \end{eqnarray}$$

これをフーリエ級数といい，フーリエ級数の係数$a_n$と$b_n$は次のように表すことができます。

$$a_0 = \frac{1}{L} \int_0^{2L} f(x)dx$$$$a_n = \frac{1}{L} \int_0^{2L} f(x)\cos(\frac{n\pi x}{L})dx$$$$b_n = \frac{1}{L} \int_0^{2L} f(x)\sin(\frac{n\pi x}{L})dx$$

データ解析の視点からいうと，これらは$f(x)$のcosやsinに対する回帰係数そのものです（復習：環境情報論第8回）。

発展的な注：これらの公式は，f(x)と三角関数を各$x$で掛けた後それらを足し合わせること，すなわち「内積を取る」ことに対応しています。各係数は，「cosやsinへの射影の大きさ」であるということです。回帰係数を計算する際に必要な「共分散」も結局は内積だったことを思い出しましょう（環境情報論第11回を参照）。三角関数は互いに直交しているゆえ，重回帰を考える必要はありません。

いま，以下のような$L = 1$の矩形波をフーリエ級数展開することを考えます。

$$f(x) = \left\{ \begin{array}{ll} -1 & (2m-1 \leq x \leq 2m)\\ 1 & (2m \leq x \leq 2m+1) \end{array} \right. \ \ \ ただしmは整数$$
In [5]:
square_wave = np.ones_like(x)
for ii in range(0, x.shape[0]):
        if int((x[ii] >0) * (int(x[ii])%2 == 1) + (x[ii] <0) * (int(x[ii])%2 == 0)):
               square_wave[ii] = -1

plt.plot(x, square_wave, 'k')
plt.show()

矩形波のフーリエ級数を第$n$項まで求める関数を定義します。

In [6]:
def fourier_square(n):
    a = np.zeros(n+1)
    b = np.zeros(n+1)

    # a_0 = 0は明らかなので，a_n, b_n (n>0)を計算する
    for ii in range(1, n+1):
        # cosとsinの関数を定義
        cosine = np.cos(np.real(ii)*np.pi*x)
        sine = np.sin(np.real(ii)*np.pi*x)
        # cosとsinに対する回帰係数を求める（復習：Lecture 8)
        a[ii] = np.polyfit(cosine, square_wave, 1)[0] # np.polyfit(x, y, 1)の出力のうち，第ゼロ成分が回帰係数
        b[ii] = np.polyfit(sine, square_wave, 1)[0]
        
    # 矩形波の近似解（フーリエ級数第n項目まで）
    square_approx = np.zeros_like(x)
    for ii in range(0, n+1):
        square_approx = square_approx + a[ii]*np.cos(np.real(ii)*np.pi*x) + b[ii]*np.sin(np.real(ii)*np.pi*x)

    # 厳密解と近似解を描画してみる
    plt.plot(x, square_wave, 'k') # 厳密解を黒
    plt.plot(x, square_approx, 'b') # 近似解を青
    plt.show()
    
    # フーリエ係数を出力
    return a, b

試しに，第10項まで含めた近似でどこまでうまくいくか計算してみましょう。

In [7]:
[a, b] = fourier_square(10)

第10項まででも，なかなか良い近似になっていることがわかります。

面白半分で，これを第100項まで増やしてみましょう。

In [8]:
[a, b] = fourier_square(100)

あら不思議，ほぼ完璧な近似になっていることがわかります！

次に，マクローリン展開のときと同様に，$b_n$のグラフを書いておきましょう(実はこの例では，全ての$n$について$a_n$=0)。

In [9]:
n = np.arange(0, 15)
plt.plot(n, b[n], 'r')
plt.scatter(n, b[n])
plt.show()

このように，$n$が奇数のときのみ$b_n \neq 0$で，かつ$n$が大きくなればなるほど，項の重みが小さくなることがわかります。

繰り返しますが，フーリエ級数は色々な周期を持つ三角関数を重ね合わせて周期関数を近似する方法でした。

そのような視点でこのグラフを見ると，矩形波には長周期成分（=低周波成分ともいいます）がより多く含まれているということができます。

パワースペクトル

一般には，$a_n$や$b_n$は正の値とは限らず，負の値やゼロをとることもあります。そこで，三角関数の合成を用いることで

$$a_n \cos(\frac{n \pi x}{L}) + b_n \sin(\frac{n \pi x}{L}) = \sqrt{a_n^2 + b_n^2}\sin(\frac{n \pi x}{L} + \phi)$$

（ただし$\phi = \tan^{-1}\frac{a_n}{b_n}$）

と変形することで，三角関数の振幅の2乗$(a_n^2 + b_n^2)$を周波数(=周期の逆数。この例だと$n/2L$)ごとにプロットしたグラフを，パワースペクトルといいます。

In [10]:
plt.plot(n/2, a[n]**2 + b[n]**2, 'r') # いまはL = 1の矩形波について考えている
plt.xlabel("Frequency")
plt.ylabel("Power")
plt.show()

先ほどの周期が2（周波数が1/2，つまり$x=1$進むごとに1/2回分波が揺れている）の矩形波には，同じく周波数1/2を持つ正弦波の「パワー」（＝分散，あるいは波のエネルギーのこと）が最も多く含まれていることがわかります。次に多く含まれているのは周波数が3/2の波，その次は周波数が5/2の波です。

課題1（締切：次回の授業開始前まで）

この講義資料のようなNotebook形式で次の問に答え，html形式またはpdf形式に保存して提出して下さい。NotebookにはMarkdownセルを用いて，それぞれのセルについて（この講義資料のように）説明を加えて下さい。

A問題. fourier_square(n)を色々な$n$で試すことで，フーリエ級数の項数を増やすにつれ矩形波の近似が良くなっていくことを確かめてください。

B問題. 上の例では三角関数への回帰係数を数値的に求めることによって$a_n$や$b_n$を求めましたが，矩形波のように簡単な関数の場合は，$a_n$や$b_n$を手計算で求めることもできます。上記に示した積分の公式を用いて，$a_n$と$b_n$を計算してください。また，$b_n$がプログラムで求めた$b_n$に（数値誤差の範囲内で）一致することを確かめてください。

C問題. 上の例では矩形波のフーリエ級数を求めて描画しましたが，同様のことを三角波で行うことによって，fourier_triangle(n)を実装してください。

D問題. B問題と同様のことを，三角波についても行ってください。

提出方法：課題提出専用メールアドレスkohyama.coursework/at/gmail.com（普段のメールと違います！！/at/を@に変えてください）に対して，件名を「気象情報解析特論1，〇〇〇〇」として（〇〇〇〇にはご自身の氏名を入れてください），課題ファイルをメールに添付して送信してください。

課題に関する諸注意

Markdownセルは，pythonコードについての注釈を入れるためのセルです。Markdownセルでは，たとえば#を打つと見出しになったり，アスタリスクを入れると斜体や太字になったりします。詳しくは例えばこちら：https://qiita.com/tbpgr/items/989c6badefff69377da7

Markdownセルを入れるには，

・Anacondaの場合：セルにカーソルが合っている状態で，上のメニューバーで「Code」と書いてあるところを「Markdown」に変えます。

・Google Colabの場合：上のメニューバーで「+テキスト」をクリックします。

Notebook形式をhtml形式やpdf形式としてダウンロードするには，たとえば

・Anacondaの場合：File>Download as...>から，htmlまたはpdfを選択します。

・Google Colabの場合：「ファイル>印刷」を選択し，プリンターの設定画面が出たら左下の方のメニューで「プレビューで開く」を選び，プレビューでpdf形式に保存します（Safariで不具合が報告されているので，Google Chrome推奨）。

その他わからなければ，遠慮なく質問してください。質問は居室に来ていただくか，通常のメールアドレスtsubasa/at/is.ocha.ac.jpの方にいただけると助かります。

成績評価について

毎週課題を出します。ただし，全部の問題を解く必要はありません。

A問題 これが毎週きちんと解けていれば単位はあげます。

B問題 ここまでだいたいできていればA評価以上の成績はあげます。

C問題 Pythonを頑張りたい方，良い成績を狙いたい方，あるいは神山研の学生さんはどうぞ。きちんと挑戦した跡があれば，正誤は問いません。

D問題 DはD進のDです。神山研でD進しましょう。