はじめに

FIRフィルタは畳み込み積分で実現するフィルタです。そのまま計算する場合、処理は非常にシンプルですが、計算量も非常に多く、フィルタカーネル長が長くなればなるほど非現実的な処理時間になってしまいます。この問題に対して、DFTを応用して高速化する手法が知られています。

DFTによる畳み込みの定理

以前紹介した通り、DFTの性質として、原信号領域の畳み込み積分は、変換領域における積になる畳み込みの定理が存在します。

畳み込みを行いたい信号列をそれぞれDFTでスペクトラムへ変換し、スペクトラム同士で積を計算して、その結果を原信号領域へIDFTで戻すと、原信号領域で巡回畳み込み積分を行なった結果と等しくなります。

巡回畳み込み積分とは、畳み込み積分範囲を信号長の長い方に合わせてゼロ値で拡張した後、どちらでも良いのでどちらかの信号をローテーションさせながら計算した畳み込み積分です。

信号長は有限のため、畳み込み積分を計算する際のはみ出した部分については、巡回したもの、もしくは、対象信号を周期信号と仮定して周期的に同じ信号が続くものとして計算します。

DFTで畳み込み積分を計算する

長さ $L$ の信号 $A$ に対して、長さ $K$ のフィルタカーネル $B$ で畳み込み積分を計算するとします。
巡回畳み込み積分では、フィルタカーネルが対象信号からはみ出した部分は巡回したものとして、いわば嘘の信号に対して畳み込み積分を行います。この嘘の信号に対する畳み込み積分の結果は、畳み込み積分の結果 $L$ 長の端点、長さ $K-1$ 分だけ発生します。DFTによる畳み込み積分は信号の過去方向へ遡った結果のため、巡回部分は計算結果の先頭部分です。そのため、巡回部分は畳み込み積分の結果の先頭から $K - 1$ 部分となります。

計算の概略は以下の通りです。

フィルタカーネルの末尾にゼロ値を挿入して長さ $L$ の信号 $B’$ にする
信号 $A$ と $B’$ をそれぞれDFTでスペクトラムへ変換する
$DFT(A) \times DFT(B’)$の結果をIDFTで原信号領域へ戻す
巡回畳み込み部分である先頭 $K-1$ を取り除く

少々面倒なのは、DFTをFFTで実装する場合、信号長が2の冪乗でなければならない点です。信号 $A$ の長さ $L$ が2の冪乗ではない場合、要は信号全体の末尾部分がちょうど2の冪乗長ではない場合ですが、これは通常の畳み込み積分と同じ様に拡張するか、計算結果を $K-1$ 分だけ捨ててしまいます。

また、信号長が非常に長い場合に、メモリ使用量の都合でそのままではFFTで計算できない可能性もあるため、短い窓長に分割してDFTを施し、結果を連結する、という操作も必要になるかもしれません。

Pythonによる実装実験

__comvolve関数がDFTによる畳み込み積分の実装です。この計算結果とNumpy.convolveによる計算結果を比較します。

入力信号はノコギリ波で、カーネルフィルタはなんとなくで決めた移動平均フィルタのカーネルです。

import scipy.fftpack
import numpy

def __convolve(signal, kernel):
    fft = scipy.fftpack.fft
    ifft = scipy.fftpack.ifft
    kernel_ = numpy.r_[kernel, numpy.zeros(len(signal) - len(kernel))]
    filtered = numpy.real(ifft(fft(signal) * fft(kernel_)))
    return filtered[(len(kernel) - 1):]

def __make_saw_wave(length, freq):
    a = length / freq
    t = numpy.linspace(0, length - 1, length)
    saw_wave = 2 * (t / a - numpy.floor(t / a + 1 / 2))
    return saw_wave

import matplotlib.pyplot as plt

if __name__ == "__main__":
    saw_wave = __make_saw_wave(256, 8)
    ones = numpy.ones(8)
    filtered_a = __convolve(saw_wave, ones)
    filtered_b = numpy.convolve(saw_wave, ones, mode="valid")

    plt.title("Result of convolution by DFT")
    plt.plot(saw_wave, label="SAW wave")
    plt.plot(filtered_a, label="by DFT", linewidth=4)
    plt.plot(filtered_b, label="by Numpy", linewidth=4, linestyle="dashed")
    plt.legend()
    plt.show()
    plt.close()
    exit()

結果は以下の通りです。信号端の境界処理について厳密に検証できていませんが、非巡回部分については期待通りの結果になっていそうです。

おわりに

DFTの畳み込みの定理は、窓関数の周波数応答を理解するための理論的なツールとしてだけではなく、畳み込み積分の高速化という工学的な実用面での応用にも役立つ重要な概念である事がわかりました。

この手法により、より対象範囲を広げた自己相関関数の計算も高速化して行うことができますし、IIRリバーブの高速化にも応用可能です。