高速フーリエ変換・ラプラス変換・Z変換

事前に読む

/icons/hr.icon

高速フーリエ変換

離散フーリエ変換$ F(k) = \sum_{n = 0}^{N - 1}f[n]e^{-j\frac{2\pi}{N}kn} (k = 0, 1, \cdots, N - 1) を計算する時、一つの$ kについて係数を計算するのに$ N回の和を取るので、全ての$ kについて係数を計算するためには計算量が$ O(N^2)となり計算量が多い。

また、一般に計算機では積の計算は負荷の高い処理であり、可能な限り積の計算を省略したい。

そこで、高速フーリエ変換をおこなうアルゴリズムの一つである（基数2の）Cooley-Tukeyのアルゴリズム（クーリー - テューキーのアルゴリズム）を考える。

なお、今回は周波数を偶奇で分けることを考えるが、この様に周波数に対して分割をおこなうことを周波数間引きといい、逆に時間に対して分割をおこなうことを時間間引きという。

まずクーリー - チューキーのアルゴリズムでは、計算量や計算負荷を削減するために、複素指数関数の特性である周期性と対称性を利用する。

周期性とは複素平面において角度が一周すると一周前と同じ値を取る特性であり、対称性とは複素平面上の単位円の中心に対して点対称の位置においても元の点と符号だけが異なる同じ値を取る特性である。

https://scrapbox.io/files/5fd51a061067f50022df01b5.png $ N = 8の複素平面

ここで、関数$ f[n] が一周期内で取る値の数$ Nを$ 2の冪乗であるとする。

離散フーリエ変換の複素指数関数部を

$ W_N = e^{-j\frac{2\pi}{N}}\cdots(1)

としたとき、この$ (1)式を回転因子（回転子）という。

このとき、離散フーリエ変換を書き直すと、

$ F[k] = \sum_{n = 0}^{N - 1}f[n]W_N^{nk}

ここで、回転因子の対称性を利用するためにサンプリング点$ kを偶数要素と奇数要素に分割すると、

$ \begin{aligned} F[2k] &= \sum_{n=0}^{N-1}f[n]W_N^{2nk} \\ &= \sum_{n=0}^{\frac{N}{2}-1}f[n]W_N^{2nk} + \sum_{n=0}^{\frac{N}{2}-1}f[n + \frac{N}{2}]W_N^{2(n + \frac{N}{2})k} \\ &= \sum_{n=0}^{\frac{N}{2}-1}f[n]W_{\frac{N}{2}}^{nk} + \sum_{n=0}^{\frac{N}{2}-1}f[n + \frac{N}{2}]W_{\frac{N}{2}}^{(n + \frac{N}{2})k} \quad (\because W_N^{2k} = W_{\frac{N}{2}}^k) \\ &= \sum_{n=0}^{\frac{N}{2}-1}(f[n] + f[n + \frac{N}{2}])W_{\frac{N}{2}}^{nk} \quad (\because W_{\frac{N}{2}}^{n+\frac{N}{2}} = W_{\frac{N}{2}}^{n}) \end{aligned}

$ \begin{aligned} F[2k + 1] &= \sum_{n=0}^{N-1}f[n]W_N^{n(2k+1)} \\ &= \sum_{n=0}^{\frac{N}{2}-1}f[n]W_N^{n(2k+1)} + \sum_{n=0}^{\frac{N}{2}-1}f[n + \frac{N}{2}]W_N^{(n + \frac{N}{2})(2k+1)} \\ &= \sum_{n=0}^{\frac{N}{2}-1}(f[n] - f[n + \frac{N}{2}])W_{\frac{N}{2}}^{nk}W_N^n \end{aligned}

が得られる。

これは、出力であるサンプリングの点$ k の偶数要素は要素数$ \frac{N}{2} の$ f[n] + f[n + \frac{N}{2}] の離散フーリエ変換に等しく、奇数要素は要素数$ \frac{N}{2} の$ f[n] - f[n + \frac{N}{2}]W_N^n のフーリエ変換に等しいことを示している。

すなわち、要素数$ Nの離散フーリエ変換は、要素数$ \frac{N}{2}の離散フーリエ変換を2回計算することによって、計算することができる。

このとき、2つの離散フーリエ変換の計算量はそれぞれ$ O(\frac{N}{2})であるので、元の離散フーリエ変換の計算量は$ O(\frac{N^2}{4} + \frac{N^2}{4}) = O(\frac{N^2}{2})と半分になる。

そして、分解した離散フーリエ変換をさらに再帰的に要素数半分の離散フーリエ変換に分解していくことができ、$ Nは2の冪乗であると仮定しているので、$ \log_2 N回繰り返すことで最終的には要素数1の離散フーリエ変換になる。

したがって、各繰り返しにおける積の計算があるので、最終的な計算量は$ O(N \log_2N)とすることができる。

ここで、実際に高速フーリエ変換の計算を考える。

$ N = 2 のとき

$ F[0] = f[0]W_2^0 + f[1]W_2^0

$ F[1] = f[0]W_2^0 + f[1]W_2^{(0 + \frac{2}{2})(0+1)} = f[0]W_2^0 + f[1]W_2^1

である。

ここで、記述を簡単化するために行列で表現することを考えると、$ F[k]=X_{2,k}, f[n]=x_n として、

$ \begin{pmatrix}X_{2,0}\\X_{2,1}\end{pmatrix}= \begin{pmatrix} W_2^0&W_2^0\\ W_2^0&W_2^1\end{pmatrix} \begin{pmatrix}x_0\\x_1\end{pmatrix} = \begin{pmatrix} 1&1\\ 1&-1\end{pmatrix} \begin{pmatrix}x_0\\x_1\end{pmatrix}

となる。

また、$ N = 4 のとき$ F[k]=X_{4,k},f[n]=x_n とすると（$ W = W_4と略記）、

$ \begin{pmatrix}X_{4,0}\\X_{4,1}\\X_{4,2}\\X_{4,3}\end{pmatrix}= \begin{pmatrix} W^0&W^0&W^0&W^0\\ W^0&W^1&W^2&W^3\\ W^0&W^2&W^4&W^6\\ W^0&W^3&W^6&W^9 \end{pmatrix} \begin{pmatrix}x_0\\x_1\\x_2\\x_3\end{pmatrix}

周波数間引きのために、$ kについて偶奇で分けて変形すると、

$ {\begin{aligned}{\begin{pmatrix}X_{4,0}\\X_{4,2}\\X_{4,1}\\X_{4,3}\end{pmatrix}}&={\begin{pmatrix} W^0&W^0&W^0&W^0\\ W^0&W^2&W^4&W^6\\ W^0&W^1&W^2&W^3\\ W^0&W^3&W^6&W^9 \end{pmatrix}}\begin{pmatrix}x_0\\x_1\\x_2\\x_3\end{pmatrix} \\ &= {\begin{pmatrix} W^0 & W^0 & W^0 & W^0 \\ W^0 & -W^0 & W^0 & -W^0 \\ W^0 & W^1 & -W^0 & -W^1 \\ W^0 & -W^1 & -W^0 & W^1\\ \end{pmatrix}}{\begin{pmatrix}x_{0}\\x_{1}\\x_{2}\\x_{3}\end{pmatrix}} \quad (\because W^{k+N} = W^{k}, W^{k+\frac{N}{2}} = -W^{k}) \\ &= {\begin{pmatrix} 1 & 1 & 1 & 1 \\ 1 & -1 & 1 & -1 \\ 1 & W^1 & -1 & -W^1 \\ 1 & -W^1 & -1 & W^1\\ \end{pmatrix}}{\begin{pmatrix}x_{0}\\x_{1}\\x_{2}\\x_{3}\end{pmatrix}} \end{aligned}}

ここで、$ N = 2における回転因子の行列を

$ W'_2= \begin{pmatrix} W_2^0&W_2^0\\ W_2^0&W_2^1\end{pmatrix} = \begin{pmatrix} 1&1\\ 1&-1\end{pmatrix}

とし、行列$ D_4を

$ D_4 = \begin{pmatrix} 1&W^1\\ 1&-W^1\end{pmatrix}

とすると、

$ \begin{pmatrix}X_{4,0}\\X_{4,2}\\X_{4,1}\\X_{4,3}\end{pmatrix}=\begin{pmatrix}W_2' & W_2' \\ D_4 & -D_4\end{pmatrix}\begin{pmatrix}x_{0}\\x_{1}\\x_{2}\\x_{3}\end{pmatrix}

と表現でき、$ N = 4の計算において$ N = 2で計算した値を利用することで、計算量や計算負荷が削減できることが分かる。

このとき、高速フーリエ変換の計算手順を図式化すると、蝶のような図になることからバタフライ演算ともいわれる。

https://scrapbox.io/files/5fd63b8774f064002247499b.png

$ N = 2のバタフライ図

https://scrapbox.io/files/5fd68b46b22e8c001d61d998.pnghttps://scrapbox.io/files/5fd63bb8af9b05001d43d1a5.png

$ N = 4のバタフライ図

ここで、時間間引きについても考える。

離散フーリエ変換における総和内の各要素を添字$ nの偶奇で分けると

$ \begin{aligned} F[k] &= \sum_{n = 0}^{N - 1}f[n]e^{-j\frac{2\pi}{N}nk} \\ &= \sum_{n = 0}^{\frac{N}{2} - 1}f[2n]e^{-j\frac{2\pi}{N}2nk} + \sum_{n = 0}^{\frac{N}{2}-1}f[2n + 1]e^{-j\frac{2\pi}{N}(2n+1)k} \\ &= \sum_{n = 0}^{\frac{N}{2}-1}f[2n]e^{-j\frac{2\pi}{\frac{N}{2}}nk} + e^{-j\frac{2\pi}{N}k}\sum_{n = 0}^{\frac{N}{2}-1}f[2n+1]e^{-j\frac{2\pi}{\frac{N}{2}}nk} \end{aligned}

ここで、$ W_N = e^{-j\frac{2\pi}{N}}より、

$ F[k] = \sum_{n=0}^{\frac{N}{2}-1}f[2n]W_{\frac{N}{2}}^{nk} + W_N^k\sum_{n=0}^{\frac{N}{2}-1}f[2n+1]W_{\frac{N}{2}}^{nk}

となる。

$ (3)式の右辺について、第一項は要素数$ \frac{N}{2}の離散フーリエ変換であり、第二項は要素数$ \frac{N}{2}の離散フーリエ変換に$ W_N^kをかけたものである。

すなわち、要素数$ Nの離散フーリエ変換は、要素数$ \frac{N}{2}の離散フーリエ変換2つと片方に$ W_N^kをかけて和を取った処理に分解できる。

よって周波数間引きと同様に、2つの離散フーリエ変換の計算量はそれぞれ$ O(\frac{N}{2})であるので、元の離散フーリエ変換の計算量は$ O(\frac{N^2}{4} + \frac{N^2}{4}) = O(\frac{N^2}{2})と半分になる。

したがって、各繰り返しにおける積の計算があるので、最終的な計算量は$ O(N \log_2N)とすることができる。

/icons/hr.icon

ラプラス変換

非周期関数$ f(x)を求めるフーリエ解析は、フーリエ変換を用いて

$ f(x) \sim \frac{1}{2\pi}\int_{-\infty}^{\infty}F(\omega)e^{j\omega x}d\omega

$ F(\omega) = \int_{-\infty}^{\infty}f(t)e^{-j\omega t}dt

で表されるが、フーリエ変換は利用頻度の高い関数などでも、積分が発散してしまう問題がある。

例えば$ f(t) = 1のとき、

$ F(\omega) = \int_{-\infty}^{\infty}f(t)e^{-j\omega t}dt = \int_{-\infty}^{\infty}e^{-j\omega t}dt = 収束しない

となり、計算ができない。

そこで、関数$ f(t)を積分できる形に変形することを考える。

まず、関数$ f(t)は$ f(t) = 0 \quad (t < 0)であるとする。

ここで、$ a > 0のとき$ t \to \inftyで非常に速く$ 0に収束する関数$ e^{-at}を考え、この関数と$ f(t)との積を考える。

つまり、関数$ f(t)に非常に速く収束する関数$ e^{-at}をかけることで、関数$ f(t)に収束性を与えるのである。

https://scrapbox.io/files/5fc91981f20908001f46eedb.jpghttps://scrapbox.io/files/5fc919833cf7510036ce2f23.jpg

関数の積$ f(t)e^{-at}のフーリエ変換を考えると、

$ \int_{-\infty}^{\infty}f(t)e^{-at}e^{-j\omega t}dt = \int_0^{\infty}f(t)e^{-st}dt \cdots (2)

$ s = a + j\omega

が得られ、この$ (2)式の右辺を

$ \mathcal{L}[f(s)] = F(s) = \int_0^{\infty}f(t)e^{-st}dt \cdots (3)

と表し、$ (3)式を関数$ fの（片側）ラプラス変換という。

このとき、$ s = a + j\omegaは複素数であり、その存在領域を複素領域、s領域、s空間などといい、変数$ tの存在する空間をt領域またはt空間という。

なお、関数$ f(t)を全ての実数について考えることで

$ \mathcal{B}[f(s)] = F(s) = \int_{-\infty}^{\infty}f(t)e^{-st}dt

を得ることができ、これを両側ラプラス変換という。

ここで、$ (3) 式のフーリエ逆変換は$ \mathcal{F^{-1}}[F(\omega)] \coloneqq \frac{1}{2\pi}\int_{-\infty}^{\infty} F(\omega)e^{j\omega t}d\omega より、

$ f(t)e^{-at} = \frac{1}{2\pi}\int_{-\infty}^{\infty}F(s)e^{j\omega t}d\omega

両辺を$ e^{-at}で割って

$ f(t) = \frac{1}{2\pi}\int_{-\infty}^{\infty}F(s)e^{at + j\omega t}d\omega

$ s = a + j\omegaなので積分変数を$ \omegaから$ sに変換すると$ \frac{ds}{d\omega} = jであり、積分範囲は$ a - j\inftyから$ a + j\inftyとなるので、

$ f(t) = \frac{1}{2\pi j}\int_{a - j\infty}^{a + j\infty}F(s)e^{st}ds \cdots (4)

となり、この$ (4)式をラプラス逆変換といい、

$ \mathcal{L}^{-1}[f(s)] = \frac{1}{2\pi j}\int_{a - j\infty}^{a + j\infty}F(s)e^{st}ds

と表す。

ここで、ラプラス変換の収束性を考える。

ラプラス変換が$ s_0 = a_0 + j\omega_0で求められたとすると、$ s = a + j\omega(a \geq a_0)となる$ sにおいて

$ |e^{-sx}| = |e^{-(a + j\omega)x|}| = |e^{-ax}e^{-j\omega x}| = e^{-ax} \leq e^{-a_0 x} = |e^{-s_0 x}|

となり、その$ sにおいてもラプラス変換は収束する。

したがって、与えられた関数$ f(t) に対して$ Re[s] < a ならラプラス変換が発散し、$ Re[s] \geq a ならラプラス変換が収束するような実数$ aがただ一つ定まる。

この$ aのことをラプラス変換の収束座標といい、$ a = -\inftyのとき全ての$ sについてラプラス変換は収束し、$ a = \inftyのとき全てのラプラス変換は発散する。

また、収束座標$ aより大きい任意の実数を$ \alpha（$ \forall\alpha \in \mathbb{R}, \alpha > a）とする時、$ f(t)が任意の有限区間で区分的に滑らかであれば、$ (4)式について

$ \frac{1}{2\pi j}\int_{\alpha-j\infty}^{\alpha+j\infty}F(s)e^{sx}ds = \frac{\lim_{a \to t-0}f(a) + \lim_{a \to t+0}f(a)}{2}

が成り立つことが知られており、これをラプラス変換の反転公式という。

ここまで、フーリエ変換を拡張してラプラス変換を導出した。

ある関数$ f(t)がフーリエ変換できない時、$ f_1(t) = f(t)e^{-st}を考えることで$ f(t)に収束性を与え、$ f_1(t)に対してフーリエ変換をおこなう、この一連の流れがラプラス変換である。

フーリエ解析とは、ある強さ（フーリエ係数）を持つ波（複素正弦波$ e^{j\omega t}）を重ね合わせることで任意の関数を表せないか、というものであったが、フーリエ変換は単調な波の和を考えているため$ t \to \inftyで$ 0に収束しない関数を変換できなかった。

そこで、ラプラス変換では増大する波（複素正弦波$ e^{st}）を導入し、指数関数的に拡散する波の和を考えることで、元の関数の発散を抑えフーリエ変換を可能とするのである。

フーリエ変換：https://scrapbox.io/files/5fca9ee6bf7541001c21289b.jpg　ラプラス変換：https://scrapbox.io/files/5fca9ee7770992001cf58fa3.jpg

/icons/hr.icon

Z変換

ラプラス変換によって、非周期で連続な関数$ f(x)に収束性を与えることでフーリエ変換が計算できる。

そこで、関数$ f が不連続な関数$ f[n] である場合についても、収束性を与えることを考える。

非周期で不連続な関数$ f[n] のフーリエ変換は、離散時間フーリエ変換より

$ F(\omega) = \sum_{n = -\infty}^{\infty}f[n]e^{-j\omega n}

で求められるが、ここでは関数$ f[n] は$ f[n] = 0 \quad (n < 0) であるとする。

ここで、ラプラス変換と同様に$ a > 0 のとき$ n \to \infty で非常に速く$ 0 に収束する関数$ e^{-an} を考え、関数の積$ f[n]e^{-an} の離散時間フーリエ変換を考えると、

$ \begin{aligned} \sum_{n = 0}^{\infty}f[n]e^{-an}e^{-j\omega n} &= \sum_{n = 0}^{\infty}f[n]e^{-(a+j\omega)n} \\ &= \sum_{n = 0}^{\infty}f[n]z^{-n} \cdots (5) \end{aligned}

$ z = e^{a+j\omega}

が得られ、この$ (5)式の右辺を

$ \mathcal{Z}[f[n]] = X(z) = \sum_{n=0}^{\infty}f[n]z^{-n} \cdots (6)

と表し、$ (6)式を関数$ fの（片側）z変換という。

このとき、$ z= e^{a + j\omega}は複素数である。

なお、関数$ f[n] を全ての実数について考えることで

$ \mathcal{Z}[f[n]] = X(z) = \sum_{n=-\infty}^{\infty}f[n]z^{-n}

を得ることができ、これを両側z変換という。

ここで、$ z変換の収束性を考える。

$ a を定数として、関数$ f[n] = a^n のとき、z変換は

$ \begin{aligned} X(z) &= \sum_{n=0}^{\infty}a^nz^{-n} \\ &= \sum_{n=0}^{\infty}(az^{-1})^n \\ &= \frac{1}{1-(az^{-1})} \cdots (7) \end{aligned}

ここで、$ (7) 式の無限等比級数が収束するためには、$ |az^{-1}| < 1 \Leftrightarrow |a| < |z|を満たす必要がある。