読者です 読者をやめる 読者になる 読者になる

fftfiltを使うときに注意すること

Matlab 2011bで、単純な低域フィルタ [0.5 1 0.5] とインパルス的な [5 1 1 1 1] をfftfiltを使って畳み込みしてみたところ、fftfiltをそのまま使った場合、尻尾が切られてしまいます。

>> fftfilt([0.5 1 0.5], [5 1 1 1 1])
ans =
    2.5000    5.5000    4.0000    2.0000    2.0000

入力信号に対してフィルタをかけたとき、出力信号の長さが入力信号と同じほうが良い場合もあるので、この「尻尾が切られる」挙動はダメというわけではありませんが、これが問題になることもあります。

では「尻尾が切られない」とはどういうことでしょうか。以下のfftconvというのは自作の畳み込み演算プログラムで、線形畳み込みと環状畳み込みが選択できるようになってます*1

まず線形畳み込みだと、出力信号の尻尾は畳み込み演算が行われる範囲全体となります。

>> fftconv([0.5 1 0.5], [5 1 1 1 1], 'linear')'
ans =
    2.5000    5.5000    4.0000    2.0000    2.0000    1.5000    0.5000

出力信号長は「入力信号1の長さ+入力信号2の長さ−1」となっています。これが教科書的な畳み込みです。

一方、環状畳み込みが欲しいときもあります。入力信号の長さからはみ出た尻尾部分が前半に加算されるというもので、fftconvで計算してみると以下のようになります。

>> fftconv([0.5 1 0.5], [5 1 1 1 1], 'circular')'
ans =
    4.0000    6.0000    4.0000    2.0000    2.0000

線形畳み込みのときにはみ出た[1.5 0.5]が、出力信号[2.5 5.5 4.0 2.0 2.0]の前半に回り込んで加算されることで[2.5+1.5 5.5+0.5 4.0 2.0 2.0]となるのです。*2

線形にしても環状にしても、fftfiltを「フィルタをかける」のではなく「畳み込み」に使用するのであれば、尻尾を捨ててしまう挙動は望ましくありません。簡単な対策としては、fftfiltを使う前に入力信号のゼロ詰めをしてやることで、線形畳み込みができます。

>> fftfilt([0.5 1 0.5], [5 1 1 1 1 0 0])
ans =
    2.5000    5.5000    4.0000    2.0000    2.0000    1.5000    0.5000

もしくは、自分用にFFTを使った畳み込みの関数を定義してしまってもいいでしょう。ゼロ詰めをしてからFFTをして、積を取り、IFFTしています。

function fftconv(x, y)
  X = fft([x ; zeros(length(y)-1, 1)]);
  Y = fft([y ; zeros(length(x)-1, 1)]);
  z = real(ifft(X .* Y));

*1:ちなみにOctaveには線形畳み込みを行うfftconvという関数が用意されています。

*2:線形・環状(もしくは直線状・円状)の使い分けは、たとえば[http://tosa.mri.co.jp/sounddb/tsp/index.htm:title=TSPを用いたインパルス応答の測定]などで考慮する必要が出てきます。