離散フーリエ変換

出典: フリー百科事典『ウィキペディア(Wikipedia)』
移動: 案内検索

離散フーリエ変換(りさんフーリエへんかん、英語: discrete Fourier transformDFT)とは離散群上のフーリエ変換であり、信号処理などで離散化されたデジタル信号周波数解析などによく使われる。また偏微分方程式畳み込み積分を効率的に計算するためにも使われる。離散フーリエ変換は(計算機上で)高速フーリエ変換(FFT)を使って高速に計算することができる。

N 個の複素数x_0, ..., x_{N-1} に対して DFT することで N 個の複素数列 X_0, ..., X_{N-1}が得られる。

X_k = \sum_{n=0}^{N-1} x_n e^{-i\frac{2 \pi k n}{N}} \quad \quad k = 0, \dots, N-1

ここで eネイピア数i虚数単位 (i^2 = -1)[1]で、\pi円周率である。また、この変換を \mathcal{F} という記号で表し、x= (x_0, ..., x_{N-1}), X= (X_0, ..., X_{N-1}) とおいて

\mathbf{X} = \mathcal{F}(\mathbf{x})

のように略記することが多い。

この逆変換にあたる逆離散フーリエ変換英語: inverse discrete Fourier transformIDFT)は

x_n = \frac{1}{N} \sum_{k=0}^{N-1} X_k e^{i\frac{2\pi k n}{N}} \quad \quad n = 0,\dots,N-1

正規化係数(DFT は 1, IDFT は 1/N)や指数の符号は単なる慣習的なものであり、上式とは異なる式を扱うことがある。DFT と IDFT の差について、それぞれの正規化係数を掛けると 1 / N になることと、指数の符号が異符号であるということがだけが重要であり、根本的には同一の変換作用素と考えられる。DFT と IDFT の正規化係数を両方とも \frac{1}{\sqrt{N}} にすると、両方ともユニタリ作用素(ユニタリ変換)になる。理論的にはユニタリ作用素にするのが好ましいが、実用上数値計算を行うときは上式のように正規化係数を1つにまとめて、スケーリングを一度に行うことが多い。

性質[編集]

離散フーリエ変換はフーリエ変換を離散群上で定義したものであるので、フーリエ変換と同様の性質を持つ。

完全性[編集]

完全性 DFT は逆変換を持ち、次のような線形写像である。

\mathcal{F}:\mathbf{C}^N \to \mathbf{C}^N

但し、C複素数集合である。つまり、任意の N 次元複素ベクトルは、その DFT と IDFT としてN次元複素ベクトルを持つということである(N ≥ 0 とした)。

直交性[編集]

e^{i2\pi\frac{k n}{N}}直交基底である。

\sum_{n=0}^{N-1}
\left(e^{ i2\pi\frac{kn}{N}}\right)
\left(e^{-i2\pi\frac{k'n}{N}}\right)
=N~\delta_{kk'}

\delta_{kk'}クロネッカーのデルタ

パーセバルの定理[編集]

X_k, Y_k がそれぞれ x_n, y_nの DFT の時、Plancherelの定理から

\sum_{n=0}^{N-1} x_n y^*_n = \frac{1}{N} \sum_{k=0}^{N-1} X_k Y^*_k

ここで * は複素共役を表す。パーセバルの定理は Plancherelの定理の特殊なものであり、

\sum_{n=0}^{N-1} |x_n|^2 = \frac{1}{N} \sum_{k=0}^{N-1} |X_k|^2.

信号処理の観点から見ると、左辺は x_nエネルギーを示している。また、|X_k|^2パワースペクトルと呼ぶように、右辺は全周波数のエネルギーを表しているといえる。つまり DFT の変換前後でのエネルギーの関係を示していると解釈できる。

畳み込み定理と相互相関定理[編集]

x= x_ny = y_n の畳み込みx*y

(\mathbf{x*y})_k = \sum_{n=0}^{N-1} x_n y_{k-n} \quad \quad k = 0,\dots,N-1

ここでyは次のように循環するとする。

y_{-n} = y_{N-n}\quad\quad~~~~~~~~~~ n = 0, ..., N-1

循環畳み込みを DFT すると単なる積になる(畳み込み定理)。つまりz_k = (\mathbf{x*y})_k

Z_k=X_k Y_k \quad \quad~~~~~~~~~~ k = 0,\dots,N-1

となる。ただし大文字のX, Y, Zはそれぞれ小文字のx, y, zの DFT である。なお、DFT の正規化係数が 1 ではない場合は上式に定数係数が付くことになる。

畳み込み和を直接定義式を用いて計算すると O(N²) の計算量が掛かる。しかし、上式より一旦 DFT をしてから掛算をして、そして IDFT で戻す方法もある。DFT の高速アルゴリズムである FFT を介してこのように計算すると O(N log N) の計算量で済む。

xjyk相互相関 zk は以下のように定義される。

z_k=(\mathbf{x\star y})_k = \sum_{n=0}^{N-1}x_n^*\,y_{n+k}

ここでyは既述のように循環するとして、zk の DFT は

Z_k = X_k^*\,Y_k

補間三角多項式との関係[編集]

実数列の DFT[編集]

応用の上ではx0, ..., xN−1実数のことが多いが、このとき、fj = fNj* である(*は複素共役)。そのため出力f0, ..., fN−1の半分で全ての情報を持っていることになる。このとき、直流成分 f0 は純粋に実数で、ナイキスト成分 fN/2 も実数であるので、複素数出力fの最初の半分とナイキスト成分というN個実数列によって DFT は冗長性なく完全に記述ができる。

一般化[編集]

2次元での変換[編集]

デジタル画像処理では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) } \quad \quad u = 0, \dots, M-1; v = 0, \dots, N-1

そして逆変換は次のようになる。

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) } \quad \quad x = 0, \dots, M-1; y = 0, \dots, N-1

但し

  • f(x, y)は2次元信号(例えば画像)であり、fxy行成分のことである。
  • F(u, v)f(x, y)の2次元周波数スペクトラムである。ここでux成分の周波数、vy成分の周波数である。

2次元DFT は行列を用いて簡単に記述できる。

F = W f^T W

ここで

  • 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 を行列計算によって以下のように変形できる。

  1. 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} ) }
  2. 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} }
  3. 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} }
  4. F_v = W f
  5. F(u, v) = \sum_{x=0}^{M-1} F_v(x, v) e^{- 2 \pi i \frac{u x}{M} }
  6. F = W F_v^{T}
  7. F = W f^T W

以下上式の 1 - 7 を解説すると、

  1. 2次元DFTの定義
  2. 式1から e-2πiux/M を内側σから出した
  3. 内側のσは、信号f(x, y)yの次元(\mathcal{F}_y\{f(x, y)\}と書いた行)の1次元DFTであることを示している
  4. 式3で示した、\mathcal{F}_y\{f(x, y)\}の行列表現
  5. 式3の注目箇所をF_v(x, v)で書き変えた - F_v(x, v)f(x, y)x行目の行のv番目の周波数。f(x, y)の1次元DFTの行を集めたものがFvである。
  6. 式4の1次元DFTの行列表現と同様に、FvTを使って式5を表現した
  7. 式4を式6に代入した結果。ここでW対称行列であるのでW=WTとした。

F_vの行はf(x, y)x行目の行を1次元DFTしたものである。ゆえにF_vf(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}}

また、特に指定の無い場合は、

k=0,...,{N-1}

n=0,...,{N-1}

とする。

時間領域 周波数領域 備考
x_n=\frac{1}{N}\sum_{k=0}^{N-1}X_k W^{-kn} X_k=\sum_{n=0}^{N-1}x_n W^{kn} IDFT,DFTのWを使った定義
x_{n} X_{k} 定義
y_{n} Y_{k} 定義
ax_{n}+by_{n} aX_{k}+bY_{k} 線形性
x_n{e^{i\frac{2\pi{mn}}{N}}} X_{k-m} 時間軸変調、周波数軸移動
x_{n-m} W^{mk}X_{k} 時間軸移動(正)
x_{n+m} W^{-mk}X_{k} 時間軸移動(負)
\sum_{m=0}^{N-1}{x_m}{y_{n-m}} X_{k}Y_{k} 時間軸畳み込み、周波数軸積
x_{n}y_{n} {\frac{1}{N}}\sum_{m=0}^{N-1}X_{m}Y_{k-m} 時間軸積、周波数軸畳み込み
x^{*}_{n} X^{*}_{N-k} 時間軸共役、周波数軸反転
x^{*}_{N-n} X^{*}_{k} 時間軸反転、周波数軸共役
real(x_{n}) X^{even}_{k}=\frac{1}{2}(X_k+X^{*}_{N-k}) 時間軸実部、周波数軸偶関数
imag(x_{n}) X^{odd}_k=\frac{1}{2i}(X_k-X^{*}_{N-k}) 時間軸虚部、周波数軸奇関数
x^{even}_n=\frac{1}{2}(x_k+x^{*}_{N-k}) real(X_{k}) 時間軸偶関数、周波数軸実部
x^{odd}_n=\frac{1}{2i}(x_k-x^{*}_{N-k}) imag(X_k) 時間軸奇関数、周波数軸虚部
a^n\, \frac{1-a^N{W}^{kn}}{1-{a}W^k} べき乗
{N-1 \choose n}\, (1+W^{k})^{N-1}\, 組み合わせ

応用[編集]

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年にショーンハーゲ・ストラッセン法より理論的に高速なアルゴリズム(en:Furer's algorithm)が見付かった。

[編集]

  1. ^ ディジタル信号処理分野の古典である、OppenheimとSchaferの著書『ディジタル信号処理』ではiの代わりにj^2 = -1を使用している。

参考文献[編集]

  • E. O. Brigham, The Fast Fourier Transform and Its Applications (Prentice-Hall, Englewood Cliffs, NJ, 1988).
  • A. V. Oppenheim, R. W. Schafer, and J. R. Buck, Discrete-Time Signal Processing (Prentice-Hall, 1999).
  • S. W. Smith, The Scientist and Engineer's Guide to Digital Signal Processing, (California Technical Publishing, San Diego, 2nd edition, 1999)
  • A. V. Oppenheim, R. W. Schafer, 伊達玄訳, ディジタル信号処理 (コロナ社, 1978).

関連項目[編集]