フーリエ変換とラプラス変換

はじめに

こちらの参考書の第二章、Fourier and Wavelet Transformsの個人的備忘として投稿。後半の章をいきなり読んでもこの章の知識が抜けているとキャッチアップできない部分があったので、今一度復習を兼ねて学習する。
本稿はフーリエ級数から始まり、フーリエ変換、離散フーリエ変換/高速フーリエ変換、ガボール変換について説明し、最後にラプラス変換について記述する。

概要

数理物理学と工学数学の分野では、方程式をより簡単、分離可能、分析可能な座標系へ変換することが主な関心事である。この概念は、J.-B.Joseph Fourierによる1800年代初頭の熱理論研究に端を発し、正弦関数$\sin$と余弦関数$\cos$が解の関数空間における直交基底を提供することを示した。実際にフーリエ変換基底は熱方程式の固有関数であり、周波数は固有値、振幅は境界条件によって決定される。これがヒルベルト空間、演算子理論などの基盤を形成し、高速フーリエ変換（FFT）は現代の計算数学の要となった。近年では、データ駆動型の特異値分解（SVD）やウェーブレットなど、より複雑な基底が利用されるようになっている。本稿では、Fourier Transformsの多岐にわたる応用を探る。

Fourier Series and Fourier Transforms

データのフーリエ変換実装を解説する前にフーリエ級数およびフーリエ変換の導入をする。

関数およびベクトルの内積

まずは関数の内積、エルミート内積（Hermitian inner product）を導入する。変数$x\in[a,b]$について定義された2つの関数$f(x)$と$g(x)$があり、

\begin{align}
\langle f(x),g(x)\rangle=\int^b_af(x)\bar{g}(x)dx\tag{1}
\end{align}

ここで、$\bar{g}(x)$は複素共役を表す。
初見だとこの定義式は意味深に思えるが、ベクトルの内積を考えてみるとクリアになる。関数$f(x)$と$g(x)$をデータベクトルに離散化した場合、サンプリングの解像度が上がるにつれて、ベクトルの内積が関数の内積に収束するようにしたい。以下にそのイメージ図を掲載する。

左図がサンプリグング数=20、右図がサンプリング数=100のグラフだが、サンプリング数$n\rightarrow\infty$のとき、ベクトルの内積が関数の内積と解釈できるようになる。データベクトル$\boldsymbol{f}=\begin{bmatrix}f_1&f_2&\cdots&f_n\end{bmatrix}^T$と$\boldsymbol{g}=\begin{bmatrix}g_1&g_2&\cdots&g_n\end{bmatrix}^T$の内積は

\begin{align}
\langle\boldsymbol{f},\boldsymbol{g}\rangle
=\boldsymbol{g}*\boldsymbol{f}
=\sum_{k=1}^n f_k\bar{g}_k
=\sum_{k=1}^n f(x_k)\bar{g}(x_k)
\end{align}

この値の大きさはデータ点が増えれば増えるほど大きくなる。そこで、$\Delta x=(b-a)/(n-1), b=x_n, a=x_1$を導入することによる正規化を行う。

\begin{align}
\cfrac{b-a}{n-1}\langle\boldsymbol{f},\boldsymbol{g}\rangle
=\sum_{k=1}^n f(x_k)\bar{g}(x_k)\Delta x\tag{2}
\end{align}

これは式$(1)$のRiemann近似である。すると、$n\rightarrow\infty$の時、Riemann積分となり式$(2)$は式$(1)$に収束することが分かる。

この内積は関数のノルムを定義付けることもできる。

\begin{align}
\|f\|_2=(\langle f,f\rangle)^\frac{1}{2}=\sqrt{\langle f,f\rangle}
=\left(\int^b_af(x)\bar{f}(x)dx\right)^\frac{1}{2}
\end{align}

有限なノルムを持つ関数は、$L^2([a,b])$で表される自乗可積分函数の集合として定義され、これはルベーグ可積分関数としても知られる。（余談）

フーリエ級数

関数$f(x)$が周期的かつ部分的に滑らかならば、その関数はフーリエ級数で表現できる。

\begin{align}
f(x)=\frac{a_0}{2}+\sum_{k=1}^\infty(a_k\cos{(kx)}+b_k\sin(kx))\tag{3}
\end{align}

そして、係数$a_k$と$b_k$は以下の式で与えられる。

\begin{align}
\begin{cases}
a_k=\frac{1}{\pi}\displaystyle\int_{-\pi}^{\pi}f(x)\cos(kx)\,dx \\
b_k=\frac{1}{\pi}\displaystyle\int_{-\pi}^{\pi}f(x)\sin(kx)\,dx
\end{cases}
\end{align}

このフーリエ級数は関数を余弦$\cos$と正弦$\sin$の基底$\lbrace\cos{(kx)},\sin{(kx)}\rbrace^\infty_{k=0}$へ射影することによって得られる座標とみなすことができる。
また、この式は

$\|\cos{(kx)}\|_2^2=\|\sin{(kx)}\|_2^2=\pi$

\begin{align}
\|\cos{(kx)}\|_2^2
&=\int_{-\pi}^\pi\cos{(kx)}\overline{\cos}{(kx)}\,dx \\
&=\int_{-\pi}^\pi\cos{(kx)}^2\,dx \\
&=\int_{-\pi}^\pi\frac{1+\cos{(2kx)}}{2}\,dx \\
&=\pi+\frac{1}{2}\int_{-\pi}^\pi\cos{(2kx)}\,dx=\pi
\end{align}

であることを用いれば、

\begin{align}
a_k
&=\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)\cos(kx)\,dx \\
&=\frac{1}{\|\cos{(kx)}\|^2}\int_{-\pi}^{\pi}f(x)\overline{\cos}(kx)\,dx \\
&=\frac{1}{\|\cos{(kx)}\|^2}\langle f(x),\cos{(kx)}\rangle
\end{align}

と式$(1)$の内積を用いて表すことができる。（正弦に関しても同様）

次に、式$(3)$から周期$L$の関数を半開区間$[0,L)$で考える。この時のフーリエ級数は

\begin{align}
f(x)=\frac{a_0}{2}+\sum_{k=1}^\infty\left(a_k\cos{\left(\frac{2\pi kx}{L}\right)}+b_k\sin{\left(\frac{2\pi kx}{L}\right)}\right)
\end{align}

そして、係数$a_k$と$b_k$は以下の式で与えられる。

\begin{align}
\begin{cases}
a_k=\cfrac{1}{L}\displaystyle\int_{0}^{L}f(x)\cos{\left(\frac{2\pi kx}{L}\right)}\,dx \\
b_k=\cfrac{1}{L}\displaystyle\int_{0}^{L}f(x)\sin{\left(\frac{2\pi kx}{L}\right)}\,dx
\end{cases}
\end{align}

ここまで、正弦と余弦を用いて関数を拡張してきたが、Eulerの公式を用いれば、指数関数で拡張することも出来る。複素数$c_k:=\alpha_k+i\beta_k, \alpha_k,\beta_k\in\mathbb{R}$として、

\begin{align}
f(x)
&=\sum_{k=-\infty}^\infty c_ke^{ikx} \tag{4}\\
&=\sum_{k=-\infty}^\infty\left(\alpha_k+i\beta_k\right)\left(\cos{(kx)}+i\sin{(kx)}\right) \\
&=(\alpha_0+i\beta_0) 
+\sum_{k=1}^\infty \left[(\alpha_{-k}+\alpha_k)\cos{(kx)}+(\beta_{-k}-\beta_k)\sin{(kx)}\right]
+i\sum_{k=1}^\infty \left[(\beta_{-k}+\beta_k)\cos{(kx)}-(\alpha_{-k}-\alpha_k)\sin{(kx)}\right]
\end{align}

したがって、整数$k$に対して関数$\psi_k=e^{ikx}$を考えた時に、これは半開区間$[0,2\pi)$で周期的な複素関数基底を与える。この関数が直交することは内積を計算することですぐにわかる。

\begin{align}
\langle\psi_j,\psi_k\rangle
&=\int_{-\pi}^\pi e^{ijx}e^{-ikx}\,dx \\
&=
\begin{cases}
0\quad &\text{if}\,i\neq k,\\
2\pi\quad &\text{if}\,i=k,
\end{cases} \\
&=2\pi\delta_{jk}
\end{align}

ここで、$\delta_{jk}$はKroneckerのデルタ関数である。
つまり、フーリエ級数は関数$f(x)$を正弦と余弦で張られる無限次元の直交関数空間への座標変換を行なっていることになる。式$(4)$は$c_k=(1/(2\pi))\langle f(x),\psi_k(x)\rangle$として、

\begin{align}
f(x)
=\sum_{k=-\infty}^\infty c_ke^{ikx}
=\cfrac{1}{2\pi}\sum_{k=-\infty}^\infty\langle f(x),\psi_k(x)\rangle\psi_k(x)
\end{align}

非連続ハット関数のフーリエ級数表現

pythonで非連続矩形ハット関数のフーリエ級数を計算してみる。元になる関数は半開区間$\left[0,L\right)$で以下のように与える。

\begin{align}
f(x)
=
\begin{cases}
0\quad\text{for}\ x\in\left[0,L/4\right),\\
1\quad\text{for}\ x\in\left[L/4,3L/4\right),\\
0\quad\text{for}\ x\in\left[3L/4,L\right)
\end{cases}
\end{align}

勿論、無限級数を全て足しあげる事は不可能なので、適当な項数で打ち切ることになるのだが、不連続点付近で振幅が乱れる過渡現象が観察できる。これはGibbs現象と一般的に知られる。

import numpy as np
import matplotlib.pyplot as plt

### Define domain
dx = 10**-3
L = 4
x = np.arange(0, L, dx)
n = len(x)
nquart = int(np.floor(n/4))

### Define discontinuous hat function
f = np.zeros_like(x)
f[nquart:3*nquart] = 1

### Compute Fourier series
A0 = np.sum(f * np.ones_like(x)) * dx
fFS = A0 / 2

A = np.zeros(1000)
B = np.zeros(1000)

### Set up the plot
fig, ax = plt.subplots()

### Define points of interest for plotting
plot_points = [1, 10, 100, 1000]

for k in range(1000):
    A[k] = np.sum(f * np.cos(2 * np.pi * (k + 1) * x / L)) * dx
    B[k] = np.sum(f * np.sin(2 * np.pi * (k + 1) * x / L)) * dx
    fFS += A[k] * np.cos(2 * (k + 1) * np.pi * x / L) + B[k] * np.sin(2 * (k + 1) * np.pi * x / L)
    if (k + 1) in plot_points:
        ax.plot(x, fFS, label=f'n = {k + 1}')

### Show the plot
plt.xlabel('x')
plt.ylabel('fFS')
plt.title('Fourier Series Approximation of Discontinuous Hat Function')
plt.legend()
plt.grid(True)
plt.show()

実行結果で以下のような図を得る。

$L=4$としているので、この関数の不連続点は$x=1,3$だが、たしかにこの付近で振幅の過渡現象が見られる。

フーリエ変換

フーリエ級数は周期関数に対して定義されたものであり、その定義式の外で関数は無限に同じ形を繰り返している事になる。
そこで、フーリエ変換が登場する。任意の関数は周期無限の周期関数と捉えることができる。半開区間$x\in\left[-L,L\right)$上でのフーリエ級数について、$L\rightarrow\infty$の極限を考えてみる。

\begin{align}
f(x)
=\frac{a_0}{2}+\sum_{k=1}^\infty\left[a_k\cos{\left(\frac{\pi kx}{L}\right)}+b_k\sin{\left(\frac{\pi kx}{L}\right)}\right]
=\sum_{k=-\infty}^\infty c_ke^{ik\pi x/L}
\end{align}

係数$c_k$は以下の通り。

\begin{align}
c_k
=\cfrac{1}{2L}\langle f(x),\psi_k\rangle
=\cfrac{1}{2L}\int_{-L}^L f(x)e^{-ik\pi x/L}\,dx
\end{align}

ここで、関数$f(x)$は周波数$\omega_k:=k\pi/L$の正弦と余弦で表現されている事に留意する。さて、ここで極限$L\rightarrow\infty$を取ってみる。この時、離散周波数は連続周波数へ近づき、$\Delta\omega:=\pi/L\rightarrow0$として、

\begin{align}
f(x)
&=\lim_{L\to\infty}
\sum_{k=-\infty}^\infty c_ke^{ik\pi x/L} \\
&=\lim_{L\to\infty}
\sum_{k=-\infty}^{\infty}
\cfrac{1}{2L}\int_{-L}^L f(x)e^{-ik\pi x/L}\,dx\,e^{ik\pi x/L} \\
&=\lim_{\Delta \omega \to 0}
\sum_{k=-\infty}^{\infty}
\frac{\Delta \omega}{2\pi}
\underbrace{
\int_{-\pi/\Delta \omega}^{\pi/\Delta \omega}
f(\xi)e^{-ik\Delta \omega \xi}\,d\xi
}_{\langle f(x),\psi_k(x)\rangle}\,e^{ik\Delta\omega x}\quad\because\Delta\omega=\pi/L
\end{align}

極限を取った際、内積$\langle f(x),\psi_k(x)\rangle$は$f(x)$のフーリエ変換となり、$\hat{f}(\omega)\triangleq\mathcal{F}(f(\omega))$と記述される。これを用いれば上記の式もまたリーマン近似となり、以下のように計算できる。

\begin{align}
f(x)
&=\lim_{\Delta \omega \to 0}
\sum_{k=-\infty}^{\infty}
\frac{\Delta \omega}{2\pi}
\langle f(x),\psi_k(x)\rangle
\,e^{ik\Delta\omega x} \\
&=\lim_{\Delta \omega \to 0}
\sum_{k=-\infty}^{\infty}
\frac{\Delta \omega}{2\pi}
\hat{f}(\omega)e^{ik\Delta\omega x} \\
&=\cfrac{1}{2\pi}\int_{-\infty}^\infty\hat{f}(\omega)e^{i\omega x}\,dx \\
&=\mathcal{F}^{-1}(\hat{f}(\omega)) \\\\
\Rightarrow&\,\begin{cases}
\hat{f}(\omega)=\mathcal{F}(f(x))=\displaystyle\int_{-\infty}^\infty f(x)e^{-i\omega x}\,dx \\
f(\omega)=\mathcal{F}^{-1}(\hat{f}(x))=\cfrac{1}{2\pi}\displaystyle\int_{-\infty}^\infty \hat{f}(x)e^{i\omega x}\,dx \\
\end{cases}
\end{align}

このフーリエ変換およびフーリエ逆変換は関数$f(x)$がルベーグ可積分、即ち$L^1(\mathbb{R})$空間に属する時に収束する。

フーリエ変換は線形性や導関数の特性など、扱いやすい性質をもつため、データ分析やコンピュータ科学の分野で慣れ親しまれている。以下の項ではこれら特性について紹介する。

導関数

導関数のフーリエ変換は次のように計算される。

\begin{align}
\mathcal{F}\left(\frac{d}{dx}f(x)\right)
&=\int_{-\infty}^\infty f'(x)e^{-i\omega x}\,dx \\
&=\left[f(x)e^{-i\omega x}\right]_{-\infty}^\infty
-\int_{-\infty}^\infty f(x)\left(-i\omega e^{-i\omega x}\right)\,dx \\
\end{align}

ここで、境界項$\left[f(x)e^{-i\omega x}\right]_{-\infty}^\infty$に関して、$f(x)\in L^1(\mathbb{R}),\lim f(x)=0$を仮定すれば、これは$0$となるので、

\begin{align}
\left[f(x)e^{-i\omega x}\right]_{-\infty}^\infty
-\int_{-\infty}^\infty f(x)\left(-i\omega e^{-i\omega x}\right)\,dx
&=i\omega\int_{-\infty}^\infty f(x)e^{-i\omega x}\,dx \\
&=i\omega\mathcal{F}(f(x)) \\
i.e.\,\mathcal{F}\left(\frac{d}{dx}f(x)\right)
&=i\omega\mathcal{F}(f(x))\tag{5}
\end{align}

よって、$n$階導関数は、

\begin{align}
\mathcal{F}\left(\frac{d^n}{dx^n}f(x)\right)
=i^n\omega^n\mathcal{F}(f(x))
\end{align}

となる。これはフーリエ変換の非常に重要な特性であり、偏微分方程式（Partial Differential EquationPartial Differential Equation, PDEs）を常微分方程式（Ordinary Differential Equation, ODE）へ変換することができるようになる。

\begin{align}
u_{tt} = c u_{xx} \quad \xrightarrow{\mathcal{F}} \quad \hat{u}_{tt} = -c \omega^2  \hat{u}
\end{align}

\begin{align}
\text{(PDE)} \quad \quad \quad \quad \quad \text{(ODE)}
\end{align}

線形性

フーリエ変換は線形作用素であり、以下の条件を満足する。

\begin{align}
\mathcal{F}(\alpha f(x)+\beta g(x))=\alpha\mathcal{F}(f(x))+\beta\mathcal{F}(g(x))
\end{align}

\begin{align}
\mathcal{F}^{-1}(\alpha \hat{f}(\omega)+\beta \hat{g}(\omega))=\alpha\mathcal{F}^{-1}(\hat{f}(\omega))+\beta\mathcal{F}^{-1}(\hat{g}(\omega))
\end{align}

パーセバルの定理（Parseval's theorem）

\begin{align}
\int_{-\infty}^\infty |f(x)|^2\,dx
=\frac{1}{2\pi}\int_{-\infty}^\infty |\hat{f}(\omega)|^2\,d\omega
\end{align}

この式は、$\omega=2\pi t$と置換すれば以下のようにも書き換えられる。（角周波→周波）

\begin{align}
\int_{-\infty}^\infty |f(x)|^2\,dx
=\int_{-\infty}^\infty |\hat{f}(2\pi t)|^2\,dt
\end{align}

即ち、フーリエ変換がエネルギー保存あるいはユニタリ性（≒内積を保存）を持つことが分かる。この性質は近似などで、与えられた切り捨てで誤差を制限する術を提供する。

畳み込み

まず、2つの関数$f(x)$と$g(x)$の畳み込みを定義する。

\begin{align}
(f*g)(x)=\int_{-\infty}^\infty f(x-\xi)g(\xi)\,d\xi
\end{align}

例えば、この式において$\zeta:=x-\xi$などとすると、

\begin{align}
(f*g)(x)
&=\int_{\infty}^{-\infty} f(\zeta)g(x-\zeta)\,d(-\zeta) \\
&=\int_{-\infty}^{\infty} g(x-\zeta)f(\zeta)\,d\zeta
=(g*f)(x)
\end{align}

ここで、$\hat{f}=\mathcal{F}(f),\ \hat{g}=\mathcal{F}(g)$と置くと、

\begin{align}
(f*g)(x)=(g*f)(x)
&=\int_{-\infty}^\infty \underbrace{f(x-\xi)}_{\mathcal{F}^{-1}(\hat{f})}g(\xi)\,d\xi \\
&=\int_{-\infty}^{\infty} \left(\frac{1}{2\pi}\int_{-\infty}^{\infty}\hat{f}(\omega)e^{i\omega(x-\xi)}\,d\omega\right)g(\xi)\,d\xi \\
&=\frac{1}{2\pi}\int_{-\infty}^\infty \int_{-\infty}^\infty\hat{f}(\omega)g(\xi)e^{i\omega(x-\xi)}\,d\omega\,d\xi \\
&=\frac{1}{2\pi}\int_{-\infty}^\infty \hat{f}(\omega)e^{i\omega x}\underbrace{\left(\int_{-\infty}^\infty g(\xi)e^{-i\omega\xi}\,d\xi\right)}_{\hat{g}}\,d\omega \\
&=\frac{1}{2\pi}\int_{-\infty}^\infty \hat{f}(\omega)\hat{g}(\omega)e^{i\omega x}\,d\omega \\
&=\mathcal{F}^{-1}(\hat{f}\hat{g})(x) \\
i.e.\ (f*g)(x)&=(g*f)(x)=\mathcal{F}^{-1}(\hat{f}\hat{g})(x)
\end{align}

この左辺はユークリッド空間における関数の畳み込みを意味し、右辺は周波数領域での関数の乗算を意味しており、これらが同じであることを示している。

簡単な具体例

フーリエ変換は、時間領域の関数$f(x)$を周波数領域の関数$\hat{f}(\omega)$へ変換する手法である。ここで、具体的な例として単純な余弦波$f(x)=\cos{(\omega_0x)}$を考えてみる。

\begin{align}
\mathcal{F}(f(x))
&=\int_{-\infty}^\infty\cos{(\omega_0x)}e^{-i\omega x}\,dx \\
&=\cfrac{1}{2}\int_{-\infty}^\infty\left(e^{i\omega_0x}+e^{-i\omega_0x}\right)e^{-i\omega x}\,dx \\
&=\cfrac{1}{2}\left(
\int_{-\infty}^\infty e^{i\omega_0x}e^{-i\omega x}\,dx+
\int_{-\infty}^\infty e^{-i\omega_0x}e^{-i\omega x}\,dx
\right) \\
&=\cfrac{1}{2}\left(\langle\psi_{\omega_0},\psi_{-\omega}\rangle
+\langle\psi_{-\omega_0},\psi_{-\omega}\rangle\right) \\
&=\pi\left(\delta_{\omega_0,-\omega}+\delta_{\omega_0,\omega}\right)
\end{align}

このように、余弦波のフーリエ変換は2つのクロネッカーのデルタ関数の和となる。これが最大値となる、即ち一番主となる角周波数$\hat{\omega}$を求めてみると、

\begin{align}
\hat{\omega}=\mathop{\rm arg~max}\limits_\omega\left(\pi\left(\delta_{\omega_0,-\omega}+\delta_{\omega_0,\omega}\right)\right)=\pm\omega_0
\end{align}

かくして、$\hat{\omega}=\pm\omega_0$を得るのだが、最初に与えた余弦波の角周波は$\omega_0$であり、偶関数であることから$-\omega_0$も同じく角周波となる。故に、フーリエ変換をすることで元の関数の周波数成分の構成が分かることが実感できた。

離散フーリエ変換と高速フーリエ変換

これまで、連続関数$f(x)$に対するフーリエ級数とフーリエ変換を見てきたが、実際のコンピュータによる計算では離散データからフーリエ変換の近似を求めることになる。離散フーリエ変換（Disrete Fourier transform, DFT）はまさにデータ$\boldsymbol{f}=\begin{bmatrix}f_1 & f_2 & \cdots f_n\end{bmatrix}^T$からフーリエ級数の近似値を求める手法として必要不可欠なものである。しかしながら、計算量が$\mathcal{O}(n^2)$と計算コストが高く、その改善手法として登場したのが高速フーリエ変換（Fast Fourier transform）である。計算量$\mathcal{O}(n\log{n})$

余談になるが、このFFTはIBMのW.Cooleyとプリンストン大学のW.Tukeyによって1965年に開発された手法だが、彼らは一般化された数式を提案しただけであり、構想自体は1805年にGaussが発明していた（紙とペンによる高速フーリエ変換だったらしい笑）。しかもそのFFTの発明もFourierによるフーリエ変換の発表（1807年）2年前だったという、ガウスの奇人ぶりが一層際立つ逸話があるみたいだ。

離散フーリエ変換（DFT）

FFTの導入の前にDFTの簡単な定式化から始める。式は以下の通り。

\begin{align}
\hat{f}_k=\sum_{j=0}^{n-1}f_je^{-i2\pi jk/n}\tag{6}
\end{align}

そして、逆離散フーリエ変換（iDFT）は、

\begin{align}
f_k=\cfrac{1}{n}\sum_{j=0}^{n-1}\hat{f}_je^{i2\pi jk/n}
\end{align}

これは連続フーリエ変換同様に、データ$\boldsymbol{f}$の点を周波数領域$\hat{\boldsymbol{f}}$へ写す線形作用素となる。

\begin{align}
\{f_1,f_2,\cdots,f_n\}\quad\xrightarrow{DFT}\{\hat{f}_1,\hat{f}_2,\cdots,\hat{f}_n\}
\end{align}

$n$個の点が与えられた時、DFTはデータを基本周波数の整数倍の周波数$\omega_n=e^{-2\pi i/n}$を持つ正弦$\sin$と余弦$\cos$を用いてデータを表現する。そしてこれは行列の積演算として次のように計算できる。

\begin{align}
\begin{bmatrix}
\hat{f}_1 \\ \hat{f}_2 \\ \hat{f}_3 \\ \vdots \\ \hat{f}_n
\end{bmatrix}
=
\begin{bmatrix}
1 & 1 & 1 & \cdots & 1 \\
1 & \omega_n & \omega_n^2 & \cdots & \omega_n^{n-1} \\
1 & \omega_n^2 & \omega_n^4 & \cdots & \omega_n^{2(n-1)} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \cdots & \omega_n^{(n-1)^2}
\end{bmatrix}
\begin{bmatrix}
f_1 \\ f_2 \\ f_3 \\ \vdots \\ f_{n}
\end{bmatrix}
\end{align}

\begin{align}
\Rightarrow
\hat{\boldsymbol{f}}=\boldsymbol{Ff}
\end{align}

DFT行列$\boldsymbol{F}$は複素ユニタリVandermonde行列であり、左辺の出力$\hat{\boldsymbol{f}}$は入力$\boldsymbol{f}$のフーリエ係数（=magnitude）と位相を含有し、物理的な解釈にとても有用となる。

先述した連続フーリエ変換について成立するパーセバルの定理（エネルギー保存側）が離散フーリエ変換だと、ユニタリVandermonde行列による線形写像になった、と考えることもできる。

DFT行列の実部の可視化

import numpy as np
import matplotlib.pyplot as plt

n = 256
w = np.exp(-1j * 2 * np.pi / n)

### original DFT
DFT = np.zeros((n, n), dtype=complex)
for i in range(n):
    for k in range(n):
        DFT[i, k] = w**(i * k)
        
### Fast FT
J, K = np.meshgrid(np.arange(n), np.arange(n))
DFT_fast = np.power(w, J * K)
DFT_real = np.real(DFT_fast)

### show real part of the DFT matrix
plt.imshow(DFT_real, cmap='jet', aspect='equal')
plt.colorbar()
plt.title('Real Part of DFT Matrix')
plt.show()

出力結果は以下の通り。

基本周波数$\omega_n=e^{-2\pi i/n}$で生成したDFT行列$\boldsymbol{F}$を上で描画したが、階層的で高い対称性を持つマルチスケール構造が観察できる。各行および列は周波数が増加する余弦波（=実部）となっている。

高速フーリエ変換（FFT）

先述の通り、DFTは行列演算となるので$\mathcal{O}(n^2)$の計算コストが求められる。例えば、一般的なサンプリングレート$44.1{\rm kHz}$のラジオを10秒間取得すると、ベクトル$\boldsymbol{f}$は$44.1\times10^3\times10=4.41\times10^5$次元となる。これのDFT計算ともなれば、二乗になるので、約2000億回もの演算を行うことになる。これに対し、FFTでは$44.1\times10^3\times\lg{(44.1\times10^3)}\simeq6.8\times10^5$と3000倍ほど高速化できる。

参考書の実例と数値が違っています。
FFT $\mathcal{O}(n\log{n})=6\times10^6$で、DFTに対し3万倍高速化、と計算されています。

FFTの基本的なアイデアは計算量$\mathcal{O}(n^2)$のDFTがより効率的に計算可能であることから由来する。例えば、データ数$n=1024=2^{10}$のDTF行列を例に取ると、基本周波数$\omega:=\omega_{1024}$として、

\begin{align}
\hat{\boldsymbol{f}}
&=\boldsymbol{F}_{1024}\boldsymbol{f} \\
&=
\begin{bmatrix}
1 & 1 & 1 & \cdots & 1 \\
1 & \omega & \omega^2 & \cdots & \omega^{1023} \\
1 & \omega^2 & \omega^4 & \cdots & \omega^{2\cdot1023} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & \omega^{1023} & \omega^{2\cdot1023} & \cdots & \omega^{1023^2}
\end{bmatrix}
\begin{bmatrix}
f_1 \\ f_2 \\ f_3 \\ \vdots \\ f_{1024}
\end{bmatrix}
\end{align}

ここで、入力$\boldsymbol{f}$の列を偶数-奇数で分ける。

\begin{align}
\hat{\boldsymbol{f}}
&=
\begin{bmatrix}
1 & 1 & 1 & \cdots & 1 \\
1 & \omega & \omega^2 & \cdots & \omega^{1023} \\
1 & \omega^2 & \omega^4 & \cdots & \omega^{2\cdot1023} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & \omega^{1023} & \omega^{2\cdot1023} & \cdots & \omega^{1023^2}
\end{bmatrix}
\begin{bmatrix}
f_1 \\ f_2 \\ f_3 \\ \vdots \\ f_{1024}
\end{bmatrix} \\
&=
\left[
\begin{array}{cccc:cccc}
1 & 1 & \cdots & 1 & 1 & 1 & \cdots & 1 \\
\omega & \omega^3 & \cdots & \omega^{1023} & 1 & \omega^2 & \cdots & \omega^{1022} \\
\vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\
\omega^{1023} & \omega^{3\cdot1023} & \cdots & \omega^{1023^2} & 1 & \omega^{2\cdot1023} & \cdots & \omega^{1022\cdot1023} 
\end{array}
\right]
\left[
\begin{array}{c}
f_2 \\ f_4 \\ \vdots \\ f_{1024} \\ \hdashline f_1 \\ f_3 \\ \vdots \\ f_{1023}
\end{array}
\right] \\
&=
\left[
\begin{array}{ccc:ccc}
1 & \cdots & 1 & 1 & \cdots & 1 \\
\vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\
\omega^{511} & \cdots & \omega^{1023\cdot511} & 1 & \cdots & \omega^{1022\cdot511} \\
\hdashline
\omega^{512} & \cdots & \omega^{1023\cdot512} & 1 & \cdots & \omega^{1022\cdot512} \\
\vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\
\omega^{1023} & \cdots & \omega^{1023^2} & 1 & \cdots & \omega^{1022\cdot1023}
\end{array}
\right]
\left[
\begin{array}{c}
f_2 \\ f_4 \\ \vdots \\ f_{1024} \\ \hdashline f_1 \\ f_3 \\ \vdots \\ f_{1023}
\end{array}
\right] \\
&=
\left[
\begin{array}{c:c}
\boldsymbol{L} & \boldsymbol{F}' \\
\hdashline
\boldsymbol{M} & \boldsymbol{N}
\end{array}
\right]
\left[
\begin{array}{c}
\boldsymbol{f}_{\rm even} \\\hdashline \boldsymbol{f}_{\rm odd}
\end{array}
\right]
\end{align}

ここで、基本周波数$\omega_{n}=e^{-i2\pi/n}$であることから、$\tilde{\omega}:=\omega_{512}\ i.e.\ \omega^2=\tilde{\omega}$となるので、

\begin{align}
\boldsymbol{F}'
&=
\begin{bmatrix}
1 & 1 & 1 & \cdots & 1 \\
1 & \omega^2 & \omega^4 & \cdots & \omega^{2\cdot511} \\
1 & \omega^4 & \omega^8 & \cdots & \omega^{4\cdot511} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & \omega^{2\cdot511} & \omega^{4\cdot511} & \cdots & \omega^{1022\cdot511}
\end{bmatrix} \\
&=
\begin{bmatrix}
    1 & 1 & 1 & \cdots & 1 \\
    1 & \tilde{\omega} & \tilde{\omega}^2 & \cdots & \tilde{\omega}^{511} \\
    1 & \tilde{\omega}^2 & \tilde{\omega}^4 & \cdots & \tilde{\omega}^{2\cdot511} \\
    \vdots & \vdots & \vdots & \ddots & \vdots \\
    1 & \tilde{\omega}^{511} & \tilde{\omega}^{2\cdot511} & \cdots & \tilde{\omega}^{511^2}
    \end{bmatrix}
=\boldsymbol{F}_{512}
\end{align}

同じ要領で行列$\boldsymbol{L}$,$\boldsymbol{M}$,$\boldsymbol{N}$を$\tilde{\omega}$で表現すると、

\begin{align}
\boldsymbol{L}
&=
\begin{bmatrix}
    1 & 1 & 1 & \cdots & 1 \\
    \omega & \omega\tilde{\omega} & \omega\tilde{\omega}^2 & \cdots & \omega\tilde{\omega}^{511} \\
    \omega^2 & \omega^2\tilde{\omega}^2 & \omega^2\tilde{\omega}^4 & \cdots & \omega^{2}\tilde{\omega}^{2\cdot511} \\
    \vdots & \vdots & \vdots & \ddots & \vdots \\
    \omega^{511} & \omega^{511}\tilde{\omega}^{511} & \omega^{511}\tilde{\omega}^{2\cdot511} & \cdots & \omega^{511}\tilde{\omega}^{511^2}
\end{bmatrix} \\
\boldsymbol{M}
&=
\omega^{512}
\begin{bmatrix}
    1 & 1 & 1 & \cdots & 1 \\
    \omega & \omega\tilde{\omega} & \omega\tilde{\omega}^2 & \cdots & \omega\tilde{\omega}^{511} \\
    \omega^2 & \omega^2\tilde{\omega}^2 & \omega^2\tilde{\omega}^4 & \cdots & \omega^{2}\tilde{\omega}^{2\cdot511} \\
    \vdots & \vdots & \vdots & \ddots & \vdots \\
    \omega^{511} & \omega^{511}\tilde{\omega}^{511} & \omega^{511}\tilde{\omega}^{2\cdot511} & \cdots & \omega^{511}\tilde{\omega}^{511^2}
\end{bmatrix}
=-\boldsymbol{L}
\quad\because\omega^{512}=-1
\\
\boldsymbol{N}
&=
\begin{bmatrix}
    1 & \tilde{\omega}^{512} & \tilde{\omega}^{2\cdot512} & \cdots & \tilde{\omega}^{511\cdot512} \\
    1 & \tilde{\omega}^{512}\tilde{\omega} & \tilde{\omega}^{2\cdot512}\tilde{\omega}^2 & \cdots & \tilde{\omega}^{511\cdot512}\tilde{\omega}^{511} \\
    1 & \tilde{\omega}^{512}\tilde{\omega}^2 & \tilde{\omega}^{2\cdot512}\tilde{\omega}^4 & \cdots & \tilde{\omega}^{511\cdot512}\tilde{\omega}^{2\cdot511} \\
    \vdots & \vdots & \vdots & \ddots & \vdots \\
    1 & \tilde{\omega}^{512}\tilde{\omega}^{511} & \tilde{\omega}^{2\cdot512}\tilde{\omega}^{511} & \cdots & \tilde{\omega}^{511\cdot512}\tilde{\omega}^{511^2}
    \end{bmatrix}
=\boldsymbol{F}_{512}
\quad\because\tilde{\omega}^{512}=1
\end{align}

さらに、複素対角行列$\boldsymbol{D}_{512}\in\mathbb{C}^{512\times 512}$を導入すると、

\begin{align}
\boldsymbol{D}_{512}
&=
\begin{bmatrix}
1 & 0 & 0 & \cdots & 0 \\
0 & \omega & 0 & \cdots & 0 \\
0 & 0 & \omega^2 & \cdots & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
0 & 0 & 0 & \cdots & \omega^{511}
\end{bmatrix}
\quad\Rightarrow\quad
\begin{cases}
\boldsymbol{L}=\boldsymbol{D}_{512}\boldsymbol{F}_{512} \\\\
\boldsymbol{M}=-\boldsymbol{D}_{512}\boldsymbol{F}_{512}
\end{cases} \\
i.e.\ 
\hat{\boldsymbol{f}}
&=
\left[
\begin{array}{c:c}
\boldsymbol{L} & \boldsymbol{F}' \\
\hdashline
\boldsymbol{M} & \boldsymbol{N}
\end{array}
\right]
\left[
\begin{array}{c}
\boldsymbol{f}_{\rm even} \\\hdashline \boldsymbol{f}_{\rm odd}
\end{array}
\right] \\
&=
\begin{bmatrix}
\boldsymbol{D}_{512} & \boldsymbol{I}_{512} \\
-\boldsymbol{D}_{512} & \boldsymbol{I}_{512}
\end{bmatrix}
\begin{bmatrix}
\boldsymbol{F}_{512} & \\
 & \boldsymbol{F}_{512}
\end{bmatrix}
\begin{bmatrix}
\boldsymbol{f}_{\rm even} \\ \boldsymbol{f}_{\rm odd}
\end{bmatrix}
\end{align}

結論は同じになるが、参考書では、$f_{\rm even}$と$f_{\rm odd}$の選び方が誤って入れ替わっており、行列$\boldsymbol{F}_{1024}$の分解が左右逆になっています。

このように、$n=2^p$の時はドミノ倒しの要領で$\boldsymbol{F}_{2^p}\rightarrow\boldsymbol{F}_{2^{p-1}}\rightarrow\boldsymbol{F}_{2^{p-2}}\cdots$と繰り返し表現できる。$n\neq2^p$の時は行列のサイズが２の指数乗になるまでゼロパディングすれば良い。

FFTはデータ$\boldsymbol{f}$の行を偶奇分けすることで、より小さいDFT計算に持ち込み高速化を実現している。

ノイズフィルター

本項ではFFTの簡単な例として信号のノイズキャンセリングを取り扱う。
まず、関数$f(t)$を以下のように定める。

\begin{align}
f(t)=\sin{(2\pi f_1t)}+\sin{(2\pi f_2t)}
\end{align}

ここで、周波数$f_1=50,\ f_2=120$とする。この信号に対して強烈なガウシアンホワイトノイズを加える。

import numpy as np
import matplotlib.pyplot as plt

### Create a simple signal with two frequencies
dt = 10**-3
t = np.arange(0,1,dt)
f_clean = np.sin(2*np.pi*50*t) + np.sin(2*np.pi*120*t)
noise = 2.5*np.random.randn(len(t))

### Plot f and f_clean
plt.figure(figsize=(12, 6))
plt.plot(t, noise, label='Noisy', color='red')
plt.plot(t, f_clean, label='Clean', color='blue')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.legend()
plt.show()

出力結果は以下の通り。

Power Spectol Density, PSDは信号の各周波成分エネルギーの分布を示す指標であり、これはFFTの結果を正規化、つまり$\hat{\boldsymbol{f}}$の二乗で割ることによって得られる。

### Compute the Fast Fourier Transform
n = len(t)
f = f_clean + noise
fhat = np.fft.fft(f, n)
PSD = (fhat * np.conj(fhat) / n).real
freq = (1/(dt*n)) * np.arange(n)
L = np.arange(1, np.floor(n/2), dtype='int')

### Plot
plt.figure(figsize=(12, 6))
plt.plot(freq[:n // 2], PSD[:n // 2])
plt.xlabel('Frequency [Hz]')
plt.ylabel('Power/Frequency [dB/Hz]')
plt.grid()
plt.show()

出力結果は以下の通り。

この結果を見ると、ノイズを含む信号が周波数$50{\rm Hz}$と$120{\rm Hz}$でピークを持つことがわかる。つまり、ある閾値を決定すれば、それ以下の周波数帯を除去することでノイズキャンセルが実現できる。

### Use the PSD to filter out noise
indices = PSD > 100
PSDclean = PSD * indices
fhat = indices * fhat
ffilt = np.fft.ifft(fhat).real

### Plot
plt.figure(figsize=(12, 6))
plt.plot(t, ffilt, label='Filtered', color='orange')
plt.plot(t, f_clean, label='Clean', color='blue')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.legend()
plt.show()

出力結果は以下の通り。

PSDで周波数をフィルタイング後に逆フーリエ変換で元の信号に戻した系列と元の系列を見比べるとかなりマッチしていることが観察できる。

スペクトル微分

次の例としてスペクトル微分による高速かつ正確な微分計算を実演する。イメージとしては下記図のようになる。

式$(5)$で述べた通り、連続関数$f$の導関数のフーリエ変換は、

\begin{align}
\mathcal{F}\left(\cfrac{df}{dx}\right)=i\omega\mathcal{F}\left(f\right)
\end{align}

となる。すなわち、スペクトル微分の基本的なアイデアは対象の関数$f$を一度周波数領域へ持ち込み、$i\omega$を掛けることで微分を再現し、逆フーリエ変換することで導関数$f'$を求めるというものである。

ここでは、離散データベクトルに対するスペクトル微分を考える。DFTの定義は式$(6)$から、

\begin{align}
\hat{f}_k=\sum_{j=0}^{n-1}f_je^{-i2\pi jk/n}
\end{align}

ここで、$\kappa:=2\pi k/n$と置くと

\begin{align}
\hat{f}_k&=\sum_{j=0}^{n-1}f_je^{-ij\kappa} \\
\Rightarrow
\hat{f'}_k&=i\kappa\sum_{j=0}^{n-1}f_je^{-ij\kappa}\tag{7}
\end{align}

このように、離散データベクトルの数値微分は、ベクトル$\hat{\boldsymbol{f}}$のDFTの各成分に$i\kappa$を掛けることで得ることができる。
例として実際に解析的な微分が可能な関数$f(x)$を次のように取ってみる。

\begin{align}
f(x)&=\cos{(x)}e^{-\frac{x^2}{25}} \\
\Rightarrow
\cfrac{df}{dx}&=-\sin{(x)}e^{-\frac{x^2}{25}}-\frac{2}{25}xf(x)
\end{align}

これを1.解析微分、2.有限差分法、3.スペクトル微分の3手法で比較してみる。ちなみに、有限（前進）差分法は以下のように差分商でどう関数を近似する手法である。

\begin{align}
\frac{df}{dx}(x_k)\simeq\frac{f(x_{k+1})-f(x_k)}{\Delta x}
\end{align}

import numpy as np
import matplotlib.pyplot as plt

### Create a function and its derivative
n = 128
L = 30
dx = L / n
x = np.linspace(-L/2, L/2 - dx, n)
f = np.cos(x) * np.exp(-x**2/25)
df = - np.sin(x) * np.exp(-x**2/25) - (2/25) * x * f

### Derivative using finite forward difference method
dfFD = np.zeros_like(f)
dfFD[:-1] = (f[1:] - f[:-1]) / dx
dfFD[-1] = (f[0] - f[-1]) / dx 

### Derivative using FFT
fhat = np.fft.fft(f)
kappa = (2 * np.pi / L) * np.arange(-n//2, n//2)
kappa = np.fft.fftshift(kappa)
dfhat = kappa * fhat * (1j)
dfFFT = (np.fft.ifft(dfhat)).real

### Plot
plt.figure(figsize=(12, 6))
plt.plot(x, df, label='True Derivative', color='black')
plt.plot(x, dfFD, label='Finite Difference', color='red', linestyle='--')
plt.plot(x, dfFFT, label='FFT Derivative', color='blue', linestyle='--')
plt.xlabel('x')
plt.ylabel('df/dx')
plt.legend()
plt.show()

np.fft.fftshiftについて
np.fft.fftでは配列$\lbrace\hat{f}_k\rbrace_{k=0}^{n-1}$が得られるが、fftshiftはこれをゼロ周波数成分を中心に移動させる。すなわち、$n=8$の場合

\begin{align}
\lbrace\hat{f}_0,\hat{f}_1,\hat{f}_2,\hat{f}_3,\hat{f}_4,\hat{f}_5,\hat{f}_6,\hat{f}_{7}\rbrace
\implies
\lbrace\hat{f}_4,\hat{f}_5,\hat{f}_6,\hat{f}_7,\hat{f}_0,\hat{f}_1,\hat{f}_2,\hat{f}_3\rbrace
\end{align}

のように配列の順序を変える。
上のコードで、fhatは左の昇順の配列を格納しているが、$\kappa$を生成したnp.arangeは右のゼロ中心の配列$\lbrace-4\sim3\rbrace$の様になり、逆にこちらの方にnp.fft.fftshiftを適用し、$\lbrace0\sim3,-4\sim-1\rbrace$とすることで整合性を取っている。

出力結果は以下の通り。

目視による評価でも確かにスペクトル微分の方が解析微分の結果に近いことがわかる。次に定量評価として、横軸に$n$、縦軸に点全体のRMSEを取り、$n$が増加するにつれて、スペクトル微分の方が有限差分法よりも優位であることを見てみる。

### Calculate Error
def compute_errors(n):
    L = 30
    dx = L / n
    x = np.linspace(-L/2, L/2 - dx, n)
    f = np.cos(x) * np.exp(-x**2/25)
    df = - np.sin(x) * np.exp(-x**2/25) - (2/25) * x * f

    dfFD = np.zeros_like(f)
    dfFD[:-1] = (f[1:] - f[:-1]) / dx
    dfFD[-1] = (f[0] - f[-1]) / dx

    fhat = np.fft.fft(f)
    kappa = (2 * np.pi / L) * np.arange(-n//2, n//2)
    kappa = np.fft.fftshift(kappa)
    dfhat = kappa * fhat * (1j)
    dfFFT = (np.fft.ifft(dfhat)).real

    error_fd = np.abs(dfFD - df)
    error_fft = np.abs(dfFFT - df)
    rmse_fd = np.sqrt(np.mean(error_fd**2))
    rmse_fft = np.sqrt(np.mean(error_fft**2))

    return rmse_fd, rmse_fft


### The range of n
n_values = np.logspace(1, 5, num=20, dtype=int)
print(n_values)
rmse_fd_values = []
rmse_fft_values = []

### Calculate error for each n
for n in n_values:
    rmse_fd, rmse_fft = compute_errors(n)
    rmse_fd_values.append(rmse_fd)
    rmse_fft_values.append(rmse_fft)

### Plot
plt.figure(figsize=(12, 6))
plt.plot(n_values, rmse_fd_values, label='Finite Difference', color='red', marker='o')
plt.plot(n_values, rmse_fft_values, label='Spectral Derivative', color='blue', marker='o')
plt.xscale('log')
plt.xlabel('n')
plt.yscale('log')
plt.ylabel('Error')
plt.grid(True, which="both", ls="--", lw=0.3)  
plt.legend()
plt.show()

出力結果は以下の通り。

両手法ともに$n$が増加、即ち$\Delta x$が小さくなるに従い、誤差は小さくなっているが、スペクトル微分の方がより顕著に収束していくのがわかる。

偏微分方程式の変形

先に述べたように、フーリエ変換は偏微分方程式を常微分方程式へと変換する。本項では、その有用性の例として熱方程式、一方向波動方程式、Burger's方程式を挙げる。

熱方程式

一次元の熱方程式は次のように表せられる。

\begin{align}
\cfrac{\partial u}{\partial t}&=\alpha^2\cfrac{\partial^2 u}{\partial x^2} \\
\Rightarrow
u_t&=\alpha^2u_{xx} \\
\end{align}

$u(t,x)$は時間$t$、位置$x$における温度分布、$\alpha$は熱伝導率である。
ここで、両辺$x$についてフーリエ変換すると

\begin{align}
\hat{u}_t=-\alpha^2\omega^2\hat{u}
\end{align}

これは、周波数$\omega$を固定したと見れば、$t$についての常微分方程式となっていることがわかる。これの解は以下のようになる。

導出はこちら。

\begin{align}
&&
\hat{u}_t&=-\alpha^2\omega^2\hat{u} \\
&\Rightarrow&
\cfrac{\hat{u}_t}{\hat{u}}&=-\alpha^2\omega^2 \\
&\Rightarrow&
\log{(\hat{u})}&=-\alpha^2\omega^2t+C \\
&\Rightarrow&
\hat{u}&=Ce^{-\alpha^2\omega^2t},\quad
C=\hat{u}_{t=0}
\end{align}

\begin{align}
\hat{u}=\hat{u}_0e^{-\alpha^2\omega^2t}
\end{align}

ここで、初期条件$\hat{u}_0=\hat{u}(0,\omega)$は元の方程式の初期条件$u(0,x)$をフーリエ変換すれば得られる。高周波数領域、つまり、$\omega$が大きい値を取る場合には、指数的に減衰率が高くなることが明らかである。

元の関数$u(t,x)$は$\hat{u}(t,\omega)$を逆フーリエ変換すれば良いので、

\begin{align}
&&
u(t,x)&=\mathcal{F}^{-1}\left(\hat{u}(t,\omega)\right) \\
&\Rightarrow&
&=\mathcal{F}^{-1}\left(\hat{u}_0e^{-\alpha^2\omega^2t}\right) \\
&\Rightarrow&
&=u_0*\mathcal{F}^{-1}\left(e^{-\alpha^2\omega^2t}\right) 
\quad\because\mathcal{F}^{-1}(\hat{f}\hat{g})=f*g\\
&\Rightarrow&
&=u_0*\cfrac{1}{\sqrt{4\pi\alpha^2t}}e^{-x^2/(4\alpha^2t)} \\
\end{align}

$\mathcal{F}^{-1}\left(e^{-\alpha^2\omega^2t}\right)$の計算はこちら。

\begin{align}
\mathcal{F}^{-1}\left(e^{-\alpha^2\omega^2t}\right)
&=\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{-\alpha^2\omega^2t} e^{i\omega x}\,d\omega \\
&=\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{-\alpha^2\omega^2t+i\omega x}\,d\omega \\
&=\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{-\alpha^2t(\omega-ix/(2\alpha t))^2-x^2/(4\alpha^2t)}\,d\omega \\
&=\frac{1}{2\pi}\cdot e^{-x^2/(4\alpha^2t)}\cdot\frac{1}{\sqrt{\alpha^2t}}\cdot\sqrt{\pi} \\
&=\cfrac{1}{\sqrt{4\pi\alpha^2t}}e^{-x^2/(4\alpha^2t)}
\end{align}

この結果は、熱方程式の基本解（グリーン関数）に対応する。
さて、次はこの偏微分方程式をFFTを用いて高速かつ正確に数値計算として求めてみる。離散関数についての偏微分方程式は、式$(7)$と同じ要領で、

\begin{align}
\hat{u}_t=-\alpha^2\kappa^2\hat{u}
\end{align}

となることに留意して、周波数領域でscipy.integrate.odeintを用いて常微分方程式を数値的に解き、結果を逆フーリエ変換して位置領域へ戻す。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

### Parameter Setting
a = 1
L = 100
N = 1000
dx = L / N
x = np.arange(-L/2, L/2, dx)

### Definition of discrete wavenumber
kappa = 2 * np.pi * np.fft.fftfreq(N, d=dx)

### Initial conditions
u0 = np.zeros_like(x)
u0[int((L/2 - L/10)/dx):int((L/2 + L/10)/dx)] = 1
u0hat = np.fft.fft(u0)

### The right-hand side of the differential equation
def rhsHeat(uhat_ri, t, kappa, a):
    N = len(kappa)
    uhat = uhat_ri[:N] + 1j * uhat_ri[N:]
    d_uhat = -a**2 * kappa**2 * uhat
    return np.concatenate((d_uhat.real, d_uhat.imag))

### Simulation in Fourier space
dt = 0.1
t = np.arange(0, 10, dt)
u0hat_ri = np.concatenate((u0hat.real, u0hat.imag))
uhat_ri = odeint(rhsHeat, u0hat_ri, t, args=(kappa, a))

### Convert results to real space
uhat = uhat_ri[:, :N] + 1j * uhat_ri[:, N:]
u = np.fft.ifft(uhat, axis=1).real

### Plot
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, projection='3d')
T, X = np.meshgrid(t, x)

### Apply color map
norm = plt.Normalize(u.min(), u.max())

for i in range(0, len(t), 10):
    for j in range(len(x) - 1):
        ax.plot(
            [X[j, i], X[j+1, i]],
            [T[j, i], T[j+1, i]],
            [u[i, j], u[i, j+1]],
            color=plt.cm.jet(norm(u[i, j]))
        )

ax.set_xlabel('x')
ax.set_ylabel('t')
ax.set_zlabel('u(t, x)', labelpad=-1)

plt.show()

出力結果は以下の通り。

図2は、温度分布$u(t, x)$が時間的に変化していく様子を、示したものである。初期条件の尖った角（=不連続点）が最も高い周波数帯に対応するため、これらの点から急速に拡散することがわかる。最終的に温度が一定の分布になるまで、低周波数の変動も減衰する。

一方向波動方程式

2つ目の例として、以下の一方向波動方程式を例にとる。

\begin{align}
&&\frac{\partial u}{\partial t}+c\cfrac{\partial u}{\partial x}= 0 \\
&\Rightarrow&
u_t+cu_x=0
\end{align}

初期条件$u(0,x)$は右方向へ速度$c$で伝搬するものとすれば、すぐに$u(t,x)=x-ct$の解を得る。
先ほどの熱方程式のコードにて、周波数領域での微分方程式を生成する関数と初期条件を置き換えるだけで、再現できる。

### Initial conditions: Gaussian pulse
sigma = L / 10  
u0 = np.exp(-x**2 / (2 * sigma**2))
u0hat = np.fft.fft(u0)

### The right-hand side of the differential equation
def rhsWave(uhat_ri, t, kappa, c):
    N = len(kappa)
    uhat = uhat_ri[:N] + 1j * uhat_ri[N:]
    d_uhat = -c * 1j * kappa * uhat
    return np.concatenate((d_uhat.real, d_uhat.imag))

出力結果は以下の通り。

Burger's方程式

最後に、流体中の衝撃波を生じさせる非線形対流・拡散の1次元の例である、非線形Burger's方程式を考える。

\begin{align}
&&
\cfrac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}=\nu\frac{\partial^2 u}{\partial x^2} \\
&\Rightarrow&
u_t+uu_x=\nu u_{xx}
\end{align}

非線形項$uu_{x}$は移流項（対流項）と呼ばれ、波の急激な変化やショック波の形成に寄与する。$\nu$は粘性係数であり、$\nu u_{xx}$は粘性項（拡散項）と呼ばれる。

Burger's方程式は、その非線形性、つまり$u,\ u_x,\ u_{xx}$を用いるため、各時間ごとに空間領域↔︎周波数領域の変換が必要となる。

### Parameter Setting
nu = 0.001
L = 20
N = 1000
dx = L / N
x = np.arange(-L/2, L/2, dx)

### Initial conditions
u0 = 1/np.cosh(x)

### The right-hand side of the differential equation
def rhsBurgers(u, t, kappa, nu):
    uhat = np.fft.fft(u)
    d_uhat = 1j * kappa * uhat
    dd_uhat = -np.power(kappa, 2) * uhat
    d_u = np.fft.ifft(d_uhat)
    dd_u = np.fft.ifft(dd_uhat)
    du_dt = -u * d_u + nu * dd_u
    return du_dt.real

### Simulation in Fourier space
dt = 0.025
t = np.arange(0,100*dt,dt)
u = odeint(rhsBurgers,u0,t,args=(kappa,nu))

出力結果は以下の通り。

Gabor変換とスペクトログラム

これまで扱ってきたフーリエ変換は信号の周波数情報のみを解析し、時間情報については一切得られない。フーリエ変換は周期的かつ定常な信号にしか特徴付けることが出来ない。非定常な信号に関しては時間情報と周波数情報を同時に特徴付ける必要がある。

Gabor変換（短時間フーリエ変換;STFT）は窓関数を用いて時間的局所ごとに周波数解析を行うことで、時間と周波数の両方の情報を保持しながら信号を解析する手法である。

\begin{align}
\mathcal{G}(f)(t,\omega)
=\hat{f}_\mathcal{G}(t,\omega)
=\int_{-\infty}^\infty f(\tau)e^{-i\omega\tau}\bar{g}(\tau-t)\,d\tau
\end{align}

ここで、上式は次のように関数$g_{t,\omega}$を定義することで、エルミート内積で表現することができる。

\begin{align}
g:=e^{i\omega\tau}g(\tau-t) \\
\Rightarrow
\mathcal{G}(f)(t,\omega)=\langle f,g_{t,\omega}\rangle
\end{align}

また、$g(t)$は窓関数であり、ガウス窓などがよく用いられる。

\begin{align}
g(t)=e^{-(t-\tau)^2/a^2}
\end{align}

パラメータ$a$が大きい程、$g\rightarrow 1$なので時間幅が狭くなり、周波数分解能が上がる。
逆Gabor変換の式は以下の通り。

\begin{align}
f(t)
&=\mathcal{G}^{-1}(\hat{f}_g(t,\omega)) \\
&=\cfrac{1}{2\pi\|g\|^2}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}
\hat{f}_g(\tau,\omega)g(t-\tau)e^{i\omega t}\,d\omega\,d\tau
\end{align}

離散Gabor変換

Gabor変換はフーリエ変換同様に離散データに対しても適用できる。Gabor変換の場合は、時間と周波数両方を離散化する必要がある。

\begin{align}
\nu&=j\Delta\omega \\
\tau&=k\Delta t
\end{align}

この離散化を窓関数に適用すると、

\begin{align}
g_{j,k}
=e^{i2\pi j\Delta\omega t}g(t-k\Delta t)
\end{align}

となり、離散Gabor関数は以下の通り。

\begin{align}
\hat{f}_{j,k}=\int_{-\infty}^{\infty}f(\tau)\bar{g}_{j,k}(\tau)\,d\tau
\end{align}

この積分は離散関数$f$と$\bar{g}_{j,k}$において有限リーマン和を用いることで近似することができる。

2次チャープ

簡単な例として、周波数が二乗で増幅する余弦関数を考えてみる。

\begin{align}
f(t)=\cos{(2\pi t\omega(t))}
\quad\text{where}\quad
\omega(t)=\omega_0+(\omega_1-\omega_0)t^2/t_1^2
\end{align}

この定義から明らかだが、周波数$\omega$は一定ではなく、時間$t:0\rightarrow t_1$で$\omega_0\rightarrow\omega_1$と変化する。
この関数に対して、PSDを適用してみると以下のようになる。

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [12, 8]
plt.rcParams.update({'font.size': 18})

### Create a cosine function whose frequency varies.
dt = 0.001
t = np.arange(0,2,dt)
f0 = 50
f1 = 250
t1 = 2
x = np.cos(2*np.pi*t*(f0 + (f1-f0)*np.power(t,2)/(t1**2)))

### Compute the Fast Fourier Transform
n = len(t)
xhat = np.fft.fft(x, n)
PSD = (xhat * np.conj(xhat) / n).real
freq = (1/(dt*n)) * np.arange(n)
L = np.arange(1, np.floor(n/2), dtype='int')

### Plot
plt.figure(figsize=(12, 6))
plt.plot(freq[:n // 2], PSD[:n // 2])
plt.xlabel('Frequency [Hz]')
plt.ylabel('Power/Frequency [dB/Hz]')
plt.show()

出力結果は以下の通り。

初期周波数$\omega(0)=50$とコード内では定義したので、PSDの結果も$50{\rm Hz}$でピークが立っていることが分かるが、時間に関する情報が全く得られない。
これに対して、スペクトログラムによる結果を可視化してみる。

plt.specgram(x, NFFT=128, Fs=1/dt, noverlap=120,cmap='jet')
cbar = plt.colorbar()
cbar.set_label('Power/Frequency [dB/Hz]')
plt.xlabel('Time [s]')
plt.ylabel('Frequency [Hz]')
plt.show()

出力結果は以下の通り。

PSDではy軸になっていた${\rm Power/Frequency [dB/Hz]}$がz軸方向に与えられカラーマップで表現されている。この結果から周波数$\omega$の時間変化がよく表現されていることが分かった。

ベートベン ピアノソナタ第8番 「悲愴」 第一楽章

ベードーベンの三大ピアノソナタと言えば、「第8番 悲愴」、「第14番 月光」、「第23番 熱情」である。ピアノソナタは複数の楽章で構成された楽曲構造を指し、それぞれ有名な楽章は、第二楽章、第三楽章、第三楽章かなと個人的には思う。参考書では悲愴の第一が取り上げられているが、これは和音が多く使用されており、スペクトログラムのしがいがあるからだろう。このような複雑な音楽に対してもスペクトログラムは有効に働き、ユーザーが録音した音楽の一部を元に曲のタイトルやアーティスト情報を特定するShazamアルゴリズムなどにも利用されている。

データはこちらからVBR MP3をダウンロードして、実際にスペクトログラムを適用してみる。

import IPython.display
import librosa

import matplotlib.pyplot as plt

### Original mp3 file
FS = 24000
y, sr = librosa.load('../data/Beethoven__Piano_Sonata_Pathetique__Arthur_Rubenstein.mp3', sr=FS)

### Extract the first 30 seconds
y = y[:int(sr * 30)]

### Plot
plt.figure(figsize=(20,5))
librosa.display.waveshow(y, sr=sr)

### You can hear the original music
IPython.display.Audio(y, rate=FS)

### Plot the result of spectrogram
plt.specgram(y, NFFT=8192, Fs=FS, noverlap=2048, cmap='RdGy')
cbar = plt.colorbar()
cbar.set_label('Power/Frequency [dB/Hz]')
plt.xlabel('Time [s]')
plt.ylabel('Frequency [Hz]')
plt.show()

出力結果は以下の通り。

不確定性原理

先述したが、周波数解析では時間領域と周波数領域の分解能を両立することに制約があるという、根本的な不確定性原理が存在する。つまり、時間窓を短くすると時間分解能は高くなるが、周波数分解能が低くなり、逆に時間窓を長くすると周波数分解能は高くなるが、時間分解能が低くなる。
数学的には以下のように記述される。

\begin{align}
\left(
    \int_{-\infty}^{\infty}x^2|f(x)|^2\,dx
\right)
\left(
    \int_{-\infty}^{\infty}\omega^2|\hat{f}(\omega)|^2\,d\omega
\right)
\geq
\cfrac{1}{16\pi^2}
\end{align}

上式は関数$f(x)$が絶対連続かつ、$xf(x)$、$f'(x)$共に平方可積分、つまり、

\begin{align}
\int_{-\infty}^\infty |xf(x)|^2\,dx < \infty \\
\int_{-\infty}^\infty |f'(x)|^2\,dx < \infty
\end{align}

の時に成立する。
よく見てみると、関数$x^2|f(x)|^2$は点$x=0$における2次モーメントを表しており、特に$f(x)$がガウス関数の場合は分散となる。つまり、関数$f(x)$とそのフーリエ変換$\hat{f}(\omega)$の分散の積がある有限定数以上であるので、共に局在化させることに縛り（$=\frac{1}{16\pi^2}$）があるということになる。

これらの不確定原理はGabor極限として知られ、量子力学におけるハイゼンベルグの不確定原理へと関連する。

ラプラス変換

ラプラス変換はフーリエ変換と非常に密接な関係であり、微分方程式や制御の理論で広く用いられている。ここでは、ラプラス変換をフーリエ変換の一般化として導出し、いくつかの使用例を提示する。

フーリエ変換はルベーグ可積分、つまり無限遠でゼロに収束するような関数についてwell-definedであった。しかし指数関数、ヘヴィサイド関数、通常の余弦関数などについては、適用不可、あるいは不安定となってしまう。つまり、これらの関数を含むような偏/常微分方程式はフーリエ変換を使うことが出来ないことを意味する。ラプラス変換はフーリエ変換と似た要領で、非ルベーグ可積分関数を含む微分方程式についても有効となる。

ラプラス変換は上述のような急激に増加する関数、無限大で発散する関数に指数関数で重みをつけた片側フーリエ変換と解釈することができる。
ある関数$f(t)$に重み$e^{-\gamma t}$を付けた片側フーリエ変換は

\begin{align}
F(t):&=f(t)e^{-\gamma t}H(t)
\quad\text{where }H(t)\text{ is Heaviside func.} \\
\Rightarrow
\hat{F}(\omega)&=\mathcal{F}(F(t)) \\
&=\int_{-\infty}^\infty F(t)e^{-i\omega t}\,dt \\
&=\int_0^\infty f(t)e^{-\gamma t}e^{-i\omega t}\,dt \\
&=\int_0^\infty f(t)e^{-(\gamma+i\omega)t}\,dt \\
&=\int_0^\infty f(t)e^{-st}\,dt
=\bar{f}(s)
\end{align}

ここで、ラプラス変数$s$は$s:=\gamma+i\omega$とおいたので、複素変数であることに留意する。
逆ラプラス変換$\bar{f}\rightarrow f$についても逆フーリエ変換から導出することができる。

\begin{align}
F(t)
&=\mathcal{F}^{-1}(\hat{F}(\omega)) \\
&=\cfrac{1}{2\pi}\int_{-\infty}^\infty \hat{F}(\omega)e^{i\omega t}\,d\omega \\
\Rightarrow
e^{\gamma t}F(t)
&=\cfrac{1}{2\pi}\int_{\infty}^\infty e^{\gamma t}\hat{F}(\omega)e^{i\omega t}\,d\omega \\
&=\cfrac{1}{2\pi}\int_{\infty}^\infty \hat{F}(\omega)e^{(\gamma+i\omega)t}\,d\omega \\
\end{align}

ここで、ラプラス変数$s=\gamma+i\omega$であることから、$1/i,ds=d\omega$なので

\begin{align}
e^{\gamma t}F(t)
&=\cfrac{1}{2\pi i}\int_{\gamma-i\infty}^{\gamma+i\infty} \hat{F}(\omega)e^{st}\,ds \\
&=\cfrac{1}{2\pi i}\int_{\gamma-i\infty}^{\gamma+i\infty} \bar{f}(s)e^{st}\,ds
\end{align}

故に、まとめるとラプラス変換/逆変換は以下のようになる。

\begin{align}
\bar{f}(s)
&=\mathcal{L}(f(t))
=\int_0^{\infty}f(t)e^{-st}\,dt \\
f(s)
&=\mathcal{L}^{-1}(\bar{f}(s))
=\cfrac{1}{2\pi i}\int_{\gamma-i\infty}^{\gamma+i\infty}\bar{f}(s)e^{st}\,ds
\end{align}

ラプラス変換そのものは関数$f(t)$を時間領域から複素数領域に変換する手法であり、これが有効なのは関数$f(t)$が半無限区間$t>0$で定義されている場合である。そして、逆ラプラス変換により元の時間領域に戻す際にも、元の関数$f(t)$が$t>0$で定義されていることが前提となる。

関数の微分

導関数のラプラス変換は次のように計算できる。

\begin{align}
\mathcal{L}\left(\frac{d}{dt}f(t)\right)
&=\int_0^\infty f'(t)e^{-st}\,dt \\
&=-f(0)-\int_0^\infty f(t)\left(-se^{-st}\right)\,dt \\
&=-f(0)+s\bar{f}(s)
\end{align}

故に、$n$階導関数は

\begin{align}
\mathcal{L}\left(\frac{d^n}{dt^n}f(t)\right)
&=-f^{(n-1)}(0)-sf^{(n-2)}(0)-\cdots-s^{n-1}f(0)+s^n\bar{f}(s) \\
&=-\sum_{k=0}^{n-1}s^kf^{(n-k-1)}(0)+s^n\bar{f}(s)
\end{align}

この性質を用いて、偏微分方程式→微分方程式→代数式へと変換することができる。例として線形二階減衰調和振動子（=バネに取り付けられた質量がダンパーによって減衰するシステム）を考えてみる。

\begin{align}
&&
\ddot{x}+a\dot{x}+bx=0 \\
&\Rightarrow&
(-v_0-sx_0+s^2\bar{x}(s))+a(-x_0+s\bar{x}(s))+b\bar{x}(s)=0 \\
&\Rightarrow&
\underbrace{(s^2+as+b)\bar{x}(s)}_{特性多項式}=\underbrace{sx_0+v_0+ax_0}_{初期条件} \\
&\Rightarrow&
\bar{x}(s)=\frac{sx_0+v_0+ax_0}{s^2+as+b}
\end{align}

このように、$\bar{x}(s)$を代数的に求めれば、逆ラプラス変換で元の解$x(t)$を得ることができる。
実際、$a=5,\ b=4,\ x_0=2,\ v_0=-5$などとして解いてみると、

\begin{align}
\bar{x}(s)
&=\frac{2s+5}{s^2+5s+4}
&=\frac{1}{s+1}+\frac{1}{s+4} \\
\Rightarrow
x(t)&=e^{-t}+e^{-4t}
\end{align}

ここで、$\mathcal{L}(e^{\lambda t})=1/(s-\lambda)$を用いた。

おわりに

次回はみんな大好きウェーブレット変換について記述する。