1.1 離散 Fourier 変換

 通常 FFT は,離散 Fourier 変換を高速に行う技法です.この 離散 Fourier 変換は通常の Fourier 変換とは異なり,少し癖のある 変換です.FFT を安全に利用するためにはまず,離散 Fourier 変換の 性質を知っていなければなりません.

1.1.1 DFT の定義

 離散 Fourier 変換 - Discrete Fourier Transform は, 通常の Fourier 変換の無限区間積分を有限の和で書き換えたもので, 時間領域,周波数領域ともに離散化された Fourier 変換のことです. 因みに,Fourier 級数展開は,周波数領域でのみ離散化された変換に 相当します.

 ここでは以後,時間領域の離散データ (a_0, a_1, ..., a_N-1) から 周波数領域への離散データ (A_0, A_1, ..., A_N-1) への 長さ N の DFT を

          N-1       jk           -2πi / N    
    A  =  Σ   a   W   ,   W  = e             
     k    j=0  j   N       N                 

で定義します.この変換の逆変換は

          1  N-1      -jk    
    a  = --- Σ   A   W       
     k    N  j=0  j   N      

となります.

 DFT の定義は,書物によって異なるので注意してください. 少し調べた限りでは,規格化定数 (1/N, sqrt(1/N)) の付け方と, 指数関数の符号のとり方は,考えられるすべての組合せがありました. FFT の使用での多くのトラブルはこのことに起因しているようです.

 一般に,DFT では,負の周波数データ A_-k や,添字が N より大きい データは除外されます.なぜならば,W^N = 1 という周期性から, 負の周波数データ A_-k は A_N-k に一致し,A_N+k は A_k に 一致するので,連続な N 個のデータさえあれば,あとは周期的に 拡張できるからです.そのことから,負の引数のデータを扱う場合に, DFT の定義として,データの添字を 0 から N-1 ではなく, -N/2 から N/2-1 までとすることがよくあります.この場合, データを適当に並べ替える(あるいは符号を付け変える)ことが 必要になります.

 また,対称性をよくするために,DFT を次のように拡張することが あります.

          N-1      (j+δ1)(k+δ2)   
    A  =  Σ   a   W               
     k    j=0  j   N              

定数 δ1, δ2, は,片方を 1/2 にすることが多く,これは,Odd DFT と いわれるものです.この拡張 DFT は,通常の DFT の変換前と変換後の データに適当な係数を掛けたものにすぎないのですが,実対称な データなどの変換(離散コサイン変換など)では,データの対称性がよくなり, FFT アルゴリズムが直接適用できるようになります.

1.1.2 DFT と通常の Fourier 変換

 DFT を実際に応用するときに注意しなければいけないことがあります. それは,DFT と普通の Fourier 変換とは似た性質はあるが,異なるもの であるということです.違いの一つは,DFT は離散的だということです. 連続な関数の Fourier 変換を離散的な DFT で近似するとき,離散化誤差が 発生します.一般に,m 階微分が有限な関数を離散化すると,単位区間の データ数 N に対して,DFT の離散化誤差は,ほぼ N^{-m-1} に比例した 大きさになります.また,離散化する関数が無限階微分可能で,ある条件を 満たすならば,離散化誤差はデータ数 N に対して e^{-C N} となり,N を 少し大きくすると指数関数的に急激に減少します.さらに,離散化する 関数が整関数で,その Fourier 変換の周波数成分が fmax 以上を含まない ならば,単位区間のデータ数 N を 2fmax 以上にすれば,離散化誤差は 完全にゼロになります(いわゆる標本化定理).要するに,関数が十分滑らかで あれば,標本数を多くすれば離散化誤差は十分小さくなるということです. もし,関数が滑らかでなく不連続な場合は,いくら N を大きくしても ほとんど近似はよくならず,普通の Fourier 変換と DFT とはまったく 別物になってしまいます[参考文献].

 DFT と Fourier 変換とのもう一つの違いは,通常の Fourier 変換は 無限区間の積分なのに対して DFT は有限区間だということです. 通常の Fourier 変換を DFT で置き換える場合,積分を有限で打ち切ら なければならず,当然誤差(打ち切り誤差)が発生します.さらに, DFT は離散的な関数の Fourier 級数展開に相当するので,DFT される関数は 単なる長さ N の関数ではなく,周期的に拡張された関数とみなされます. DFT を使って周波数解析を行う場合,長さ N のデータは,この N の 周期で無限に続くデータに改ざんされてしまいます.この周期的に拡張 された関数は一般に,つなぎ目で不連続となるため,DFT には 大きな離散化誤差が発生します.この問題は,窓関数を掛けたりするなどの データを滑らかな周期関数にする操作で大幅に改善できます.また, 打ち切り誤差は,変換する関数が減少関数ならば,区間の長さを十分に 大きくとることで小さくできます.もし,関数の性質がわかっていれば, DFT を保存するような線形の加速法が使える場合もあります.

 離散化の簡単な具体例を次に示します.通常の Fourier 変換

             1   -∞        -iωx        
    F(ω) =  --- ∫   f(x)  e      dx    
            2π   ∞                     

を積分区間 ±T で打ち切り,長さ N の DFT で離散化すると

    〜          Δx      k  N-1                 jk    
    F(Δω・k) =  --- (-1)   Σ   f(-T + Δx・j)  W       
               2π         j=0                N      


          -2πi / N         2T          π     
   W  = e          , Δx = ---- , Δω = ----   
    N                      N           T     

となります.離散的な近似関数 F〜() には離散化誤差と打ち切り誤差の両方が 含まれます.さらに注意しなければいけないのは,F〜(ω) は周期 πN/T の 周期関数になるということです.したがって,F〜(Δω・k) での離散の 間隔 Δω は自動的に π/T となり,近似として有効な ω の範囲は一周期分 のみで -πN/(2T) から +πN/(2T) ということになります.このことから, 離散幅は任意には選べず,いわゆる不確定性関係というものが存在します. 周波数の離散幅 Δω と時間幅 T の積は常に π となり,時間の離散幅 Δx と 有効な周波数幅 πN/(2T) の積も常に π となります.要するに,DFT では短い 時間のデータで細かい周波数の分解能を出すのは不可能で,時間幅か周波数の 分解能のどちらかを犠牲にしなければいけないということです.

1.1.3 DFT の性質

 ここでの離散 Fourier 変換は,複素数データから複素数データへの 変換です.しかし,実際の応用では,実数のデータがほとんどです. そこで,実離散Fourier変換/逆変換の性質を少しだけ考えてみます. もし,a_j が実数ならば,A_k は複素共役対称

           _     
    A    = A     
     N-k    k    

となることが簡単にわかると思います.とくに,A_0, A_N/2 は実数に なります.この変換後のデータは半分が冗長になります.逆に,a_j が 複素共役対称ならば,A_k は実数となります.さらに,a_j が実対称 (a_j = a_N-j)ならば,A_k も実対称となり,a_j が実反対称 (a_j = -a_N-j)ならば,A_k は純虚数で反対称となります.これらは, タイプ気領セ競灰汽ぅ麒儡后のセ競汽ぅ麒儡垢箸い錣譴襪發里任后

Fourier 変換の対称性
変換前 変換後
実数 複素共役対称
純虚数 複素共役反対称
複素共役対称 実数
複素共役反対称 純虚数
対称 対称
反対称 反対称
実対称 実対称
実反対称 純虚数反対称

これらの DFT の対称性は,通常の Fourier 変換での対称性と全く同じです.

 DFT に関連する重要な演算に,畳み込み(Digital Convolution)があります. ここでは,離散データ a_j と h_j の長さ N の巡回畳み込みを

          N-1             
    y  =  Σ   a   h       
     k    j=0  j   k-j    

により定義します.直観的には,a_j は信号で h_j はフィルタの係数に 相当し,y_j はフィルタされた信号に相当します.ここでの畳み込みは 巡回で,a_j と h_j は長さ N で周期的に拡張されています.通常の Fourier 変換の性質と同様に,巡回畳み込みは DFT を行うと単なる積に 変換されます.すなわち,y_j の DFT されたデータ Y_k は a_j, h_j の DFT されたデータ A_k, H_k を用いて

                    
    Y  =  A   H     
     k     k   k    

と表されます.したがって,巡回畳み込みは,Y_k を逆 DFT すれば計算で きることになり,FFT による高速算法が利用できます.この方法は, 巡回畳み込みを計算する最も強力な方法として知られています.しかし, 比較的短い長さの巡回畳み込みの計算には,一次元畳み込みをより短い長さの 多次元畳み込みに変換する直接算法が適している場合があります.この方法は, 後に示す Prime Factor 型 FFT の添字の変換で実現できます [参考文献].


目次へ戻る