音声処理まわりの知見を調べた。インパルス応答の畳み込みはやっぱりFFTを使うのが普通らしい。畳み込みをO(N^2)からO(NlogN)に落とすにはそれしかないわな。

なんか波形の区間を切り取ってFFTかけるという話をするとどうも窓関数の話が出てくる。しかし、スペクトラムが見たいのではなく畳み込みをしたいという目的でFFTなら窓関数は要らないのではないか。なんか計算してみた感じやっぱりそんな感じだった。

周波数領域で掛け算を行うと時間領域では畳み込みになるのは常識として知っていたが、詳細な挙動がどうなるのか良く把握していなかったので手を動かして実験した。N=4のとき、信号$s_0,..s_3$とインパルス応答$h_0,…h_3$を畳みこもうとするとこんな感じになる。

$\mathbf{\hat{s}}=\mathbf{Fs}$

$\mathbf{\hat{h}}=\mathbf{Fh}$

$\mathbf{F^{-1}}\left(\begin{array}{c}\hat{s}_0\hat{h}_0 \\ \hat{s}_1\hat{h}_1 \\ \hat{s}_2\hat{h}_2 \\ \hat{s}_3\hat{h}_3\end{array}\right)=\left(\begin{array}{c}s_0h_0 + s_1h_3 + s_2h_2 + s_3h_1 \\ s_0h_1 + s_1h_0 + s_2h_3 + s_3h_2 \\ s_0h_2 + s_1h_1 + s_2h_0 + s_3h_3 \\ s_0h_3 + s_1h_2 + s_2h_1 + s_3h_0\end{array}\right)$

残念ながら、周期的になるのでこのように余計な項が入ってしまう。上記の式で言うと、最右辺の右上の部分は丸っと要らない。$s_0h_0$の後ろの$s_1h_3 + s_2h_2 + s_3h_1$とかを削りたい。

この余計な部分を除くには何をすればいいかと考えれば、hとsの後ろ半分の要素を0にすればいいことが分かる。逆に言えば、使いたいインパルス応答があるならNをそれの2倍に設定して後ろ半分を0で埋めればいい。そうすると上の式から平行四辺形みたいな部分が残る。これで順次切り取ってはFFTを使った畳み込みを行って足していけばよい。

窓の大きさを$N$、音声データ全体の長さを$L$とすると、FFTを行う回数は$O(L/N)$で1回やるごとに$O(N)$回の掛け算と足し算を行う。ので全体の計算量は$O(L/N(N\log N+N))=O(L\log N)$になる。

ところで、インパルス応答が治まるように$N$を定めるのだとすると、最低でもインパルス応答の長さだけの遅延が起きることになる。これは場合によっては許容できない。そこで、インパルス応答の長さを$M$とし、$C$分割して窓の長さを$M/C$とすることを考える。窓の切り取りは$L/(M/C)=CL/M$回あり、1回ごとに大きさ$O(M/C)$のFFTを$O(C)$回、$O(C\times(M/C))=O(M)$回の掛け算と足し算がある。計算量はこうなる。

$O(CL/M(C\times(M/C)\log (M/C)+M))=O(CL\log (M/C))$

$C$倍より少し軽いくらいのようだ。ここで$C=M$とすると普通にやったときと等価になり2乗のオーダーになる。やはり高速な応答を求めるほど負荷は高くなることが分かる。

FFTに使うライブラリも調べたのだが、muFFTというのが良さげだった。MITライセンスだし。このライブラリは以下のベンチマークから知った。

https://github.com/project-gemmi/benchmarking-fft

PostgreSQL用のFFT(?)がやたら高性能なの面白いな…

---

インパルス応答をシミュレーションで出せないかと思い、理論的にシミュレーションする方法を調べた。

https://www.jstage.jst.go.jp/article/jasj/66/9/66_KJ00006579568/_pdf/-char/ja

密度(スカラー場)と粒子速度(ベクトル場)で考えれば良いようだ。

• 質量保存則: $\frac{d\rho}{dt}+\rho_0\frac{d\mathbf{u}}{d\mathbf{x}}=0$
• 運動量保存則: $\rho_0\frac{d\mathbf{u}}{dt}+\frac{dp}{d\mathbf{x}}=0$

シミュレーションの手段としては電磁波とかと同じようにFDTD法とかやろうと思えば行けるらしい。FDTD法によるシミュレーション、未だにちゃんと実装したことないので一度やっときたい。