フーリエ変換を理解するまでの数学的メモ
理解が難しく敬遠していたフーリエ変換だが、この機会に挑戦してみることにした。参考書は次の本。この手の本は読み手のくじける心を非常に分かっているのでおすすめする。
またその後に、同じ作者（渋谷道雄氏）が書いたこちらの本も読むと、より理解しやすい。 実際にプログラミングでフーリエ変換を体験してみたい方にはおすすめ。
フーリエ変換は大学で学ぶ数学だが、実は高校生の数学レベルでもなんとかなる。「三角関数と積分、微分」だけでほとんど理解できるのだ。
ここではあくまで自分的なメモとして、これらの書籍のおさらいを残しておく。PythonでFFTを使った実践的なプログラミングをやりたい方は、 Pythonではじめてのフーリエ変換・逆フーリエ変換 を参考に。
フーリエ変換を理解するための、三角関数
半径1の円を単位円とする。単位円において円周の長さが1の時の角度を、1ラジアンと表す。1円周 = 2πrの公式より、単位円における1円周は2πとなるので、ラジアンで表せば1円周 = 2π [rad]である。
sinは正弦、cosは余弦と呼ぶ。cosθはsinθがπ/2だけずれた関係にあり、直交している。直行とは、お互いを掛け合わせた関数の積分の結果がゼロになることである。
x=cosθ, y=sinθ でθを変数にした媒介変数表示で表現できる。
ピタゴラス三平方の定理 \(x^2 + y^2 = r^2\) より
$$cos^2θ + sin^2θ = 1$$
円運動をラジアンと時間で表したものを角速度と呼び、ω[rad/s]で表す。たとえば単位円における1秒間に1周する円運動は、2π [rad/s]と表現されるだろう。これは1Hzと言い換えることができる。
では角速度が2倍の4π [rad/s]だったら何Hzになるだろうか？一秒間に2πの円周を2周回るから2Hzと言える。
これにより角速度は角周波数とも言える。
三角関数の加法定理
$$ \sin{(α+β)} = \sin{α}\cos{β} + \cos{α}\sin{β} \tag{1}$$ $$ \cos{(α+β)} = \cos{α}\cos{β} - \sin{α}\sin{β} \tag{3}$$ $$ \sin{α}\cos{β} = \frac{1}{2}\\{\sin{(α + β)} + \sin{(α - β)}\\} \tag{5}$$ $$ \cos{α}\cos{β} = \frac{1}{2}\\{\cos{(α + β)} + \cos{(α - β)}\\} \tag{7}$$上の加法定理により次の積和の公式が導き出される。たとえば(1)+(2)より(5)となる。
$$A = α + β \tag{8}$$ $$ α = \frac{A + B}{2} \tag{10}$$ $$\sin{α}+\sin{β}=2\sin{\frac{α+β}{2}}\cos{\frac{α-β}{2}} \tag{12}$$さらに
$$\cos{α}+\cos{β}=2\cos{\frac{α+β}{2}}\cos{\frac{α-β}{2}} \tag{14}$$ $$(\sin{x})^{\prime} = \cos{x}\tag{16}$$とすると次のように展開できる。
$$\int \sin{x}dx = -\cos{x}\tag{18}$$ $$\int_{0}^{2π} \sin{x}\cos{x}dx\tag{20}$$これらを(1)-(2)に代入することで和積の公式を導くことができる。
$$\int_{0}^{2π} \sin{x}\sin{nx}dx\tag{21}$$ $$\int_{0}^{2π} \sin{x}\cos{nx}dx\tag{22}$$ $$\int_{0}^{2π} \sin^{2}{x}dx = π\tag{23}$$ $$ f(x) = a {\cos x} + b {\sin x} \tag{22}$$三角関数の微分積分
$$r=\sqrt{a^{2}+b^{2}}\tag{23}$$ $$F(x)=\frac{1}{2}a_0+a_1\cos{x}+a_2\cos{2x}+a_3\cos{3x}+ ... + a_n\cos{nx}+ ...$$ $$F(x)=\frac{1}{2}a_0 + \sum_{n=1}^{\infty}(a_n{\cos nx} + b_n{\sin nx})\tag{25}$$ $$f(t) = \sum_{n=1}^{\infty}\frac{1}{n}{\sin nt}\tag{26}$$たとえばsinxを積分してみよう。
0からπの範囲では答えが2となるが、0から2πの範囲では0となる。これはグラフをみても理解できる。
直交する関数の積の定積分
直交とは2つの関数の積の定積分がゼロになる関数である。
たとえば次の式を計算してみる。
$$f(t) = \sum_{n=1}^{\infty}\frac{1}{(2n-1)}{\sin (2n-1)t }\tag{27}$$
答えはゼロとなったのでsinxとcosxは直交していると言える。 同様に、次の式も直交する。
$$f(t) = \sum_{n=1}^{\infty}{-1}^{(n+1)}\frac{1}{(2n-1)^{2}}{\sin (2n-1)t }\tag{28}$$ $$\int_{0}^{2π} \cos^{2}{nx}dx = π\tag{23}$$
しかし、sinxとsinxの積の定積分はゼロにはならない。
$$a_n = \frac{1}{π}\int_{0}^{2π}F(x)\cos{nx}dx\tag{24}$$これらの性質はとても便利で、のちのフーリエ変換におけるフィルタリングの役割となる。
直交する関数の合成と位相差
$$b_n = \frac{1}{π}\int_{0}^{2π}F(x)\sin{nx}dx\tag{25}$$a=1, b=1 の時と、a=-1, b=1の時を比較してみる。
2つの関数の波形は振幅、形がまったく同じだが、位相だけがずれている。aとbの値を変えることでいろいろな位相ずれを表現することができる。\({\sin(x + θ)}\)でも位相を変えることができるが、sinとcosの合成関数の方が都合が良い。
また、ベクトルととらえればピタゴラスの定理より、直交関数の合成の大きさrは
$$a_0 = \frac{1}{2π}\int_{0}^{2π}F(x)dx\tag{26}$$となる。
これらにより、(22)の式は次の周波数スペクトルグラフを書くことができる。
フーリエ級数展開
ある特定の周波数を\( f(x) = a\cos{nx} + b\sin{nx}\)で表現できるならば、どんな複雑な周波数もそれらの合成、和であると考えることができる。
$$r_n=\sqrt{a_n^{2}+b_n^{2}}\tag{23}$$$$+b_1\sin{x}+b_2\sin{2x}+b_3\sin{3x}+ ... + b_n\sin{nx}+ ... \tag{24}$$
これにより次のフーリエ級数展開の公式が導き出される。
$$F(x)=\frac{1}{2}a_0 + \sum_{n=1}^{\infty}(a_n{\cos nx} + b_n{\sin nx})\tag{25}$$
\( \frac{1}{2}a_0 \) は波形全体を上下に移動させる役目を持っている。
Pythonでいろいろな波形の合成
ちょっと休憩、いろいろな波形をPythonでプログラミングして遊んでみた。 最大値を1にする正規化は行っていないが、いろいろな波形がsinとcosで作れることが分かった。
Sawtooth wave
$$f(t) = \sum_{n=1}^{\infty}\frac{1}{n}{\sin nt}\tag{26}$$
def sawtooth():
t = np.arange(0, 4 * np.pi, 0.01)
wave = 0
for n in range(1, 40):
a = 1 / n
wave += a * np.sin(n * t)
graph.graph(x=t, y=wave,
title='Sawtooth wave')Square wave
$$f(t) = \sum_{n=1}^{\infty}\frac{1}{(2n-1)}{\sin (2n-1)t }\tag{27}$$
2n-1は奇数を表す。
def square():
t = np.arange(0, 4 * np.pi, 0.01)
wave = 0
for n in range(1, 40):
if n % 2 == 1:
a = 1 / n
wave += a * np.sin(n * t)
graph.graph(x=t, y=wave,
title='Square wave')Triangle wave
$$f(t) = \sum_{n=1}^{\infty}{-1}^{(n+1)}\frac{1}{(2n-1)^{2}}{\sin (2n-1)t }\tag{28}$$
矩形波よりも高音へ向かう減衰率が大きいのが特徴。\({-1}^{(n+1)}\)は1つ置きに和減を変更するための表現。
def triangle():
t = np.arange(0, 4 * np.pi, 0.01)
wave = 0
i = 0
for n in range(1, 40):
if n % 2 == 1:
a = 1 / (n * n)
if i % 2 == 0:
wave += a * np.sin(n * t)
else:
wave -= a * np.sin(n * t)
i += 1
graph.graph(x=t, y=wave,
title='Triangle wave')フーリエ解析
元の波形（関数）から各周波数がどのくらいの大きさで含んでいるかを求めるのがフーリエ解析である。
周期が不明な非周期関数の場合はたとえば１秒などで区切る。そして一番低い周波数から計算可能な一番高い周波数まで、すべての周波数成分を1つづつ考えていく。周波数成分を抽出するにはフィルターを使うことになるが、数学的には三角関数の直交の性質を利用してフィルタリングする。
たとえばフーリエ級数の式(24)(25)から\(a_n\cos{nx}\)だけを残したい場合は、F(x)全体にcosnxを掛けて積分すれば良いことになる。nの値が等しいcosどうしやsinどうしは直交しないのだった。逆にそれ以外は直交することになるので積分するとゼロとなってしまい、結果としてcosnxが残ることになる。
ところでcosnx * cosnxの積分は、
$$\int_{0}^{2π} \cos^{2}{nx}dx = π\tag{23}$$
となる。
図のようにπを区切りにひっくり返して収めることで、横幅がπで高さが1の長方形の面積と考ることができる。式(23)は\(a_n\)が1なので、\(1 * π = π\)と解釈する。
これにより求める式である積分の式をπで割ることでフーリエ係数\(a_n, b_n\)を求めることができる。
$$a_n = \frac{1}{π}\int_{0}^{2π}F(x)\cos{nx}dx\tag{24}$$
$$b_n = \frac{1}{π}\int_{0}^{2π}F(x)\sin{nx}dx\tag{25}$$
また\(a_0\)を求めるには、F(x)のほとんどはsinやcosで出来ているのでsinやcosを積分しても0であることを利用する。
$$a_0 = \frac{1}{2π}\int_{0}^{2π}F(x)dx\tag{26}$$
最後に周波数成分を求めるにはフーリエ係数より、ピタゴラスの定理（三平方の定理）をつかって次のようになる。
$$r_n=\sqrt{a_n^{2}+b_n^{2}}\tag{23}$$
以上、フーリエ変換を使って次のスペクトラムは作成されている。ただし、実際にはもっと高速なFFTアルゴリズムが使われる。
