離散フーリエ変換
離散フーリエ変換(りさんフーリエへんかん、discrete Fourier transform、DFT)とは次式で定義される変換で、フーリエ変換に類似したものであり、信号処理などで離散化されたデジタル信号の周波数解析などによく使われる。また偏微分方程式や畳み込み積分の数値計算を効率的に行うためにも使われる。離散フーリエ変換は(計算機上で)高速フーリエ変換(FFT)を使って高速に計算することができる。
離散フーリエ変換とは、複素関数 f(x)を複素関数F(t)に写す写像であって、次の式で定義されるものを言う。
F(t)= \sum_{x=0}^{N-1} f(x) \exp\left(-i\frac{2 \pi t x}{N} \right)
ここで、Nは任意の自然数、 e はネイピア数、i は虚数単位 (i^2 = -1)で、\piは円周率である。このとき、{x = 0,\dots,N-1}を標本点という。また、この変換を \mathcal{F}_N という記号で表し、
のように略記することが多い。
この逆変換にあたる逆離散フーリエ変換(inverse discrete Fourier transform、IDFT)は
f(x) = \frac{1}{N} \sum_{t=0}^{N-1} F(t) \exp\left(i\frac{2\pi x t}{N} \right)
正規化係数(DFT は 1, IDFT は 1/N)や指数の符号は単なる慣習的なものであり、上式とは異なる式を扱うことがある。DFT と IDFT の差について、それぞれの正規化係数を掛けると 1 / N になることと、指数の符号が異符号であるということがだけが重要であり、根本的には同一の変換作用素と考えられる。DFT と IDFT の正規化係数を両方とも 1/\sqrt{N} にすると、両方ともユニタリ作用素(ユニタリ変換)になる。理論的にはユニタリ作用素にするのが好ましいが、実用上数値計算を行うときは上式のように正規化係数を1つにまとめて、スケーリングを一度に行うことが多い。
離散フーリエ変換はフーリエ変換に類似した変換であるので、フーリエ変換と類似した性質を持つ。
離散フーリエ変換においては、有限個の標本点しか使わないため、ある関数を離散フーリエ変換し、それを逆変換した場合に、標本点以外で元の関数と一致するとは限らない。
すなわち、複素関数fに対して、
F(t)= \sum_{x=0}^{N-1} f(x) \exp\left(-i\frac{2 \pi t x}{N}\right)
により離散フーリエ変換を行い、それを逆変換したものをgとすると
f(x)=g(x) \quad x = 0,\dots,N-1
は言えるが、その他の点でf(x)=g(x)が言えるとは限らない。これを高周波の問題、あるいはエイリアシング(aliasing)という。
h(t,x)=\exp(i2\pi t x/N) は内積
(f,g)=\sum_{n=0}^{N-1} {f(n)}^{*} g(n)
に関し、t,t'が整数のとき直交基底である。
(h(t,x),h(t',x))=\sum_{n=0}^{N-1}
\exp \left(-i2\pi\frac{tn}{N}\right)
\exp \left(i2\pi\frac{t'n}{N}\right)
=N~\delta_{tt'}
\delta_{tt'}はクロネッカーのデルタ
二つの(一般には複素数値の)関数f(x) と g(x) の畳み込みf*gは次のように定義される。
(f*g)(x) = \sum_{n=0}^{N-1} f(n) g(x-n)
ただし、fとgは次のような周期性を持つとする。
\begin{align}
f(x+N)&=f(x) \\
g(x+N)&=g(x)
\end{align}
周期関数の畳み込みを離散フーリエ変換したものは、それぞれの離散フーリエ変換の積になる(畳み込み定理)。
\mathcal{F}_N(f*g)=\mathcal{F}_N(f)\mathcal{F}_N(g)
畳み込み和を直接定義式を用いて計算すると O(N²) の計算量が掛かる。しかし、上式より一旦 DFT をしてから掛算をして、そして IDFT で戻す方法もある。DFT の高速アルゴリズムである FFT を介してこのように計算すると O(N log N) の計算量で済む。
また、二つの(一般には複素数値の)関数f(x)とg(x)の相互相関f\star gは以下のように定義される。
f,gが上記の周期性を持てば、
さらにfの関数値が実数であれば、\mathcal{F}_N(f)(-u)=\overline{\mathcal{F}_N(f)(u)}となる.ここで上線をつけた\bar{z}はzの複素共役を表す。
応用上は、実数値関数を対象とすることが多いが、fが実数値関数であるときには、
(\bar{z}はzの複素共役)。そのため出力(\mathcal{F}_Nf)(t)の半分(t>0)で全ての情報を持っていることになる。
離散フーリエ変換が持つ多くの性質は、\omega_{N} = e^{\frac{i2\pi}{N} が1の原始N乗根(primitive root)であることのみに依存している。そのため、単位元の原始N乗根 \omega_N が存在するならば、複素数以外の環や体においても同様に離散フーリエ変換を定義することができる。離散フーリエ変換 (一般)を参照のこと。
デジタルでは2次元変換が画像の周波数成分を解析するのに使われる。
変換は以下のように定義される。
F(u, v) = \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x,y) e^{- 2 \pi i \left(\frac{u x}{M} +\frac{v y}{N} \right) }
そして逆変換は次のようになる。
f(x, y) = \frac{1}{M N} \sum_{u=0}^{M-1} \sum_{v=0}^{N-1} F(u, v) e^{2 \pi i \left(\frac{u x}{M} +\frac{v y}{N} \right) }
f(x, y)は2次元信号(例えば)であり、fのx列y行成分のことである。
F(u, v)はf(x, y)の2次元周波数スペクトラムである。ここでuはx成分の周波数、vはy成分の周波数である。
2次元DFT は行列を用いて簡単に記述できる。
F = \begin{bmatrix}
F(0, 0) & F(1, 0) & \ldots & F(M-1, 0) \\
F(0, 1) & F(1, 1) & \ldots & F(M-1, 1) \\
\vdots & \vdots & \ddots & \vdots \\
F(0, N-1) & F(1, N-1) & \ldots & F(M-1, N-1)
\end{bmatrix}
W = \begin{bmatrix}
\omega_n^{0 \cdot 0} & \omega_n^{0 \cdot 1} & \ldots & \omega_n^{0 \cdot (N-1)} \\
\omega_n^{1 \cdot 0} & \omega_n^{1 \cdot 1} & \ldots & \omega_n^{1 \cdot (N-1)} \\
\vdots & \vdots & \ddots & \vdots \\
\omega_n^{(N-1) \cdot 0} & \omega_n^{(N-1) \cdot 1} & \ldots & \omega_n^{(N-1) \cdot (N-1)} \\
\end{bmatrix} (注:これはユニタリ行列)
f = \begin{bmatrix}
f(0, 0) & f(1, 0) & \ldots & f(M-1, 0) \\
f(0, 1) & f(1, 1) & \ldots & f(M-1, 1) \\
\vdots & \vdots & \ddots & \vdots \\
f(0, N-1) & f(1, N-1) & \ldots & f(M-1, N-1)
\end{bmatrix}
2次元DFT を行列計算によって以下のように変形できる。
F(u, v) = \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x,y) e^{ -2 \pi i ( \frac{u x}{M} + \frac{v y}{N} ) }
F(u, v) = \sum_{x=0}^{M-1} \left[ \sum_{y=0}^{N-1} f(x,y) e^{ -2 \pi i \frac{v y}{N} } \right] e^{ -2 \pi i \frac{u x}{M} }
F(u, v) = \sum_{x=0}^{M-1} \underbrace{ \left[ \sum_{y=0}^{N-1} f(x,y) e^{ -2 \pi i \frac{v y}{N} } \right] }_{ \mathcal{F}_y f(x, y) } e^{ -2 \pi i \frac{u x}{M} }
F_v = W f
F(u, v) = \sum_{x=0}^{M-1} F_v(x, v) e^{- 2 \pi i \frac{u x}{M} }
F = W F_v^{T}
F = W f^T W
以下上式の 1 - 7 を解説すると、
式1から e-2πiux/M を内側σから出した
内側のσは、信号f(x, y)のyの次元(\mathcal{F}_y\{f(x, y)\}と書いた行)の1次元DFTであることを示している
式3で示した、\mathcal{F}_y\{f(x, y)\}の行列表現
式3の注目箇所をF_v(x, v)で書き変えた - F_v(x, v)はf(x, y)のx行目の行のv番目の周波数。f(x, y)の1次元DFTの行を集めたものがFvである。
式4の1次元DFTの行列表現と同様に、FvTを使って式5を表現した
式4を式6に代入した結果。ここでWは対称行列であるのでW=WTとした。
F_vの行はf(x, y)のx行目の行を1次元DFTしたものである。ゆえにF_vはf(x, y)の各行の1次元DFTした結果の行ベクトルを集めたものになる。F=WFvTにおける、FvTを後から掛ける作用素はF_vの列の1次元DFTしたものと同じと考えて良い。
つまり、2次元DFT(2次元フーリエ変換も同様だが)はf(x, y)を、各行ごとに1次元DFTし、その結果をさらに各列ごとに1次元DFTする事と等価である。ここで、1次元DFTの計算はFFTのアルゴリズムで高速に計算できる。そのため実用上は2次元DFTも、2次元FFTとして計算される。
表中で、
W=e^{-i\frac{2\pi}{N}
DFTは多くの広い分野で利用されている。ここでは、その中の一部を示しているだけに過ぎない。これらの応用は高速フーリエ変換(FFT)とその逆変換(IFFT)で高速に計算できることを前提としていて、定義通りにDFTを計算しているのではない。
信号x(t)を解析するのに使われる。ここでtは時間で[0,T]の範囲をとるものとする。例えば、音声信号の場合は、x(t)は時刻tでの空気の圧力を表現することになる。
この信号はN個の等間隔の点で標本化されて、x0, x1, x2, ... , xN-1 の実数列になる。但し標本化の間隔を Δx(=T/N) とすると xk=x(k Δx) である。これのDFTである f0, ..., fN−1 をFFTで計算できる。ただし標本化定理からこれの半分(Nが偶数とすると、fN/2 + 1, ..., fN−1)は冗長であるので捨てるか無視する。
DFT から得られる |fk|/N は信号の周波数 j/T 成分の振幅の半分の値であり、振幅スペクトルと呼ばれる。また、この偏角である arg(fj) はこの周波数成分の位相を表す。また、|fk|2はパワースペクトルと呼ばれ、この周波数成分のエネルギーを表している。
fkは信号x(t)のフーリエ級数の係数に相当するものと考えることができる。そのために、無限に広がるフーリエ級数計算を、有限のサンプル点に対しての DFT を使って近似するという形になる。連続信号の場合はこれをスペクトル推定(spectral estimation)と呼び、色々な推定法がある。(例えば、DFT が有限サンプル点を扱うことに起因するスペクトル漏れの弊害を少しでも和らげるために窓関数(窓))を使うことがよく行われる。)
信号処理の分野では周波数領域(フーリエ変換)で処理をすることも少なくない。例えば、画像の非可逆圧縮や音声圧縮技術などでは離散フーリエ変換の考えが用いられている。信号に対して DFT (実装上ではFFT) を行い、人間が知覚しづらい周波数成分の情報を取り除くことで、正味の情報量を削減(圧縮)する。最も単純な方法では、IDFTの際にその取り除いた周波数情報(フーリエ係数)を0にすることにより、通常のIDFTを実現する。実装の単純化(計算の効率化)のために、実数演算のみで可能な離散コサイン変換を使って2次元情報を圧縮した例がJPEGである。
大きな数や多項式の乗算アルゴリズムにおいて、DFTを使う高速なアルゴリズムとして1971年にショーンハーゲ・ストラッセン法が発見された。数字や多項式の係数の列を畳み込み演算の対象のベクトルと見なす。するとそれぞれのベクトルの DFT を造り、その結果同士のベクトルを要素ごとに乗算した新たなベクトルを作り、それを逆変換することにより掛算の計算結果が得られる(つまり畳み込み定理を使う)。2007年に理論上はショーンハーゲ・ストラッセン法よりも高速なアルゴリズム(Fürer's algorithm)が発見された。
離散フーリエ変換 F_n=\dfrac{1}{N} \sum^{N-1}_{k=0}f_k e^{i \dfrac{2n \pi}{N}k}
離散フーリエ逆変換 f_k=\sum^{N-1}_{n=0} F_n e^{i \dfrac{2n \pi}{N} k}とするところもある(図解雑学 フーリエ変換)。
これは、フーリエ変換をF(\omega)=\int^{\infty}_{-\infty} f(x) e^{-i \omega x}dx、フーリエ逆変換をf(x)=\dfrac{1}{2 \pi} \int^{\infty}_{-\infty} F(\omega) e^{i \omega x} d\omegaとしたときの式である。