先日Common LispでFFTが使えるようになったので、それを使って畳み込みをします。畳み込みには環状畳み込み（円状畳み込みや循環畳み込みとも呼びます）と線形畳み込みの2種類がありますので、両方とも作ってみようと思います。
ちょっとした道具
実際の畳み込みのプログラムに入る前に、Matlabを参考にしてnextpow2を作ります。これは与えられた値より大きい最小の2のべき乗の指数を計算するものです。たとえば5の次の2のべき乗は8=23なので3、12の次の2のべき乗は16=24なので4が返ります。（ceilingは多値を返すので実際には余りも出てきますが、ここでは1つ目の返り値だけを使います。）
(defun nextpow2 (n) "Find the exponent of the next smallest power of two larger than n." (ceiling (log n 2))) (nextpow2 5) ;=> 3 (nextpow2 12) ;=> 4
これを使うと、たとえば任意のベクトルを渡すと、2のべき乗の長さの先頭部分にそのベクトル内容を入れて残りはゼロ詰めしたものを返すような関数が作れます。
(defun zeropad-pow2 (vec) "Return the array of length of 2^n with the contents of <vec> in the earlier part." (let ((xx (make-array (expt 2 (nextpow2 (length vec)))))) (loop for k from 0 to (1- (length vec)) do (setf (aref xx k) (aref vec k))) xx)) (zeropad-pow2 #(1 2 3)) ;=> #(1 2 3 0) (zeropad-pow2 #(1 2 3 4 5)) ;=> #(1 2 3 4 5 0 0 0)
環状畳み込み（circular convolution）
さて、ここからが本題です。ある程度の長さのあるベクトルどうしを畳み込みする際には、高速フーリエ変換を経由すると高速に実行できます。すなわち、vec1とvec2の畳み込みをするときに、それぞれをフーリエ変換してy1とy2を計算してから、y1とy2の要素ごとの積をとり、その結果を逆フーリエ変換するという手順です。なおFFTの制約もあったりして、vec1とvec2は同じ長さで、かつ2のべき乗の長さになっている必要があります。
(defun fftconv-circular (vec1 vec2) "Circular convolution using FFT. The two vectors must be of equal length as well as 2^n in length." (coerce (mapcar #'realpart (coerce (bordeaux-fft:sifft (let ((y1 (bordeaux-fft:sfft vec1)) (y2 (bordeaux-fft:sfft vec2))) (loop for e1 across y1 and e2 across y2 collect (* e1 e2)))) 'list)) 'simple-vector)) (fftconv-circular #(1 2 3 4) #(1 1 0 0)) ;=> #(5.0d0 3.0d0 5.0d0 7.0d0)
上記のコードでは、逆フーリエ変換後に実部だけを取ってくるために一旦ベクトルをリストにしてmapcarをかけています。なんとなく効率よくない気がするので、ベクトルのままループで実部を取るものも作ってみました。
(defun fftconv-circular (vec1 vec2) "Circular convolution using FFT. The two vectors must be of equal length as well as 2^n in length." (let ((z (bordeaux-fft:sifft (let ((y1 (bordeaux-fft:sfft vec1)) (y2 (bordeaux-fft:sfft vec2))) (loop for e1 across y1 and e2 across y2 collect (* e1 e2))))) (res (make-array (length vec1)))) (loop for k from 0 to (1- (length vec1)) do (setf (aref res k) (realpart (aref z k)))) res))
どちらのほうが実行速度が速いのか、以下のように乱数列どうしの畳み込みで調べてみました。Intel i5 10400で実行した結果、どちらも14～15秒くらいでした。リストに変換してmapcarするよりも、ループを使って実部をとるほうが3～4%ほど速いことが分かりましたので、今後はそのやりかたで行こうと思います。
（比較のためにMatlab 2022aのSignal Processing Toolboxに入っているcconvを使ってtic; for n=1:100; cconv(rand(2^18,1), rand(2^18,1)); end; tocとしたところ10秒くらいでした。Julia 1.8.1でもDSP.jl 0.7.7を使って@time for n=1:100; conv(rand(2^18), rand(2^18)); endをしたら2.5秒でした。やはりJulia速い……）
(defun make-random-array (N) (let ((res (make-array N))) (loop for k from 0 to (1- (array-total-size res)) do (setf (row-major-aref res k) (random 1.0d0))) res)) (time (dotimes (k 100) (let ((x (fftconv-circular (make-random-array (expt 2 18)) (make-random-array (expt 2 18))))))))
線形畳み込み（linear convolution）
線形畳み込みでは、vec1の長さ + vec2の長さ - 1の長さの結果が得られます。たとえば[1 2 3 4]と[1 1]の畳み込み結果は[1 3 5 7 4]になります。bordeaux-fftでは2のべき乗の長さのベクトルしかFFTがかけられないのと、フーリエ領域で要素ごとの積を計算するときには2つのベクトルの長さが一致していなければいけないことから、vec1の長さ + vec2の長さ - 1を超える最小の2のべき乗のベクトルを作り、ゼロ詰めを行います。[1 2 3 4]と[1 1]の例では畳み込み後の長さが4+2-1=5になるので、xx1とxx2はそれぞれ[1 2 3 4 0 0 0 0]と[1 1 0 0 0 0 0 0]というベクトルになります。これらを環状畳み込みに渡してあげると[1 3 5 7 4 0 0 0]という結果がもどってくるので、先頭の5要素を抜き出してあげるというわけです。
(defun fftconv-linear (vec1 vec2) "Linear convolution using FFT." (let ((xx1 (make-array (expt 2 (nextpow2 (+ (length vec1) (length vec2) -1))))) (xx2 (make-array (expt 2 (nextpow2 (+ (length vec1) (length vec2) -1))))) (res (make-array (+ (length vec1) (length vec2) -1)))) (loop for k from 0 to (1- (length vec1)) do (setf (aref xx1 k) (aref vec1 k))) (loop for k from 0 to (1- (length vec2)) do (setf (aref xx2 k) (aref vec2 k))) (setq foo (fftconv-circular xx1 xx2)) (loop for k from 0 to (1- (length res)) do (setf (aref res k) (aref foo k))) res)) (fftconv-linear #(1 2 3 4) #(1 1)) ;=> #(1.0d0 3.0d0 5.0d0 7.0d0 4.0d0)
