フーリエ展開による畳み込み積分の高速化

フーリエ展開というと、「複数の波が合わさったデータを、周波数が1の成分、2の成分、…に分解してくやつでしょ」ってイメージだけど、それだけでなく、畳み込み積分に使えるらしい。

畳み込み積分

合成積（畳み込み）の意味と応用３つ | 高校数学の美しい物語

わかりやすさのため、1次元の離散的な畳み込みについての例を示す。（実際は連続的な関数でもよい）

2つの数列 $A,B$ があったとする。前から、 $0,1,2,...$ と添字がついているとする。

$\begin{align} A &= \{5, 2, 3, 8, 1\} \\ B &= \{4, 7, 6, 0\} \end{align}$

この2つを合成して、新たな1つの数列 $C=A \ast B$ を作る。

合成の仕方は、 $C=\{c_0,c_1,...\}$ を0で初期化し、全ての $(a_i, b_j)$ の組について $c_{i+j}$ に $a_i b_j$ を加算していく。

$c_k$ を主眼に言い換えると、 $k=i+j$ となるような全ての $(a_i,b_j)$ の組の $a_ib_j$ の和となる。 $\displaystyle c_k = \sum_{i=0}^{k}a_ib_{k-i}$

   i   0   1   2   3   4   5   6   7
   A   5   2   3   8   1
   B   4   7   6   2

a0*B  20  35  30  10
a1*B       8  14  12   4
a2*B          12  21  18   6
a3*B              32  56  48  16
a4*B                   4   7   6   2
------------------------------------
   C  20  43  56  75  82  61  22   2

この合成には、様々な意味を持たせることができる。

多項式の掛け合わせ
- 2つの多項式 $f(x)=a_0+a_1x^1+a_2x^2+..., \ \ g(x)=b_0+b_1x^1+b_2x^2+...$
- 積 $f(x)g(x) = c_0+c_1x^1+c_2x^2+...$ になる
信号処理
- 例えば $A$ を「夏休み $i$ 日目に勉強した英単語の数」、 $B$ を「勉強から $i$ 日後に覚えている単語数の割合（0～1）」とすると、 $C$ は「夏休み $i$ 日目に覚えている英単語の数」となる
確率変数の和

計算量については、これをそのまま実践すると $|A||B|$ 個の要素組について計算することになる。 $|A|=|B|=O(N)$ とすると、 $O(N^2)$ かかる。

FFTによる高速化

$A$ をフーリエ変換した結果を $\mathcal{F}[A]$ で表すとすると、以下が成り立つ。 $\circ$ は要素毎の積（アダマール積）を示す。

$\mathcal{F}[A \ast B] = \mathcal{F}[A] \circ \mathcal{F}[B]$

つまり、「 $A,B$ それぞれを（同じ長さで）フーリエ変換」→「要素同士の積を計算」→「逆フーリエ変換で元に戻す」とすると、畳み込みできてしまうらしい。

証明はまだよく理解できていない。他に記事は沢山あるのでそちら参照。中身を追ってみても、フーリエ変換の結果は複素数列となるため具体的に何を指すかが掴みづらい。ここでは、出来るということを知っておくにとどめる。

フーリエ変換・逆変換のコストは、数列長を $N$ として $O(N \log{N})$ 、アダマール積を求めるのは $O(N)$ なので、全体でも $O(N \log{N})$ となる。 $N$ は $C$ の長さ分だけ用意しておけばよいので、 $N=|A|+|B|$ となる。実際には $N$ は2の累乗とするのが計算効率が良いので、適当に後ろをゼロ埋めする。

        
              import numpy as np
 
f = np.array([5, 2, 3, 8, 1])
g = np.array([4, 7, 6, 2])
fft_len = 8  # |f|+|g|-1 以上の2の冪数
 
# f,gをフーリエ変換
Ff = np.fft.rfft(f, fft_len)
print(Ff)
Fg = np.fft.rfft(g, fft_len)
print(Fg)
 
# 各要素について掛け合わせる
Ffg = Ff * Fg
print(Ffg)
 
# 逆変換後、整数に丸める
fg = np.fft.irfft(Ffg, fft_len)
fg = np.rint(fg)
print(fg)

result

F(f)
[19.         +0.j         -0.24264069-10.07106781j
  3.         +6.j          8.24264069 -4.07106781j
 -1.         +0.j        ]

F(g)
[19.         +0.j          7.53553391-12.36396103j
 -2.         -5.j          0.46446609 -0.36396103j
  1.         +0.j        ]
  
F(f)*F(g)
[ 361.         +0.j         -126.34671709-72.89087297j
   24.        -27.j            2.34671709 -4.89087297j
   -1.         +0.j        ]

f*g
[20. 43. 56. 75. 82. 61. 22.  2.]  : 上記の計算結果と一致している

fftとrfftについて

np.fft.fftとnp.fft.rfft - pLog

numpy.fftには、fftとrfftという2種類の関数が存在する。（逆変換もifftとirfftがある）

fftによるフーリエ変換は、 $N$ 要素の数列を入れると $N$ 要素の変換結果が返ってくる。

これ、後ろ半分は負の周波数(?)の情報が格納されているらしい。

フーリエ変換は元の数列が実数であれば、正と負の周波数に同じような数字（複素共役）が現れる性質があるため、通常、後ろ半分の情報は無駄である。

この計算を省略するのがrfftで、前半 $\frac{N+1}{2}$ についてのみ計算することで、半分とまではいかないものの、速くなる。