From Fourier to DCT

信号を周波数の世界で見ると、時間領域では見えなかった構造が見える。画像処理でいえば、エッジは高周波、なだらかな領域は低周波。この「周波数で分解して見る」道具がフーリエ変換で、そこから派生したDFT、FFT、DCTが画像圧縮や信号処理の基盤になっている。

スペクトルとフーリエ級数展開

周期信号 $S(t)$ は三角関数の重ね合わせで表現できる。これがフーリエ級数展開。

$$S(t) = a[0] + \sum_{k=1}^{\infty} \left\{ a[k]\cos k\omega t + b[k]\sin k\omega t \right\}$$

$a[0]$ は直流成分（信号の平均値）。$a[1], b[1]$ が基本波成分（周波数 $\omega$ の成分）。$k \geq 2$ の $a[k], b[k]$ が高調波成分。スペクトルはこの各周波数成分の振幅を周波数軸上にプロットしたもの。

複素フーリエ展開にすると見通しが良くなる。

$$S(t) = \sum_{k=-\infty}^{\infty} C[k]\, e^{jk\omega t}$$

$j = \sqrt{-1}$。オイラーの公式 $e^{j\theta} = \cos\theta + j\sin\theta$ で三角関数を複素指数関数にまとめた形。$C[k]$ が複素フーリエ係数で、振幅と位相の情報を同時に持っている。

フーリエ級数の面白いところは、矩形波みたいな不連続な信号でもsin/cosの重ね合わせで近似できること。基本波に高調波を足していくと、だんだん元の信号に近づいていく。

離散フーリエ変換（DFT）

連続信号をコンピュータで扱うにはサンプリング（離散化）が必要。アナログの周期信号をサンプリング周期 $\tau$ で $N$ 個のデータとして取り出す。周期 $T = N\tau$。

$$S[k] = \frac{1}{N}\sum_{i=0}^{N-1} s[i]\, e^{-j\frac{2\pi}{N}ki}$$

これが離散フーリエ変換。連続のフーリエ変換を離散データに適用したもの。$s[i]$ が時間領域のサンプル、$S[k]$ が周波数領域のスペクトル。

DFTでは信号の周期とサンプリング区間を一致させる必要がある。一致してないとスペクトルが漏れる（スペクトルリーケージ）。これを緩和するためにハミング窓やハニング窓みたいな窓関数が使われる。

高速フーリエ変換（FFT）

DFTをそのまま計算すると、データが $N$ 個のとき $N^2$ 回の乗算が必要。$N = 1024$ だと約100万回。これをなんとかしたのがFFT。

FFTのカギは回転因子 $W$ の対称性。

$$W = e^{-j\frac{2\pi}{N}}$$

$W^n$ は複素平面上の単位円を $N$ 等分した点。$N = 8$ なら8点が円上に等間隔に並ぶ。ここで重要な性質が、

$$W^n = -W^{n - N/2}$$

半周回ったら符号が反転する。この対称性を使って、$N$ 点のDFTを $N/2$ 点のDFT 2つに分解できる。再帰的に分割していくと、計算量が $O(N^2)$ から $O(N \log N)$ に落ちる。$N = 1024$ なら100万回が約1万回。100倍速い。

FFTを行うためのデータ個数は $2^n$ 個でなくてはならない。足りなければゼロパディングで埋める。

直交変換

$N$ 個の要素からなる $N$ 個のベクトル $\boldsymbol{\Phi}$ が、

$$\langle \boldsymbol{\Phi}_m, \boldsymbol{\Phi}_n \rangle = \sum_{i=0}^{N-1} \Phi_m[i]\, \Phi_n[i] = \begin{cases} 0 & m \neq n \\ 1 & m = n \end{cases}$$

を満たすとき、$\boldsymbol{\Phi}$ は正規直交基底（orthonormal basis）という。

信号ベクトル $\mathbf{f}$ を基底で展開すると、

$$\mathbf{f} = \boldsymbol{\Phi}\, \mathbf{F} \qquad \Longleftrightarrow \qquad \mathbf{F} = \boldsymbol{\Phi}^T \mathbf{f}$$

$\boldsymbol{\Phi}^{-1} = \boldsymbol{\Phi}^T$ が成り立つとき、$\boldsymbol{\Phi}$ は直交行列と呼ばれる。逆変換が単に転置を取るだけで済むから計算が楽。

特定の係数 $F[n]$ の大きさを求めたい場合は、基底 $\boldsymbol{\Phi}_n$ との内積を取ればいい。直交性のおかげで他の成分はキャンセルされて、$F[n]$ だけが残る。これが「周波数成分を取り出す」操作の数学的な正体。

DFTの基底関数は、

$$\Phi_k[i] = \frac{1}{\sqrt{N}}\, e^{j\frac{2\pi}{N}ki}$$

これを使った変換が離散フーリエ変換。直交変換の一種として位置づけられる。

空間周波数

画像の文脈では、時間ではなく空間方向の周波数を考える。長さ $L$ を周期とする1次元画像の正弦波信号は、

$$f(x) = 0.5 + A\sin\!\left(\frac{2\pi}{L}x\right)$$

空間周波数 $\mu = L^{-1}$ は単位長さあたりの正弦波の個数を表す。

離散化すると、全体の画素数 $M$、周期 $L$、正弦波の山の数 $k$ として、

$$f[m] = 0.5 + A\sin\!\left(k\frac{2\pi}{M}m\right)$$

空間周波数 $k$ を変えると波の細かさが変わる。$k = 1$ はゆるやかな波、$k$ が大きくなるほど細かい波になる。

2次元に拡張すると、

$$f[m, n] = 0.5 + A\sin\!\left(k\frac{2\pi}{M}m + l\frac{2\pi}{N}n\right)$$

$k$ が水平方向の空間周波数、$l$ が垂直方向の空間周波数。画像のフーリエスペクトルで中心から離れた点ほど高い空間周波数に対応するのはこの構造のため。

離散コサイン変換（DCT）

DCTはフーリエ変換の親戚で、コサイン関数だけを基底に使う変換。JPEG圧縮の心臓部。

DCT係数は $C[0]$ がDC係数（直流成分、画像ブロックの平均輝度）、$C[k]$（$k \geq 1$）がAC係数（交流成分、変化の情報）に分かれる。

DFTと比べたDCTの利点は、実数だけで完結すること（虚数部がない）、自然画像のエネルギーが低周波のDC近傍に集中しやすいこと。つまり少数の係数で画像のほとんどの情報を表現できる。$N$ が大きくなるほど直流成分方向の集中度が増して、高周波係数を大胆に量子化（切り捨て）できる。これがJPEG圧縮の原理。

DCTの基底関数はこんな感じ。$k$ が大きくなるほど振動が激しくなる。

1次元のDCTだけど、画像は2次元だから実際には $8 \times 8$ ブロックに対して2次元DCTを適用する。JPEGの場合は画像を $8 \times 8$ のブロックに分割して、各ブロックに2D-DCTをかけて、量子化テーブルで高周波成分を丸めて、エントロピー符号化で圧縮する。DCTそのものは可逆だけど、量子化のステップで情報が失われるから非可逆圧縮になる。

ちなみに、1次元DCTだと時間方向のDCTは $O(N^2)$ だけど、FFTと同様の高速アルゴリズム（FCT: Fast Cosine Transform）で $O(N \log N)$ に落とせる。$N = 8$ くらいだと素朴に計算しても十分速いけど、大きなブロックサイズやリアルタイム処理では効いてくる。