1. はじめに
本記事ではフーリエ変換から高速フーリエ変換までまとめる。流れを掴むための記事であるため細かな部分は省略する場合があることを承知した上で読んでほしい。前提知識としては高校数学があることが望ましい。読み進めるうえで利用する知識は以下の通りである。
・区分級数展開(積分)
・複素平面
・数列(シグマ)
2. フーリエ変換
フーリエ変換を行う上で最も重要な知識はフーリエ級数展開である。フーリエ級数展開を読み解くことがフーリエ変換を理解するうえで最も近道となる。尚、フーリエ級数展開は以下の公式が知られている。
f(x)=a+\sum_{n=1} ^{\infty}(a_ncosnx+b_nsinnx)
ここで、$sinx$や$cosx$に掛かる係数が周波数になるのは周知の通りである。即ち、$a_n$と$b_n$の値を求めれば合成波を周波数毎に分解することが可能あり、フーリエ級数展開の目的はここにある。
さて、フーリエ級数展開の目的を知ったところで実際に$a_n$と$b_n$の値を求めるのだが、その前に三角関数の直行性の話をする。直行性の話には積和の公式を用いる。
sinnxsinmx=\frac{1}{2}{cos(n+m)x-cos(n-m)x} \\
cosnxcosmx=\frac{1}{2}{cos(n+m)x+cos(n-m)x} \\
sinnxcosmx=-\frac{1}{2}{sin(n+m)x+sin(n-m)x}
三角関数の直行性とは、n, mの値が異なるときの積分が0になることである。
\begin{align}
\int_{-\pi}^{\pi} sinnxsinmxdx&=\frac{1}{2} \int_{-\pi}^{\pi} {[cos(n+m)x+cos(n-m)x]}dx \\
&=\frac{1}{2}[\frac{1}{n+m}sin(n+m)x-\frac{1}{n-m}sin(n-m)x]_{-\pi}^{\pi} \\
&=0
\end{align}
ここで、$n=m$１のとき、
\begin{align}
\int_{-\pi}^{\pi} sinnxsinnxdx&=\frac{1}{2} \int_{-\pi}^{\pi} {[cos(n+n)x+cos(n-n)x]}dx \\
&=\frac{1}{2}[\frac{1}{2n}sin(2n)x+x]_{-\pi}^{\pi} \\
&=\pi
\end{align}
同様に、
\begin{align}
\int_{-\pi}^{\pi} cosnxcosmxdx&=\frac{1}{2} \int_{-\pi}^{\pi} {[cos(n+m)x+cos(n-m)x]}dx \\
&=0 \\
\int_{-\pi}^{\pi} cosn^{2}xdx&=\frac{1}{2} \int_{-\pi}^{\pi} [{cos(2n)x}+1]dx \\
&=\pi
\end{align}
$sinnx$は奇関数、$cosmx$は偶関数であるため積分の結果が0になるのは明白である。これらを用いて実際に係数を求める。
\begin{align}
\int_{-\pi}^{\pi} f(x)cosmxdx&=\int_{-\pi}^{\pi}acosmxdx+\sum_{n=1} ^{\infty}\int_{-\pi}^{\pi}(a_ncosnx+b_nsinnx)cosmxdx \\
&=a_m\pi ・・・①
\end{align}
\begin{align}
\int_{-\pi}^{\pi} f(x)sinmxdx&=\int_{-\pi}^{\pi}asinmxdx+\sum_{n=1} ^{\infty}\int_{-\pi}^{\pi}(a_ncosnx+b_nsinnx)sinmxdx \\
&=b_m\pi ・・・②
\end{align}
\begin{align}
\int_{-\pi}^{\pi} f(x)dx&=\int_{-\pi}^{\pi}adx \\
&=2\pi a ・・・③
\end{align}
①. ②. ③より、
\begin{align}
a_m&=\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)cosmxdx \\
b_m&=\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)sinmxdx \\
a&=\frac{1}{2\pi}\int_{-\pi}^{\pi}f(x)dx
\end{align}
$a_m$の式に$m=0$を代入する。
\begin{align}
a_0&=\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)cos0dx \\
&=\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)dx
\end{align}
これは$a$の式の2倍になっていることが伺える。よって以降$a$を$\frac{a_0}{2}$と書く。これらをまとめると、
\frac{a_0}{2}+\sum_{n=1} ^{\infty}(a_ncosnx+b_nsinnx) \\
\begin{pmatrix}
a_n&=\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)cosnxdx \\
b_n&=\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)sinnxdx
\end{pmatrix}
以上が周期$2\pi$における実数のフーリエ級数である。
次に複素数に対応させる。数学で最も美しいと言われているオイラーの公式はご存じだろうか。
e^{i\theta}=cos\theta+isin\theta
尚、$\theta= nx$とすると
\begin{align}
e^{inx}&=cosnx+isinnx \\
e^{-inx}&=cosnx-isinnx
\end{align}
であるから2つの式を$sinnx$と$cosnx$について解くと
\begin{align}
sinnx&=\frac{1}{2i}(e^{inx}-e^{-inx}) \\
cosnx&=\frac{1}{2}(e^{inx}+e^{-inx})
\end{align}
となる。この$sinnx$と$cosnx$をフーリエ級数展開に代入する。
\frac{a_0}{2}+\sum_{n=1} ^{\infty} [\frac{a_n}{2}(e^{inx}+e^{-inx})+\frac{b_n}{2i}(e^{inx}-e^{-inx})] \\
\begin{align}
&=\frac{a_0}{2}+\sum_{n=1} ^{\infty} [(\frac{a_n}{2}+\frac{b_n}{2i})e^{inx}+(\frac{a_n}{2}-\frac{b_n}{2i})e^{-inx}] \\
&=\frac{a_0}{2}+\sum_{n=1} ^{\infty} (\frac{a_n-ib_n}{2}e^{inx}+\frac{a_n+ib_n}{2}e^{-inx})
\end{align}
同様に$a_n$と$b_n$の式にも代入すると、
\begin{align}
a_n&=\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)\frac{1}{2}(e^{inx}+e^{-inx})dx \\
&=\frac{1}{2 \pi}\int_{-\pi}^{\pi}f(x)(e^{inx}+e^{-inx})dx
\end{align}
\begin{align}
b_n&=\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)\frac{1}{2i}(e^{inx}-e^{-inx})dx \\
&=\frac{1}{2i\pi}\int_{-\pi}^{\pi}f(x)(e^{inx}-e^{-inx})dx
\end{align}
ここで、
\begin{align}
a_n-ib_n&=\frac{1}{2 \pi}\int_{-\pi}^{\pi}f(x)(e^{inx}+e^{-inx})dx-\frac{1}{2\pi}\int_{-\pi}^{\pi}f(x)(e^{inx}-e^{-inx})dx \\
&=\frac{1}{2 \pi}\int_{-\pi}^{\pi}[f(x)(e^{inx}+e^{-inx})-f(x)(e^{inx}-e^{-inx}]dx \\
&=\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)e^{-inx}dx ・・・① \\
\end{align}
\begin{align}
a_n+ib_n&=\frac{1}{2 \pi}\int_{-\pi}^{\pi}f(x)(e^{inx}+e^{-inx})dx+\frac{1}{2\pi}\int_{-\pi}^{\pi}f(x)(e^{inx}-e^{-inx})dx \\
&=\frac{1}{2 \pi}\int_{-\pi}^{\pi}[f(x)(e^{inx}+e^{-inx})+f(x)(e^{inx}-e^{-inx}]dx \\
&=\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)e^{inx}dx \\
&=(a_n-ib_n)^{-1} ・・・② \\
\end{align}
\begin{align}
a_0&=\frac{1}{2 \pi}\int_{-\pi}^{\pi}f(x)(e^{0}+e^{0})dx \\
&=\frac{1}{2 \pi}\int_{-\pi}^{\pi}2f(x)dx \\
&=\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)dx \\
&=(a_n-ib_n)^{0} ・・・③
\end{align}
①. ②. ③より、$\frac{a_n-ib_n}{2}=C_n$とすると、
$$C_0+\sum_{n=1} ^{\infty} (C_ne^{inx}+C_n^{-1}e^{-inx})=\sum_{n=-\infty} ^{\infty}C_ne^{inx}$$
となる。よって、周期$2\pi$の複素フーリエ級数は以下の通りである。
\sum_{n=-\infty} ^{\infty}C_ne^{inx} \\
(C_n=\frac{1}{2\pi}\int_{-\pi}^{\pi}f(x)e^{-inx}dx)
これまで周期が$2\pi$の場合のみ考えてきた。しかし、実際に扱うデータの周期が$2\pi$であるとは限らない。よって周期を一般化する必要がある。
下の画像は$sinx$と$sin2x$のグラフである。見てわかる通り、$sinnx$は$sinx$の周期が$\frac{1}{n}$倍になる。即ち、周期$2\pi$の関数を$\frac{2\pi}{L}$倍するれば周期Lとなる。
ここで、$e^{i \theta}=cosi \theta+sini \theta$であり周期がLであるから
\sum_{n=-\infty} ^{\infty}C_ne^{i\frac{2\pi n}{L}x} \\
(C_n=\frac{1}{L}\int_{-\frac{L}{2}}^{\frac{L}{2}}f(x)e^{-i\frac{2\pi n}{L}x}dx)
さて、ようやく本番のフーリエ変換を実装する。と言っても上記で求めた一般周期のフーリエ級数を少し改造してあげればよいので身構える必要はない。一般周期のフーリエ級数を求める際に、「実際に扱うデータの周期が$2\pi$であるとは限らない。」と話をしたが、世の中のデータのと殆どはそもそも周期をもたない。フーリエ変換は周期をもたないデータに対してフーリエ級数展開を行うことである。
周期をもたないということは、周期が$\infty$と捉えることができる。つまり、Lを$\infty$に近似すればよい。
まず始めに、$\frac{2\pi n}{L}=\omega$とし、$\omega$の変化量を求める。
\begin{align}
\delta \omega &=\omega_{n+1}-\omega_n \\
&=\frac{2\pi (n+1)}{L}-\frac{2\pi n}{L} \\
&=\frac{2\pi n+2\pi}{L}-\frac{2\pi n}{L} \\
&=\frac{2\pi}{L}
\end{align}
フーリエ級数に$C_n$を代入すると、
\begin{align}
\sum_{n=-\infty} ^{\infty}[\frac{1}{L}\int_{-\frac{L}{2}}^{\frac{L}{2}}f(t)e^{-i\omega_nt}dt]e^{i\omega_nx}&=\sum_{n=-\infty} ^{\infty}[\frac{1}{2\pi}\frac{2\pi}{L}\int_{-\frac{L}{2}}^{\frac{L}{2}}f(t)e^{-i\omega_nt}dt]e^{i\omega_nx} \\
&=\frac{2\pi}{L}\sum_{n=-\infty} ^{\infty}[\frac{1}{2\pi}\int_{-\frac{L}{2}}^{\frac{L}{2}}f(t)e^{-i\omega_nt}dt]e^{i\omega_nx} \\
&=\delta \omega \sum_{n=-\infty} ^{\infty}[\frac{1}{2\pi}\int_{-\frac{L}{2}}^{\frac{L}{2}}f(t)e^{-i\omega_nt}dt]e^{i\omega_nx}
\end{align}
ここで、Lを$\infty$に近似するので$\lim_{L \to \infty}\frac{2\pi}{L}=0$となり$\delta \omega$は0に収束する。
$\delta \omega$が0に近づくのであれば、$\delta \omega \sum_{n=-\infty} ^{\infty}[\frac{1}{2\pi}\int_{-\frac{L}{2}}^{\frac{L}{2}}f(t)e^{-i\omega t}dt]e^{i\omega x}$は下記画像のような区分求積法になる。よって、
\lim_{L \to \infty}\frac{2\pi}{L}\sum_{n=-\infty} ^{\infty}[\frac{1}{2\pi}\int_{-\frac{L}{2}}^{\frac{L}{2}}f(t)e^{-i\omega t}dt]e^{i\omega x}=\int_{-\infty}^{\infty}[\frac{1}{2\pi}\int_{-\infty}^{\infty}f(t)e^{-i\omega t}dt]e^{i\omega x}d\omega
となり、
\int_{-\infty}^{\infty}F(x)e^{i\omega}d\omega \\
(F(x)=\frac{1}{2\pi}\int_{-\infty}^{\infty}f(x)e^{-i\omega}dx)
ここまで、フーリエ級数展開が元の関数に収束する前提で話してきたが、実際には収束する条件がある。その条件とは、$f(x)$が区分的に滑らかであることだ。区分的に滑らかとは、簡単に言うと$x$軸方向に連続であるか否かである。$f(x)$が連続である必要はないため、ガウス関数も区分的に滑らかであると言える。
さらに、フーリエ変換では積分範囲が無限大になっている。無限大での積分は絶対化積分である必要がある。絶対化積分とは、$|f(x)|<\infty$において積分が可能であることを指す。
これらを踏まえて、$f(x)$が$(-\infty, \infty)$で区分的に滑らかかつ$(-\infty, \infty)$において絶対化積分であるとき以下のフーリエの積分公式が成り立つ。
f(x)～\int_{-\infty}^{\infty}F(x)e^{i\omega}d\omega ・・・逆フーリエ変換 \\
F(x)=\frac{1}{2\pi}\int_{-\infty}^{\infty}f(x)e^{-i\omega}dx ・・・フーリエ変換
3. 離散フーリエ変換
さて、2章でフーリエ変換について扱ってきたが、2章の章末で述べた通りフーリエ変換を行うためには区分的に滑らかである必要がある。しかし、現代扱うデータはサンプリングレート毎に収集された離散的なデータである。このようなデータにはフーリエ変換を行うことができない。離散フーリエ変換ではこういった離散データに対応できるようフーリエ変換を改造する。
離散データに対応させる都合上、一般周期$L$に対するフーリエ変換で考えるとわかりやすい。なぜならデータ数は有限個であるため、$N$個のデータに対して周期が$N$であると考えることが可能だからだ。周期$N$に対するフーリエ級数展開は以下の通りである。
f(x)=\sum_{n=-\infty} ^{\infty}C_ne^{i\frac{2\pi n}{N}x} \\
C_n=\frac{1}{N}\int_{-\frac{N}{2}}^{\frac{N}{2}}f(x)e^{-i\frac{2\pi n}{N}x}dx
上記した通り、データ数を周期と考える。データのはじめを0番目とすると以下のように書き換えることができる。
f(x)=\sum_{n=0} ^{N-1}C_ne^{i\frac{2\pi n}{N}x} \\
C_n=\frac{1}{N}\int_{0}^{N-1}f(x)e^{-i\frac{2\pi n}{N}x}dx
ここで、データは1番目、2番目、3番目...のように1ずつ変化する。よって、x軸方向への変化量$(dx)$は1であるため、
\begin{align}
C_n&=\frac{1}{N}\int_{0}^{N-1}f(x)e^{-i\frac{2\pi n}{N}x}dx \\
&=\frac{1}{N}\sum_{n=0}^{N-1}f(n)e^{-i\frac{2\pi n}{N}x}
\end{align}
$e^{-i\frac{2\pi n}{N}x}=(e^{-i\frac{2\pi}{N}})^{nx}$であるから$e^{-i\frac{2\pi n}{N}x}=\omega^{nx}$とする。ほとんどの参考書が$C_n$を$F(x)$、と表記しているため統一して、
f(x)=\sum_{n=0} ^{N-1}F(x)\omega^{-nx} ・・・逆離散フーリエ変換 \\
F(x)=\frac{1}{N}\sum_{n=0}^{N-1}f(n)\omega^{nx} ・・・離散フーリエ変換
4. 高速フーリエ変換
本章では高速フーリエ変換について扱う。離散フーリエ変換ではx番目のデータに対してすべてのデータと計算した和を求めている。つまり計算数はデータ数の2乗になる。高速フーリエ変換ではこの離散フーリエ変換の計算回数を減らして高速化するアルゴリズムである。結論から離すと、データ数が2のべき乗である必要があるが、計算量が$\frac{N}{2}\log_2N$回となる。
ほとんどの参考書では回転因子という考えをもとに、ビット反転とバタフライ演算という手法を解説している。だが、今回はそういった類のものの説明はしない(アルゴリズム上は利用している)。
まずはじめに、高速フーリエ変換について計算しやすくするために式を変形する。離散フーリエ変換を逆離散フーリエ変換に代入すると、
\begin{align}
f(x)&=\sum_{n=0} ^{N-1}\frac{1}{N}\sum_{n=0}^{N-1}f(n)\omega^{nx}\omega^{-nx} \\
&=\frac{1}{N}\sum_{n=0} ^{N-1}\sum_{n=0}^{N-1}f(n)\omega^{nx}\omega^{-nx}
\end{align}
と表記でき、また分解すると
f(x)=\frac{1}{N}\sum_{n=0} ^{N-1}F(x)\omega^{-nx} \\
F(x)=\sum_{n=0}^{N-1}f(n)\omega^{nx}
つまり、$\frac{1}{N}$は離散フーリエ変換、逆離散フーリエ変換どちらについていても構わない。なんなら両方に$\frac{1}{\sqrt N}$を掛けていても計算結果は同じである。
次に、離散フーリエ変換について偶数番目と奇数番目に分けて計算するとある法則が出てくる。
\begin{align}
F(x)&=\sum_{n=0}^{N-1}f(n)\omega^{nx} \\
&=\sum_{n=0}^{\frac{N}{2}-1}f(2n)\omega^{2nx}+\sum_{n=0}^{\frac{N}{2}-1}f(2n+1)\omega^{(2n+1)x} \\
&=\sum_{n=0}^{\frac{N}{2}-1}f(2n)\omega^{2nx}+\sum_{n=0}^{\frac{N}{2}-1}f(2n+1)\omega^{2nx}\omega^{x} \\
&=\sum_{n=0}^{\frac{N}{2}-1}f(2n)\omega^{2nx}+\omega^{x}\sum_{n=0}^{\frac{N}{2}-1}f(2n+1)\omega^{2nx}
\end{align}
ここで、$\omega$は$e^{-i\frac{2\pi}{N}}$であるということを思い出してほしい。
\begin{align}
(e^{-i\frac{2\pi}{N}})^{2nx}&=e^{-i\frac{2\pi}{N}2nx} \\
&=e^{-i\frac{2\pi n}{N/2}x}
\end{align}
さて、法則に気づいただろうか。離散フーリエ変換を偶数番目と奇数番目に分けると偶数番目の離散フーリエ変換と、奇数番目の離散フーリエ変換に$\omega nx$を掛けたものの和になっている。これが上記した法則である。そしてもちろん、偶数番目をさらに偶数番目と奇数番目、奇数番目をさらに偶数番目と奇数番目、と分割して計算が可能である。つまり、離散フーリエ変換は再帰関数である。本章冒頭でデータ数が2のべき乗である必要があるといったのは、データ数が2個になるまで2分割し続けるからである。よって$\frac{N}{2}$回の計算を$\log_2N$回繰り返すので、計算回数が$\frac{N}{2}\log_2N$回となる。これが高速フーリエ計算のアルゴリズムの一つである。
5. おわりに
離散フーリエ変換と高速フーリエ変換は行列として計算を行うことができる。そういった場合に4章で名前だけ出した回転因子などの考え方が必要になる。気になる方は調べてみてほしい。
それでは、長くなってしまったが以上で終わりにする。最後まで読んでくれた人はありがとう。