音信号処理の代表例として、コンボリューションリバーブ *1のような入力信号にインパルス応答を畳み込む処理が挙げられます。

この畳み込み演算をリアルタイムのブロック処理で実現する*2実用的な方法としては、 Overlap-Add法とOverlap-Save法が有名です。

特に Overlap-Save 法はちゃんと実装したことがなかったため、勉強も兼ねてOverlap-Add法とOverlap-Save法の実装例を紹介します。

Overlap-Add法
Overlap-Save 法
まとめ
参考文献

Overlap-Add法

Overlap-Add法では、ブロックごとの直線畳み込みを効率的に行うため、FFTを用いた巡回畳み込みの高速演算を利用しています。

FFT用の配列長を L+M-1 以上とし、入力信号の計算用配列の x_k のL~M-1サンプルと、インパルス応答の計算用配列 h の M ~ L-1 サンプルをゼロ埋めします。このように計算用配列を用意することで、DFTを用いた巡回畳み込みの計算結果と直線畳み込みの計算結果が同じになります。そして、配列長を調整してFFTを使うことで高速に演算できます。

名前の由来でもあるように、畳み込みの結果は切り出した入力信号よりも M-1 サンプル長くなります。このオーバーラップ(重畳)成分は図のように加算することで最終的な出力信号を得ることができます。

なお、リアルタイム処理時の入力信号のブロック長 L^aが 2のべき乗で設定されていることが多いことと、FFTする信号長 L+M-1 は2のべき乗だと計算効率がよいことを踏まえて、以下の実装では計算用配列の長さを L^a+1として設定しています。

import soundfile as sf
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt


NN = 100
length_fft  = 512
length_x    = 256
length_h    = 64
length_s    = length_x * (NN-1)
length_y    = length_x * (NN-1) + length_h - 1


np.random.seed(2)
sig = np.random.normal(0, 0.5, length_s)
h   = np.random.normal(0, 0.5, length_h)
h   *= np.linspace(1.0,0.0,length_h)**2
h[0] = 1.0
y   = signal.convolve(sig, h)


#=======================================================
# Over-lap add method
y_ola  = np.zeros(length_y)
buf = np.zeros(length_fft)
prev_buf = np.zeros(length_h-1)
x_fft = np.zeros(length_fft, dtype=complex)
y_fft = np.zeros(length_fft, dtype=complex)
h_tmp = np.zeros(length_fft)
h_tmp[0:length_h] = h
h_fft = np.fft.fft(h_tmp)

for k in range(NN):
    if k < NN-1:
        buf[0:length_x] = sig[k*length_x:(k+1)*length_x]
        x_fft[:] = np.fft.fft(buf)
        y_fft[:] = x_fft * h_fft
        buf[:] = np.real(np.fft.ifft(y_fft))
        y_ola[k*length_x:(k+1)*length_x] += buf[0:length_x]

    # 前の buf[length_x:length_x+length_h-1] を加算
    y_ola[k*length_x:k*length_x+length_h-1] += prev_buf
    # 現在のbuf[length_x:length_x+length_h-1] を保存
    prev_buf[:] = buf[length_x:length_x+length_h-1] 
    buf *= 0.0


#plt.plot(y, label="y: scipy convolve")
#plt.plot(y_ola, label="y_ola")
plt.plot(np.abs(y-y_ola), label="np.abs(y-y_ola)")
plt.ylim([-1,1])
plt.legend(loc="lower left")
plt.show()

Overlap-Save 法

Over-lap Add 法よりも若干計算効率が良い方法として Overlap-Save 法があります。

こちらは巡回畳み込みの性質を利用し，直線畳み込みと巡回畳み込みの計算結果が同じになる部分だけを取り出し、計算結果が異なる部分（図の濃い網掛け部）を捨てる仕組みとなっています。

イメージとしては以下のWebサイトの解説がわかりやすいです。

50% Overlap-Save法って、どんなん？

Overlap Add, Overlap Save Visual Explanation

実装の注意点として、Overlap Save 法の入力信号計算用配列 x_k はゼロ埋めではなく、過去の入力信号を計算用配列の前方に入れています（図の薄い網掛け部）。ただし、先頭フレームでは過去の入力信号がないのでゼロ埋めで代用しています。また、最終フレームを計算する際は、現在の入力信号がないので入力信号計算用配列の後ろからLサンプル分をゼロ埋めしています。

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt


NN = 100
length_fft  = 512
length_x    = 256
length_h    = 64
length_s    = length_x * (NN-1)
length_y    = length_x * (NN-1) + length_h - 1


np.random.seed(2)
sig = np.random.normal(0, 0.5, length_s)
h   = np.random.normal(0, 0.5, length_h)
h   *= np.linspace(1.0,0.0,length_h)**2
h[0] = 1.0
y   = signal.convolve(sig, h)


#=======================================================
# Over-lap save method
y_ols  = np.zeros(length_y)
buf = np.zeros(length_fft)
x_fft = np.zeros(length_fft, dtype=complex)
y_fft = np.zeros(length_fft, dtype=complex)
h_tmp = np.zeros(length_fft)
h_tmp[0:length_h] = h
h_fft = np.fft.fft(h_tmp)

# buf[0:length_x]: xでは過去の信号を格納。yでは捨てる
# buf[length_x:length_fft] 現在の信号を格納
for k in range(NN):
    if k < NN-1:
        buf[length_x:length_fft] = sig[k*length_x:(k+1)*length_x]
    else:
        #入力信号長より長くなる成分を計算するために必要
        buf[length_x:length_fft] *= 0.0
    
    # FFTでlength_x*2の信号を循環畳み込み
    x_fft[:] = np.fft.fft(buf)
    y_fft[:] = x_fft * h_fft
        
    # 直線畳み込み結果と一致する後半の信号だけを利用
    buf[length_x:length_fft] = np.real(np.fft.ifft(y_fft))[length_x:length_fft] 
    
    if k < NN-1:
        y_ols[k*length_x:(k+1)*length_x] = buf[length_x:length_fft]
        # 次のフレーム用での過去の入力信号として、現在の入力信号をセット　（ゼロ埋めではない）
        buf[0:length_x] = sig[k*length_x:(k+1)*length_x]
    else:
        #入力信号長より長くなる成分を計算するために必要
        y_ols[k*length_x:k*length_x+length_h-1] = buf[length_x:length_x+length_h-1]
    
plt.clf()
#plt.plot(y, label="y: scipy convolve")
#plt.plot(y_ols, label="y_ols")
plt.plot(np.abs(y-y_ols), label="np.abs(y-y_ols)")
plt.ylim([-1,1])
plt.legend(loc="lower left")
plt.show()