本記事は，京都大学人工知能研究会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