離散フーリエ変換

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

離散フーリエ変換とは、複素関数 f ( x ) {\displaystyle f(x)} を複素関数 F ( t ) {\displaystyle F(t)} に写す写像であって、次の式で定義されるものを言う。

F ( k ) = n = 0 N 1 f ( x ) exp ( i 2 π k x n x N x 0 ) {\displaystyle F(k)=\sum _{n=0}^{N-1}f(x)\exp \left(-i{\frac {2\pi kx_{n}}{x_{N}-x_{0}}}\right)}

ここで、 k は波数、N は標本数、 e {\displaystyle e} ネイピア数 i {\displaystyle i} 虚数単位 ( i 2 = 1 {\displaystyle i^{2}=-1} )[注 1]で、 π {\displaystyle \pi } 円周率である。標本番号 n = { 0 , , N 1 } {\displaystyle n=\{0,\dots ,N-1\}} に対応する x n {\displaystyle x_{n}} 標本点という。

また、この変換を F N {\displaystyle {\mathcal {F}}_{N}} という記号で表し、

F = F N ( f ) {\displaystyle F={\mathcal {F}}_{N}(f)}

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

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

f ( x ) = 1 N k = 0 N 1 F ( k ) exp ( i 2 π k x N ) {\displaystyle f(x)={\frac {1}{N}}\sum _{k=0}^{N-1}F(k)\exp \left(i{\frac {2\pi kx}{N}}\right)}

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

性質

離散フーリエ変換はフーリエ変換に類似した変換であるので、フーリエ変換と類似した性質を持つ。

完全性

離散フーリエ変換においては、有限個の標本点しか使わないため、ある関数を離散フーリエ変換し、それを逆変換した場合に、標本点以外で元の関数と一致するとは限らない。

すなわち、複素関数fに対して、

F ( t ) = x = 0 N 1 f ( x ) exp ( i 2 π t x N ) {\displaystyle F(t)=\sum _{x=0}^{N-1}f(x)\exp \left(-i{\frac {2\pi tx}{N}}\right)}

により離散フーリエ変換を行い、それを逆変換したものをgとすると

f ( x ) = g ( x ) x = 0 , , N 1 {\displaystyle f(x)=g(x)\quad x=0,\dots ,N-1}

は言えるが、その他の点で f ( x ) = g ( x ) {\displaystyle f(x)=g(x)} が言えるとは限らない。これを高周波の問題、あるいはエイリアシング(aliasing)という。

選点直交性

h ( t , x ) = exp ( i 2 π t x / N ) {\displaystyle h(t,x)=\exp(i2\pi tx/N)} 内積

( f , g ) = n = 0 N 1 f ( n ) g ( n ) {\displaystyle (f,g)=\sum _{n=0}^{N-1}{f(n)}^{*}g(n)}

に関し、 t , t {\displaystyle t,t'} が整数のとき直交基底である。

( h ( t , x ) , h ( t , x ) ) = n = 0 N 1 exp ( i 2 π t n N ) exp ( i 2 π t n N ) = N   δ t t {\displaystyle (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'}}

δ t t {\displaystyle \delta _{tt'}} クロネッカーのデルタ

畳み込み定理と相互相関定理

二つの(一般には複素数値の)関数 f ( x ) {\displaystyle f(x)} g ( x ) {\displaystyle g(x)} の畳み込み f g {\displaystyle f*g} は次のように定義される。

( f g ) ( x ) = n = 0 N 1 f ( n ) g ( x n ) {\displaystyle (f*g)(x)=\sum _{n=0}^{N-1}f(n)g(x-n)}

ただし、 f {\displaystyle f} g {\displaystyle g} は次のような周期性を持つとする。

f ( x + N ) = f ( x ) g ( x + N ) = g ( x ) {\displaystyle {\begin{aligned}f(x+N)&=f(x)\\g(x+N)&=g(x)\end{aligned}}}

周期関数の畳み込みを離散フーリエ変換したものは、それぞれの離散フーリエ変換の積になる(畳み込み定理)。 つまり

F N ( f g ) = F N ( f ) F N ( g ) {\displaystyle {\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 ) {\displaystyle f(x)} g ( x ) {\displaystyle g(x)} 相互相関 f g {\displaystyle f\star g} は以下のように定義される。

h ( x ) = ( f g ) ( x ) = n = 0 N 1 f ( n ) g ( x + n ) {\displaystyle h(x)=(f\star g)(x)=\sum _{n=0}^{N-1}f(n)g(x+n)}

f , g {\displaystyle f,g} が上記の周期性を持てば、

( F N h ) ( u ) = ( F N ( f g ) ) ( u ) = F N ( f ) ( u ) F N ( g ) ( u ) {\displaystyle ({\mathcal {F}}_{N}h)(u)=({\mathcal {F}}_{N}(f\star g))(u)={\mathcal {F}}_{N}(f)(-u){\mathcal {F}}_{N}(g)(u)}

さらに f {\displaystyle f} の関数値が実数であれば、 F N ( f ) ( u ) = F N ( f ) ( u ) ¯ {\displaystyle {\mathcal {F}}_{N}(f)(-u)={\overline {{\mathcal {F}}_{N}(f)(u)}}} となる.ここで上線をつけた z ¯ {\displaystyle {\bar {z}}} z {\displaystyle z} の複素共役を表す。

補間三角多項式との関係

実数値関数の離散フーリエ変換

応用上は、実数値関数を対象とすることが多いが、 f {\displaystyle f} が実数値関数であるときには、

( F N f ) ( t ) = ( F N f ) ( t ) ¯ {\displaystyle ({\mathcal {F}}_{N}f)(-t)={\overline {({\mathcal {F}}_{N}f)(t)}}}

( z ¯ {\displaystyle {\bar {z}}} z {\displaystyle z} 複素共役)。そのため出力 ( F N f ) ( t ) {\displaystyle ({\mathcal {F}}_{N}f)(t)} の半分( t > 0 {\displaystyle t>0} )で全ての情報を持っていることになる。

一般化

離散フーリエ変換が持つ多くの性質は、 ω N = e i 2 π N {\displaystyle \omega _{N}=e^{\frac {i2\pi }{N}}} が1の原始N乗根(primitive root)であることのみに依存している。そのため、単位元の原始N乗根 ω N {\displaystyle \omega _{N}} が存在するならば、複素数以外の環や体においても同様に離散フーリエ変換を定義することができる。離散フーリエ変換 (一般)を参照のこと。

2次元での変換

デジタル画像処理では2次元変換が画像の周波数成分を解析するのに使われる。

変換は以下のように定義される。

F ( u , v ) = x = 0 M 1 y = 0 N 1 f ( x , y ) e 2 π i ( u x M + v y N ) {\displaystyle F(u,v)=\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}f(x,y)e^{-2\pi i\left({\frac {ux}{M}}+{\frac {vy}{N}}\right)}}

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

f ( x , y ) = 1 M N u = 0 M 1 v = 0 N 1 F ( u , v ) e 2 π i ( u x M + v y N ) {\displaystyle f(x,y)={\frac {1}{MN}}\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}F(u,v)e^{2\pi i\left({\frac {ux}{M}}+{\frac {vy}{N}}\right)}}

但し

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

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

F = W f T W {\displaystyle F=Wf^{T}W}

ここで

  • F = [ F ( 0 , 0 ) F ( 1 , 0 ) F ( M 1 , 0 ) F ( 0 , 1 ) F ( 1 , 1 ) F ( M 1 , 1 ) F ( 0 , N 1 ) F ( 1 , N 1 ) F ( M 1 , N 1 ) ] {\displaystyle 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 = [ ω n 0 0 ω n 0 1 ω n 0 ( N 1 ) ω n 1 0 ω n 1 1 ω n 1 ( N 1 ) ω n ( N 1 ) 0 ω n ( N 1 ) 1 ω n ( N 1 ) ( N 1 ) ] {\displaystyle 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 = [ f ( 0 , 0 ) f ( 1 , 0 ) f ( M 1 , 0 ) f ( 0 , 1 ) f ( 1 , 1 ) f ( M 1 , 1 ) f ( 0 , N 1 ) f ( 1 , N 1 ) f ( M 1 , N 1 ) ] {\displaystyle 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 ) = x = 0 M 1 y = 0 N 1 f ( x , y ) e 2 π i ( u x M + v y N ) {\displaystyle F(u,v)=\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}f(x,y)e^{-2\pi i({\frac {ux}{M}}+{\frac {vy}{N}})}}
  2. F ( u , v ) = x = 0 M 1 [ y = 0 N 1 f ( x , y ) e 2 π i v y N ] e 2 π i u x M {\displaystyle F(u,v)=\sum _{x=0}^{M-1}\left[\sum _{y=0}^{N-1}f(x,y)e^{-2\pi i{\frac {vy}{N}}}\right]e^{-2\pi i{\frac {ux}{M}}}}
  3. F ( u , v ) = x = 0 M 1 [ y = 0 N 1 f ( x , y ) e 2 π i v y N ] F y f ( x , y ) e 2 π i u x M {\displaystyle F(u,v)=\sum _{x=0}^{M-1}\underbrace {\left[\sum _{y=0}^{N-1}f(x,y)e^{-2\pi i{\frac {vy}{N}}}\right]} _{{\mathcal {F}}_{y}f(x,y)}e^{-2\pi i{\frac {ux}{M}}}}
  4. F v = W f {\displaystyle F_{v}=Wf}
  5. F ( u , v ) = x = 0 M 1 F v ( x , v ) e 2 π i u x M {\displaystyle F(u,v)=\sum _{x=0}^{M-1}F_{v}(x,v)e^{-2\pi i{\frac {ux}{M}}}}
  6. F = W F v T {\displaystyle F=WF_{v}^{T}}
  7. F = W f T W {\displaystyle F=Wf^{T}W}

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

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

F v {\displaystyle F_{v}} の行は f ( x , y ) {\displaystyle f(x,y)} x行目の行を1次元DFTしたものである。ゆえに F v {\displaystyle F_{v}} f ( x , y ) {\displaystyle f(x,y)} の各行の1次元DFTした結果の行ベクトルを集めたものになる。F=WFvTにおける、FvTを後から掛ける作用素は F v {\displaystyle F_{v}} の列の1次元DFTしたものと同じと考えて良い。

つまり、2次元DFT(2次元フーリエ変換も同様だが)は f ( x , y ) {\displaystyle f(x,y)} を、各行ごとに1次元DFTし、その結果をさらに各列ごとに1次元DFTする事と等価である。ここで、1次元DFTの計算はFFTアルゴリズムで高速に計算できる。そのため実用上は2次元DFTも、2次元FFTとして計算される。

離散フーリエ変換表

表中で、

W = e i 2 π N {\displaystyle W=e^{-i{\frac {2\pi }{N}}}}

時間領域 周波数領域 備考
f ( x ) = 1 N k = 0 N 1 F ( k ) W k x {\displaystyle f(x)={\frac {1}{N}}\sum _{k=0}^{N-1}F(k)W^{-kx}} F ( t ) = n = 0 N 1 f ( n ) W t n {\displaystyle F(t)=\sum _{n=0}^{N-1}f(n)W^{tn}} IDFT,DFTのWを使った定義
f ( x ) {\displaystyle f(x)} F ( t ) {\displaystyle F(t)} 定義
g ( x ) {\displaystyle g(x)} G ( t ) {\displaystyle G(t)} 定義
a f ( x ) + b g ( x ) {\displaystyle af(x)+bg(x)} a F ( t ) + b G ( t ) {\displaystyle aF(t)+bG(t)} 線形性
f ( x ) e i 2 π m x N {\displaystyle f(x){e^{i{\frac {2\pi {mx}}{N}}}}} F ( t m ) {\displaystyle F(t-m)} 時間軸変調、周波数軸移動
f ( x a ) {\displaystyle f(x-a)} W a t F ( t ) {\displaystyle W^{at}F(t)} 時間軸移動(正)
f ( x + a ) {\displaystyle f(x+a)} W a t F ( t ) {\displaystyle W^{-at}F(t)} 時間軸移動(負)
m = 0 N 1 f ( m ) g ( x m ) {\displaystyle \sum _{m=0}^{N-1}f(m)g(x-m)} F ( t ) G ( t ) {\displaystyle F(t)G(t)} 時間軸畳み込み、周波数軸積
f ( x ) g ( x ) {\displaystyle f(x)g(x)} 1 N m = 0 N 1 F ( m ) G ( t m ) {\displaystyle {\frac {1}{N}}\sum _{m=0}^{N-1}F(m)G(t-m)} 時間軸積、周波数軸畳み込み
f ( x ) ¯ {\displaystyle {\overline {f(x)}}} F ( N t ) ¯ {\displaystyle {\overline {F(N-t)}}} 時間軸共役、周波数軸反転
f ( N x ) ¯ {\displaystyle {\overline {f(N-x)}}} F ( t ) ¯ {\displaystyle {\overline {F(t)}}} 時間軸反転、周波数軸共役
R e ( f ( x ) ) {\displaystyle Re(f(x))} 1 2 ( F ( t ) + F ( N t ) ¯ ) {\displaystyle {\frac {1}{2}}(F(t)+{\overline {F(N-t)}})} 時間軸実部、周波数軸偶関数
I m ( f ( x ) ) {\displaystyle Im(f(x))} 1 2 i ( F ( t ) F ( N t ) ¯ ) {\displaystyle {\frac {1}{2i}}(F(t)-{\overline {F(N-t)}})} 時間軸虚部、周波数軸奇関数
1 2 ( f ( x ) + f ( N x ) ¯ ) {\displaystyle {\frac {1}{2}}(f(x)+{\overline {f(N-x)}})} R e ( F ( t ) ) {\displaystyle Re(F(t))} 時間軸偶関数、周波数軸実部
1 2 i ( f ( x ) f ( N x ) ¯ ) {\displaystyle {\frac {1}{2i}}(f(x)-{\overline {f(N-x)}})} I m ( F ( t ) ) {\displaystyle Im(F(t))} 時間軸奇関数、周波数軸虚部
a x {\displaystyle a^{x}\,} 1 a N W t x 1 a W t {\displaystyle {\frac {1-a^{N}{W}^{tx}}{1-{a}W^{t}}}} べき乗

応用

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である。

偏微分方程式

大きな数の掛け算の計算

詳細は「ショーンハーゲ・ストラッセン法」および「en:Fürer's algorithm」を参照

大きな数や多項式乗算アルゴリズム(英語版)において、DFTを使う高速なアルゴリズムとして1971年にショーンハーゲ・ストラッセン法が発見された。数字や多項式の係数の列を畳み込み演算の対象のベクトル[要曖昧さ回避]と見なす。するとそれぞれのベクトルの DFT を造り、その結果同士のベクトルを要素ごとに乗算した新たなベクトルを作り、それを逆変換することにより掛算の計算結果が得られる(つまり畳み込み定理を使う)。2007年に理論上はショーンハーゲ・ストラッセン法よりも高速なアルゴリズム(en:Fürer's algorithm)が発見された。

別な定義

離散フーリエ変換  F n = 1 N k = 0 N 1 f k e i 2 n π N k {\displaystyle F_{n}={\dfrac {1}{N}}\sum _{k=0}^{N-1}f_{k}e^{-i{\dfrac {2n\pi }{N}}k}}

離散フーリエ逆変換  f k = n = 0 N 1 F n e i 2 n π N k {\displaystyle f_{k}=\sum _{n=0}^{N-1}F_{n}e^{i{\dfrac {2n\pi }{N}}k}} とするところもある(図解雑学 フーリエ変換)。

これは、フーリエ変換を F ( ω ) = f ( x ) e i ω x d x {\displaystyle F(\omega )=\int _{-\infty }^{\infty }f(x)e^{-i\omega x}dx} 、フーリエ逆変換を f ( x ) = 1 2 π F ( ω ) e i ω x d ω {\displaystyle f(x)={\dfrac {1}{2\pi }}\int _{-\infty }^{\infty }F(\omega )e^{i\omega x}d\omega } としたときの式である。

脚注

[脚注の使い方]
  1. ^ ディジタル信号処理分野の古典である、OppenheimとSchaferの著書『ディジタル信号処理』では i {\displaystyle i} の代わりに j 2 = 1 {\displaystyle 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)
  • W. L. Briggs and van E. Henson: The DFT - An Owner's Manual for the Discrete Fourier Transform (SIAM, 1995).
  • A. V. Oppenheim, R. W. Schafer, 伊達玄訳, ディジタル信号処理 (コロナ社, 1978).

関連項目

理論
サブフィールド
テクニック
標本化