FFTで多項式のかけ算をする

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

https://qiita.com/advent-calendar/2023/kyoto-kaira

# 本記事の概要
[「pythonで学ぶ音声認識」](https://book.impress.co.jp/books/1120101083)の輪読会を行った際に、音声の特徴量を作るときに行う音声の周波数分解でFFT(高速フーリエ変換)が登場しました。そのときに多項式のかけ算をFFTで高速化できたはずという話が出たので、自分なりにまとめてみました。

# 多項式のかけ算
```math
f(x) = a_0 + a_1 x^1 + a_2 x^2 + \cdots a_{N-1} x^{N-1} 
```
```math
g(x) = b_0 + b_1 x^1 + b_2 x^2 + \cdots b_{M-1} x^{M-1}
```
としたときに、これらをかけ合わせた
```math
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$は
```math
c_n = \sum_{i+j=n} a_i b_j
```
です。
愚直にかけ合わせた場合、実装は以下のようになります。
```python
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は離散フーリエ変換を高速化したものです。
離散フーリエ変換は、次の式で定義されます。
```math
y(k) = \sum_{n=0}^{N-1} x(n) e^{-j \frac{2\pi n k}{N}}
```
また、逆離散フーリエ変換は
```math
x(n) = \frac{1}{N}\sum_{k=0}^{N-1} y(k) e^{j \frac{2\pi n k}{N}}
```
と定義されます。
逆離散フーリエ変換の式に離散フーリエ変換の定義を代入すると
```math
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$とすれば
```math
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'}})$$が成立するので、
```math
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)$とすると、
```math
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)$の離散フーリエ変換で求める必要があるのは、
```math
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$
とすると、
```math
f(x) = f_{even}(x^2) + x  f_{odd}(x^2)
```
となります。
よって
```math
f(\zeta_{N}^0), f(\zeta_{N}^1), \cdot \cdot \cdot　, f(\zeta_{N}^{N-1})
```
を求めるには
```math
f_{even}(\zeta_{N}^0), f_{even}(\zeta_{N}^2), \cdot \cdot \cdot　, f_{even}(\zeta_{N}^{2(N-1)})
```
```math
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}}}$とすると、
```math
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})
```
```math
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}}$以降は前半と同じになり、
```math
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})
```
```math
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関数を用いることができます。
```python
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で学ぶ音声認識 機械学習実践シリーズ](https://book.impress.co.jp/books/1120101083)
[FFT（高速フーリエ変換）を完全に理解する話](https://qiita.com/ageprocpp/items/0d63d4ed80de4a35fe79)