2の累乗専用のFFTを用いて任意長FFTを実装 ： チャープZ変換

概要

高速フーリエ変換（FFT）とは、離散フーリエ変換（DFT）を高速に計算するアルゴリズム。データの波を扱ううえで非常に重要なアルゴリズムなので広く実装されてきた。そのほとんどは「長さ N=N₁N₂ のDFTは、長さ N₁ のDFTを N₂ 回と長さ N₂ のDFTを N₁ 回で計算できる」というCooley-Tukey型FFTに基づいているため、合成数にしか効果が無く、さらに簡単のため扱えるデータ長が2の累乗に制限されていることも多い。

データ長を2の累乗に整えられるならこれで十分だが、そうもいかない場合も多い。データを切り落としたりゼロ詰めして長さ調整した場合、本来のDFTで計算した場合と異なる結果になってしまい、厳密さが必要な場面では問題となる。

現実には任意長FFTが計算できるライブラリがあるのでそれを使えばいいが、そうもいかないときは自力で実装しなければならない。しかしFFTを全て自力でというのは大変な作業となる。そこでチャープZ変換という手法を使い、任意長DFTを好きな長さの畳み込みに式変形してFFTライブラリで計算することで、ラッパー作成程度の感覚で任意長FFTを実装できる。真面目に実装したものと比べれば計算は遅いものの、きちんと O(N log N) で動くため長いデータに対してはDFTより速い。

仕組み

DFTの定義式

正規化の取り扱いを除けば、DFTは下式で定義される。これは NxN の行列と Nx1 のベクトルの積であり、時間計算量は O(N²) である。

X(k) = \sum_{j=0}^{N-1} x(j) W_{N}^{jk}
, \quad W_{N} = e^{-2\pi i/N}

当たり前だが、どんなFFTもこれを式変形して導かれる。

畳み込みへの式変形

jk という積を変形すると、回転因子 W が分解される。

\begin{align}
jk &= \frac{j^2}{2} + \frac{k^2}{2} - \frac{\left( k - j \right)^2}{2} \\
W_{N}^{jk} &= W_{2N}^{j^2} W_{2N}^{k^2} W_{2N}^{-(k-j)^2}
\end{align}

これを用いてDFTの定義式を変形すると、シグマ記号の部分が畳み込み演算の形になる。

X(k) = W_{2N}^{k^2} \sum_{j=0}^{N-1} \left( x(j) W_{2N}^{j^2} \right) \overline{W_{2N}^{(k-j)^2}}

ちなみに W の肩が時間（添字）の2乗で増えていくようになったが、こんなふうに時間とともに周波数が増えていく信号を「チャープ信号」と言うらしい。また、この式変形はDFT¹を拡張した「Z変換」というものにも適用できるらしい。（信号処理のことはよく知らないので自信は無い）

畳み込みの高速な計算

離散的な畳み込みは簡単に言ってしまえば「掛け算の筆算」である。例えば各桁が a[i], b[i] で表される整数の積の k 桁目を求めるところに着目すれば下式のようになる（配列の添字のため0始まり）。

c_k = \sum_{j=0}^{N-1} a_j b_{k-j} = a_0 b_k + a_1 b_{k-1} + \cdots + a_{N-1} b_{k-N+1}

畳み込みの性質に畳み込み定理という重要なものがある：元の数列同士を畳み込みした結果をフーリエ変換したものは、フーリエ変換したもの同士を単純に各項ごとに掛けたものと一致する²。これを利用すれば、素直に O(N²) で畳み込み演算をするよりも、FFTでフーリエ変換・各項の掛け算・FFTで逆フーリエ変換という手順のほうが O(N log N) で高速になる³。

ただしフーリエ変換を使った畳み込みは通常の掛け算の筆算と異なり、元の範囲外の項は循環して別の項に足し合わされてしまう。それを防ぐために、元の数列に0の項を追加しておき、循環してきた部分の値を0にする必要がある。これさえ満たせばデータの長さを好きに選んでいいので、2の累乗専用のFFTでも計算できる。

実装

配列の扱いが楽なFortranで実装してみる。（多次元配列ではC言語などと逆に左側の添字が先に回ることに注意）

2の累乗専用のFFT

ライブラリを用意できるなら不要。

fft_2x.f90

subroutine fft_2x(ary, n)
    implicit none

    complex(8), intent(inout) :: ary(0:n-1)
    integer   , intent(in)    :: n

    real(8), parameter :: PI = atan(1.0d0) * 4
    integer :: j
    integer :: n1, n2, n3
    complex(8) :: w(0:n-1)
    complex(8) :: tmp(0:n-1)

    w(:) = (/ (exp(cmplx(0, -2 * PI * j / n)), j = 0, n-1) /)

    n1 = 1
    n2 = n
    do while (n2 > 1)
        n3 = 2   ! n2 の最小の素因数（→2で固定）
        n2 = n2 / n3
        call butterfly(ary, tmp, w, n1, n2, n3)
        ary(:) = tmp(:)
        n1 = n1 * n3
    end do

contains

    subroutine butterfly(ary, tmp, w, n1, n2, n3)
        implicit none

        complex(8), intent(in)  :: ary(0:n1-1,0:n2-1,0:n3-1)
        complex(8), intent(out) :: tmp(0:n1-1,0:n3-1,0:n2-1)
        complex(8), intent(in)  :: w(0:n1*n2*n3-1)
        integer   , intent(in)  :: n1, n2, n3

        integer :: j1, j2, j3

        select case (n3)
        case (2)
            tmp(:,0,:) = ary(:,:,0) + ary(:,:,1)   ! (X[0]) = (1, 1) * (x[0])
            tmp(:,1,:) = ary(:,:,0) - ary(:,:,1)   ! (X[1])   (1,-1)   (x[1])
        end select

        do j2 = 1, n2-1
            do j3 = 1, n3-1
                tmp(:,j3,j2) = tmp(:,j3,j2) * w(j2*j3*n1)
            end do
        end do
    end subroutine
end subroutine

任意長FFT

仕組みの節とコードのコメントの通りなので、あまり説明することは無い。

fft.f90

subroutine fft(ary, n)
    implicit none

    complex(8), intent(inout) :: ary(0:n-1)
    integer   , intent(in)    :: n

    real(8), parameter :: PI = atan(1.0d0) * 4
    integer :: j, j2
    integer :: m
    complex(8) :: chirp_w(0:n-1)
    complex(8), allocatable :: a(:), b(:), c(:)

    ! nが2の累乗の場合はそのまま解く
    if (iand(n, n-1) == 0) then
        write(0, *) "n is a power of 2.  Call fft_2x() directly."
        call fft_2x(ary, n)
        return
    end if

    !--- Bluestein's FFT algorithm (chirp Z-transform) ---!

    ! chirp twiddle factors W_{2n}^{j^2}
    j2 = 0
    do j = 0, n-1
        chirp_w(j) = exp(cmplx(0, -2 * PI * j2 / 2 / n))
        j2 = modulo(j2 + 2*j + 1, 2*n)
    end do

    ! m >= 2n-1 を満たす2の累乗 m を見つける
    m = 1
    do while (m < 2 * n - 1)
        m = m * 2
    end do

    allocate(a(0:m-1))
    allocate(b(0:m-1))
    allocate(c(0:m-1))

    ! ゼロ詰めしてFFTで畳み込み演算
    a(0:n-1) = ary(0:n-1) * chirp_w(:)
    a(n:m-1) = 0
    b(0:n-1) = conjg(chirp_w(:))
    b(n:m-n) = 0
    b(m-n+1:m-1) = conjg(chirp_w(n-1:1:-1))
    call fft_2x(a, m)
    call fft_2x(b, m)
    c(:) = a(:) * b(:) / m   ! 正規化は用いるFFTライブラリに応じて修正
    c(:) = conjg(c(:))
    call fft_2x(c, m)   ! 逆変換を正変換で代用する場合、前後にデータの共役をとる補正が必要
    c(:) = conjg(c(:))

    ary(:) = c(0:n-1) * chirp_w(:)
end subroutine

テスト

フーリエ変換の厳密解が明らかな、周波数 k の正弦波（複素数）をデータに与えてFFTを使ってみる。

fft_test.f90

program fft_test
    implicit none

    real(8), parameter :: PI = atan(1.0d0) * 4
    integer :: n
    integer :: j, k, jk
    complex(8), allocatable :: ary(:)

    write(0, *) "Input samples n:"
    read(*, *) n
    write(0, *) "Input frequency k:"
    read(*, *) k

    allocate(ary(0:n-1))
    jk = 0
    do j = 0, n-1
        ary(j) = exp(cmplx(0, 2 * PI * jk / n))
        jk = modulo(jk + k, n)
    end do

    call fft(ary, n)

    ary(k) = ary(k) - n   ! 厳密解からの差を計算する
    write(*, *) "RMS Error = ", sqrt(sum(real(ary(:)) ** 2 + aimag(ary(:)) ** 2) / n)
end program

console

$ f95 *.f90
$ ./a.out
 Input samples n:
123456
 Input frequency k:
789
 RMS Error =    1.0562400625613464E-004

6桁もある答えが 10^-4 の精度で一致しているので、計算は合っていると思う。

速度の見積もり

FFTの時間計算量は O(N log N) だが、チャープZ変換を用いた場合、

• m ≧ 2n - 1 を満たす2の累乗 m は、最小のものは 2n ≦ m < 4n の範囲にある
• 畳み込みの計算には、FFTを3回実行する必要がある⁴

ので、 n がちょうど2の累乗のときに比べて単純計算で 6 ～ 12倍＋α の時間がかかる。また、時間は n でなく m の大きさに依存するので、 n の大きさが多少異なっても時間は変わらず、2の累乗を境に倍増するという階段状のグラフになると予想できる。

用いるFFTライブラリが2の累乗以外（3や5といった素因数）にも対応していれば、より小さい m を見つけることで高速化できるかもしれない（が結局2の累乗のほうが単純で速い可能性もある）。

---

1. より厳密には、無限長のデータを扱う離散時間フーリエ変換（DTFT） ↩

2. フーリエ変換の正規化を見落とすと定数倍ずれるので注意 ↩

3. 畳み込みを掛け算の筆算と言ったが、実際に巨大な整数同士を掛ける高速なアルゴリズムでは畳み込み定理を利用する ↩

4. もしも同じ長さのFFTを何度も実行するなら、配列 b(:) を事前にFFTしておき使い回すことで2回に減らせる ↩