本記事は，京都大学人工知能研究会KaiRAのAdvent Calender 11日目の記事です．

本記事の概要

「pythonで学ぶ音声認識」の輪読会を行った際に、音声の特徴量を作るときに行う音声の周波数分解でFFT(高速フーリエ変換)が登場しました。そのときに多項式のかけ算をFFTで高速化できたはずという話が出たので、自分なりにまとめてみました。

多項式のかけ算

f(x) = a_0 + a_1 x^1 + a_2 x^2 + \cdots a_{N-1} x^{N-1}

g(x) = b_0 + b_1 x^1 + b_2 x^2 + \cdots b_{M-1} x^{M-1}

としたときに、これらをかけ合わせた

h(x) = f(x) g(x) = c_0 + c_1x^1 + c_2 x^2 + \cdots c_{N+M-2} x^{N+M-2}

を求めることを考えます。
ここで$c_n$は

c_n = \sum_{i+j=n} a_i b_j

です。
愚直にかけ合わせた場合、実装は以下のようになります。

def multiply(f, g):
    h = [0] * (len(f) + len(g) -1)
    for i in range(len(f)):
        for j in range(len(g)):
            h[i+j] += f[i]*g[j]

    return h

計算量は$f(x), g(x)$の次数を$N, M$とすると、$O(NM)$となります。
これをFFTで高速化しましょう

離散フーリエ変換

FFTは離散フーリエ変換を高速化したものです。
離散フーリエ変換は、次の式で定義されます。

y(k) = \sum_{n=0}^{N-1} x(n) e^{-j \frac{2\pi n k}{N}}

また、逆離散フーリエ変換は

x(n) = \frac{1}{N}\sum_{k=0}^{N-1} y(k) e^{j \frac{2\pi n k}{N}}

と定義されます。
逆離散フーリエ変換の式に離散フーリエ変換の定義を代入すると

x(n) = \frac{1}{N}\sum_{k=0}^{N-1} (\sum_{t=0}^{N-1} x(t) e^{-j \frac{2\pi t k}{N}}) e^{j \frac{2\pi n k}{N}}

となります。
ここでさきほどの$h(x)$についてこの計算を行いましょう
$N'=N+M-1$とすれば

c_n = \frac{1}{N'}\sum_{k=0}^{N'-1}(\sum_{t=0}^{N'-1} c_t e^{-j \frac{2\pi t k}{N'}}) e^{j \frac{2\pi n k}{N'}}

となります。この括弧の中は$h(e^{-j\frac{2\pi k}{N'}})$になっています。
ここで$h(x) = f(x)g(x)$より、$$h(e^{-j\frac{2 \pi k}{N'}})=f(e^{-j\frac{2 \pi k}{N'}})g(e^{-j\frac{2 \pi k}{N'}})$$が成立するので、

c_n = \frac{1}{N'}\sum_{k=0}^{N'-1}(\sum_{t=0}^{N'-1} a_t e^{-j \frac{2\pi t k}{N'}})(\sum_{t'=0}^{N'-1} b_{t'} e^{-j \frac{2\pi t' k}{N'}}) e^{j \frac{2\pi n k}{N'}}

となります。このそれぞれの括弧の中身は$f(x), g(x)$の離散フーリエ変換になっています。
よって、$h(x)$の各係数を求めるには$f(x), g(x)$をそれぞれ離散フーリエ変換し、各係数をかけ合わせ、それを逆離散フーリエ変換すればよいことがわかりました。

FFT

$\zeta_{N} = e^{-j\frac{2\pi}{N}}$とします。
$f(x) = a_0 + a_1 x^1 + a_2 x^2 + \cdots a_{N-1} x^{N-1} $を離散フーリエ変換したものを$F(k)$とすると、

F(k) = \sum_{n=0}^{N-1} a_n e^{-j \frac{2\pi n k}{N}} 
    = \sum_{n=0}^{N-1} a_n (e^{-j \frac{2\pi n}{N}})^k = f(\zeta_{N}^n)

となります。これを$k=0 \cdot \cdot \cdot N-1$で求めるので、
関数$f(x)$の離散フーリエ変換で求める必要があるのは、

f(\zeta_{N}^0), f(\zeta_{N}^1), \cdot \cdot \cdot　, f(\zeta_{N}^{N-1})

です。ここで$f(x) = a_0 + a_1 x^1 + a_2 x^2 + \cdots a_{N-1} x^{N-1} $を偶数次と奇数次に分けます。すなわち
$f_{even}(x) = a_0 + a_2 x + \cdot \cdot \cdot$
$f_{odd}(x) = a_1 + a_3 x + \cdot \cdot \cdot$
とすると、

f(x) = f_{even}(x^2) + x  f_{odd}(x^2)

となります。
よって

f(\zeta_{N}^0), f(\zeta_{N}^1), \cdot \cdot \cdot　, f(\zeta_{N}^{N-1})

を求めるには

f_{even}(\zeta_{N}^0), f_{even}(\zeta_{N}^2), \cdot \cdot \cdot　, f_{even}(\zeta_{N}^{2(N-1)})

f_{odd}(\zeta_{N}^0), f_{odd}(\zeta_{N}^2), \cdot \cdot \cdot　, f_{odd}(\zeta_{N}^{2(N-1)})

がわかればよいです。
ここで$\zeta_N^2=e^{-j\frac{4\pi}{N}} = e^{-j\frac{2\pi}{\frac{N}{2}}}$です。
よって$\zeta_{\frac{N}{2}} = e^{-j\frac{2\pi}{\frac{N}{2}}}$とすると、

f_{even}(\zeta_{\frac{N}{2}}^0), f_{even}(\zeta_{\frac{N}{2}}^1), \cdot \cdot \cdot　, f_{even}(\zeta_{\frac{N}{2}}^{N-1})

f_{odd}(\zeta_{\frac{N}{2}}^0), f_{odd}(\zeta_{\frac{N}{2}}^1), \cdot \cdot \cdot　, f_{odd}(\zeta_{\frac{N}{2}}^{N-1})

が求まれば良いです。
ここで$\zeta_{\frac{N}{2}}^{\frac{N}{2}}=1$であるので
これらの式の$\zeta_{\frac{N}{2}}^{\frac{N}{2}}$以降は前半と同じになり、

f_{even}(\zeta_{\frac{N}{2}}^0), f_{even}(\zeta_{\frac{N}{2}}^1), \cdot \cdot \cdot　, f_{even}(\zeta_{\frac{N}{2}}^{\frac{N}{2}-1})

f_{odd}(\zeta_{\frac{N}{2}}^0), f_{odd}(\zeta_{\frac{N}{2}}^1), \cdot \cdot \cdot　, f_{odd}(\zeta_{\frac{N}{2}}^{\frac{N}{2}-1})

の偶数次と奇数次で$\frac{N}{2}$個ずつの値を求めれば良いことがわかります。これを繰り返していくと毎回大きさが2分の1になるので大きさ$N$に対し$\log N$回行えばよいことがわかります。一回で$N$個の値を計算しているので計算量は$O(N\log N)$です。

Pythonでの実装

FFTはnumpyライブラリのfft関数を用いることができます。また逆変換にはifft関数を用いることができます。

import numpy as np

def multiply(f, g):
    #fft_sizeは2のべき乗
    fft_size = 1
    while fft_size < len(f) + len(g) - 1:
        fft_size *= 2
    #f(x), g(x)を離散フーリエ変換
    fft_f = np.fft.fft(f, fft_size)
    fft_g = np.fft.fft(g, fft_size)
    #離散フーリエ変換したf(x)とg(x)をかけ合わせる
    fft_h = fft_f * fft_g
    #逆離散フーリエ変換をする
    h = np.fft.ifft(fft_h, fft_size)

    return h.real

参考

Pythonで学ぶ音声認識機械学習実践シリーズ
 FFT（高速フーリエ変換）を完全に理解する話