離散フーリエ変換を心の底からわかりたい

nrnrk

前提

かなり個人メモの位置付けが強いので、省略したりすることも多くなります。自信がついたり、ニーズを感じれば整理して記事にします。

出発点

離散フーリエ変換を通常の(連続)フーリエ変換の離散化バージョンというなんとなくの理解はあった。ただし、関数空間などの数学的な位置付けはわかっていないし、それどころか久々にやると計算もままならない。

個人的なゴール

離散フーリエ変換を三角関数を習ったことのある人に熱量を持って30分くらい語れるようになる

構想

まずは具体的な計算や大枠を捉える
その次に可視化などを踏まえて直観的な理解を目指す
ある程度把握してきたら、数学的な部分もできるだけ整理したい

Hidden comment

nrnrk

用語

(連続)フーリエ変換をFT(Fourier Transform)とし、離散フーリエ変換をDFT(Discrete Fourier Transform)と呼ぶことにする。

MATHEMATICS OF THE DISCRETE FOURIER TRANSFORM (DFT) WITH AUDIO APPLICATIONS

MATHEMATICS OF THE DISCRETE FOURIER TRANSFORM (DFT) WITH AUDIO APPLICATIONS から大枠を把握したい

定義

厳密ではないが、FTを出発点としてDFTの定義を与える。

FT は以下のように定義されている。 $x(t)$ は連続時間で定義された信号(関数)で $X(\omega)$ はそのスペクトル。周波数を角度に換算した、 $\omega$ は角周波数。

X(\omega) = \int_{-\infty}^{\infty} x(t) e^{-i \omega t} dt

これを離散化した以下がDFTと定義する。

X(\omega_k) \equiv \Sigma_{n=0}^{N-1} x(t_n) e^{-i \omega_k t_n} \quad (k = 0, 1, 2, \cdots, N - 1)

大事なポイントは、積分 -> 総和という変換が起きていること。これはリーマン積分の考えからすると違和感はない。

$t$ と $\omega$ の両方 が離散化されているのは注意しておきたい。(多分当たり前なのだけど)
変数の整理を以下にしておく。

$T$ をサンプリングの時間、 $N$ をサンプル数とする。

変数	説明
$T$	サンプリングの間隔の時間 ( $\frac{1}{T}$ が周波数)
$N$	サンプル数
$\Omega$	= $\frac{2 \pi}{N T}$ 角周波数の最小単位^[1]
$t_n$	$n$ 番目のサンプル時間 = $nT$
$\omega_k$	$k$ 番目のサンプル周波数 = $k \Omega$

$n$ と $k$ で使い分けていることからもわかるように、サンプル時間とサンプル周波数は個々に対応関係にはない。単に、離散化したときの時間の各点と離散化したときの周波数の各点を表す座標くらいの認識で考えると誤解がなさそう。

脚注

書籍では radian-frequency sampling interval とあり、直訳するとサンプリング間隔の角周波数ということだが、意味がわかりやすいように表現を変えている。離散フーリエ変換では、サンプリング間隔に対応するこの角周波数が最小の角周波数単位で、 $\omega_N = \frac{2 \pi}{N T}$ が最大の角周波数として扱うと考えることができそう。最小の方は自然だし、最大の方も"サンプリングしている範囲以上の周期性なんてわかりっこない"と思うと自然に感じられる。 ↩︎

nrnrk

今振り返ってみるとそれほど難しく考えなくても良さそうに見える。結局は 関数解析(黒田) にお世話になった。離散バージョンの方は明確な記述は見つけられなかったけど、以下であっているはず。
連続フーリエ変換: L^2(\mathbf{R}^n) \rightarrow L^2(\mathbf{R}^n)

離散時間フーリエ変換(DTFT, 周波数は連続): l^2(\mathbf{Z}) \rightarrow L^2([0, 2\pi]) (無限列)

離散フーリエ変換(DFT): \mathbf{C}^n \rightarrow \mathbf{C}^n (有限列)
と理解すると割とすっきりする。

注意点・ポイントとしては、
連続フーリエ変換は L^1(\mathbf{R}^n) で定義できるが、逆変換が収束するようにL^2(\mathbf{R}^n) で定義している。L^2(\mathbf{R}^n) であって L^1(\mathbf{R}^n) でないケースでは、フーリエ変換が収束しない。その場合は、定数 L (> 0) を導入して、積分の範囲を |x| \leq L して確実に収束させたのち、L \rightarrow \infty とすることで定義できる。ちなみに、L^1(\mathbf{R}^n) の場合だと、L^1(\mathbf{R}^n) \rightarrow C^0(\mathbf{R}) になる。
DFT の方は離散化しているので、有限列でかつ各値も有限になると考えることができる。

nrnrk

残った問題は $L^2(\mathbf{R}^n)$ と $\mathbf{C}^n$ の対応関係とそれがフーリエ変換によってどういう変換がなされるかという点か。この点がわかれば数理の観点はすっきりする。あとは、高速フーリエ変換のアルゴリズムはちゃんと理解したい。

nrnrk

以下のようになればハッピーだけど、そんなわけもなく...

L²(ℝⁿ) ---FT---> L²(ℝⁿ)
  |                 |
sampling         sampling
  ↓                 ↓
 ℂᴺ ----DFT----> ℂᴺ

単純な対応づけは難しそう。実際には帯域制限や各種適切な前処理で近づけていくというのが実践的な話っぽい。

nrnrk

ここは、完全にはすっきりしないが、どういうモヤモヤが残るかが明確になったので、次に高速フーリエ変換のアルゴリズムを理解しよう

nrnrk

まず、離散フーリエ変換は定義として、
W_N = e^{-j\frac{2\pi}{N}}を用いて、次のようにかける。
X_k = \sum_{n=0}^{N-1} x_n \cdot W_N^{kn}行列形式で書くと、
\begin{bmatrix} 
X_0 \\ 
X_1 \\ 
X_2 \\ 
\vdots \\ 
X_{N-1} 
\end{bmatrix} = 
\begin{bmatrix} 
1 & 1 & 1 & \cdots & 1 \\ 
1 & W_N^1 & W_N^2 & \cdots & W_N^{N-1} \\ 
1 & W_N^2 & W_N^4 & \cdots & W_N^{2(N-1)} \\ 
\vdots & \vdots & \vdots & \ddots & \vdots \\ 
1 & W_N^{N-1} & W_N^{2(N-1)} & \cdots & W_N^{(N-1)(N-1)} 
\end{bmatrix} 
\begin{bmatrix} 
x_0 \\ 
x_1 \\ 
x_2 \\ 
\vdots \\ 
x_{N-1} 
\end{bmatrix}のようになる。 (1 になっているところは、 W_N^0 ともかける)

行列の要素は N^2 あるが、W_N^i には周期性があるので、最大でも N 種類しかパターンはない。
さらに、ここで高速フーリエ変換の条件として、N が 2 の冪乗であることを課すことによって、この行列変換の計算量を減らすことができる。
ざっくり感覚的な話をすると、N が 2 の冪乗の場合、W_N^i には対称性がたくさん見つかる。これを利用して計算を構造化して、O(N^2) を O(N \log N) に落とすことができるというのが FFT の大枠だ。たとえば、N = 8 とすると、i = 0, 8, 16, 24, 32, 40, 48, 56 はすべて 1 になるし、i = 1, 9, 17, 25, 33, 41, 49, 57 はすべて W_N^1になる。同じように、i = 0 \dots 7 ごとにそれぞれ W_N^i の値を持つグループができる。また、i = 0 \dots 3 に対して W_N^{i+4} = - W_N^i となる。たとえば、 W_N^4 = -1 = - W_N^0 になる、という具合だ。なので、最初に書いた行列の要素は正負も無視すると、N / 2 通りだけになる。これをうまい感じにグルーピングして、まとめて計算すると、計算回数が減らせる、というのが基本の考え。

nrnrk

例を見ていけばもっとしっくりくる。

N = 2 の場合、
\begin{bmatrix} 
X_0 \\ 
X_1 
\end{bmatrix} = 
\begin{bmatrix} 
1 & 1 \\ 
1 & -1 
\end{bmatrix} 
\begin{bmatrix} 
x_0 \\ 
x_1 
\end{bmatrix}であり、これは W_2 が 1, -1 の2通りで、正負を無視すれば、絶対値 1 の1通りになる。後々わかるように、これが最小単位になる。

\begin{align*}
X_0 &= x_0 + x_1 = W_2^0 x_0 + W_2^0 x_1 \\
X_1 &= x_0 - x_1 = W_2^0 x_0 + W_2^1 x _1 = W_2^0 x_0 - W_2^0 x _1
\end{align*}
と展開できる。

nrnrk

N = 4 の場合、
\begin{bmatrix} 
X_0 \\ 
X_1 \\ 
X_2 \\ 
X_3 
\end{bmatrix} = 
\begin{bmatrix} 
1 & 1 & 1 & 1 \\ 
1 & -j & -1 & j \\ 
1 & -1 & 1 & -1 \\ 
1 & j & -1 & -j 
\end{bmatrix} 
\begin{bmatrix} 
x_0 \\ 
x_1 \\ 
x_2 \\ 
x_3 
\end{bmatrix}W_4^i のグルーピングの順番に変更すると
\begin{bmatrix} 
X_0 \\ 
X_2 \\ 
X_1 \\ 
X_3 
\end{bmatrix} = 
\begin{bmatrix} 
1 & 1 & 1 & 1 \\ 
1 & -1 & 1 & -1 \\ 
1 & -j & -1 & j \\ 
1 & j & -1 & -j 
\end{bmatrix} 
\begin{bmatrix} 
x_0 \\ 
x_1 \\ 
x_2 \\ 
x_3 
\end{bmatrix}X_0 = x_0 + x_1 + x_2 + x_3X_2 = x_0 - x_1 + x_2 - x_3X_1 = x_0 - jx_1 - x_2 + jx_3X_3 = x_0 + jx_1 - x_2 - jx_3さらに、右辺を先ほどのグルーピングを意識すると、

\begin{align*}
X_0 &= (x_0 + x_2) + (x_1 + x_3) = W_4^0 (x_0 + x_2) + W_4^0 (x_1 + x_3) \\
X_2 &= (x_0 + x_2) - (x_1 + x_3) = W_4^0 (x_0 + x_2) + W_4^2 (x_1 + x_3) = W_4^0 (x_0 + x_2) - W_4^0 (x_1 + x_3) \\
X_1 &= (x_0 - x_2) - j(x_1  - x_3) = W_4^0 (x_0 + W_4^2 x_2) + W_4^1 (x_1 + W_4^2 x_3) = W_4^0 (x_0 - W_4^0 x_2) + W_4^1 (x_1 - W_4^0 x_3) \\
X_3 &= (x_0 - x_2) + j(x_1 - x_3) = W_4^0 (x_0 + W_4^2 x_2) - W_4^1 (x_1 + W_4^2 x_3) = W_4^0 (x_0 - W_4^0 x_2) - W_4^1 (x_1 - W_4^0 x_3)
\end{align*}
となる。実は、W_N^m W_N^n = W_N^{m+n} や W_N^(i+N/2) = - W_N^i を使うと機械的に整理できる。W_4^0, W_4^1 だけで表せる点にも注意。ここで、W_N^i の "通分" を考える。これは W_N の定義から自明だが、i = 0 のときは W_N^0 = 1 であるし、i > 0 かつ N と i の公約数を p とすると、W_N^i = W_{N/p}^{i/p} となる。たとえば、W_4^2 = W_2^1 となる。
通分して整理すると次のようになる。

\begin{align*}
X_0 &= W_2^0 (x_0 + x_2) + W_2^0 (x_1 + x_3) \\
X_2 &= W_2^0 (x_0 + x_2) - W_2^0 (x_1 + x_3) \\
X_1 &= W_4^0 (x_0 - x_2) + W_4^1 (x_1 - x_3) \\
X_3 &= W_4^0 (x_0 - x_2) - W_4^1 (x_1 - x_3)
\end{align*}
掛け算は足し算よりも計算量が大きいことを踏まえて、掛け算の回数を考えると、W_2^0 (x_0 + x_2), W_2^0 (x_1 + x_3), W_4^0 (x_0 - x_2), W_4^1 (x_1 - x_3) の4回(= N回)だけになっている。素朴な離散フーリエ変換の計算の場合は、N^2 必要だったことを考えると削減されているのが明確にわかる。(W_2^0 = 1 なのでそこの掛け算は実質不要でもある)
また、(x_0 + x_2) と (x_1 + x_3) を一つの塊ととらえると、X_0, X_2 の部分は、N=2 のときと同じ構造になっている。この構造によって、計算を分割していき O(N \log N) の計算量が実現される。

この構造を意識して、再度整理をすると、

\begin{align*}
X_0 &= W_2^0 (x_0 + x_2) + W_2^0 (x_1 + x_3) \\
X_2 &= W_2^0 (x_0 + x_2) - W_2^0 (x_1 + x_3) \\
X_1 &= W_2^0 W_4^0 (x_0 - x_2) + W_2^0 W_4^1 (x_1 - x_3) \\
X_3 &= W_2^0 W_4^0 (x_0 - x_2) - W_2^0 W_4^1 (x_1 - x_3)
\end{align*}
となる。これはつまり、(x_0 + x_2), (x_0 - x_2), (x_1 + x_3), (x_1 - x_3) を計算したのち、(x_0 + x_2) と (x_1 + x_3) は N=2の構造で計算をし、(x_0 - x_2) は W_4^0 を (x_1 - x_3) は W_4^1 をかけたのち、それぞれ N=2の構造で計算をすれば良いということになっている。