この記事はLisp Advent Calendar 2022の16日目です。Lispで音を自由自在に扱うという野望に向けての次の一歩。
デジタルフィルタ
ちょっと前に高速フーリエ変換を用いて畳み込み演算の実装をしました。2信号の時間領域での畳み込み演算は周波数領域での要素ごとの乗算になるので、フーリエ変換を経由して畳み込み演算を使えば長いFIRフィルタをかけることができます。しかし、この方法はIIRフィルタをかけるのには使えません。また、フーリエ変換を経由するので短いFIRフィルタをかけるのには効率が良くありません。
そこで、今回はデジタルIIRフィルタをかける関数を作ってみました。IIRフィルタはFIRフィルタの上位互換なので、この関数を作ればFIRとIIRのどちらにも使うことができます。MatlabやOctaveのfilter関数を参考にしたかったのですが、いずれの言語でもfilterはビルトイン関数として実装されていて、OctaveのC言語のソースを読んでみても実装のコードが見つかりませんでした（探し方が悪いのか、外部ライブラリを使っているからか……）。そこで、Julia言語のDSP.jlパッケージ内のdspbase.jlを参考にしました。動けばいいので、速度は二の次です。
デジタルフィルタの式は以下のようにあらわします。nは現時刻（サンプル）で、その時刻の入力信号をx(n)、出力信号をy(n)とします。x(n-1)は1サンプル前の入力信号、y(n-2)は2サンプル前の出力信号というわけです。また、B = [b0 b1 b2 ... bM]はフィードフォワード方向のフィルタ係数で、A = [a0 a1 a2 ... aN]はフィードバック方向のフィルタ係数です。うまいことフィルタ係数を設定することで低域フィルタや高域フィルタといったフィルタ特性を作ります。*1
a0 y(n) = b0 x(n) + b1 x(n-1) + b2 x(n-2) + ... + bM x(n-M)
- a1 y(n-1) - a2 y(n-2) - ... - aN y(n-N)
Common Lisp実装
これをLisp関数で計算できるようにしたのが以下のコードです。まずassertを使って入力がベクトルであることと、a0が非零であることを確認します。
係数BとAがともに長さ1のときにはa0 y(n) = b0 x[n]つまりy(n) = b0/a0 x[n]になるので、出力信号は入力信号の定数倍です。一方でフィルタ係数の長さが2以上のときには、過去の入出力信号を参照するのでフィルタの状態を保存しつつ計算する必要があります。その分岐を(if (= 1 clen) ...でやっています。
a0が1でないときには、すべての係数をa0で除算してフィルタ係数の正規化を行います。
そして二重ループの中で、入力信号を1サンプルずつxmに取り出し、当該サンプル時刻の出力valを計算したのちに、フィルタ状態stを更新します。出力信号全体はresに溜まっていきます。
(defun filt (B A x) "Filter a signal X with filter coefficients B (feed-forward) and A (feed-back). B, A, and X are all one dimensional vectors. Length of the output signal is the same as the input signal. Apply zero-padding to the input signal if you need longer output. This function is a naive implementation of Direct Form II Transposed difference equation." (assert (and (typep B 'vector) (typep A 'vector) (typep x 'vector)) nil "All arguments must be vector.") (assert (/= 0 (aref A 0)) nil "The first feed-back coefficient must be non-zero.") (let* ((Blen (length B)) ;feed-forward filter length (Alen (length A)) ;feed-back filter length (clen (max Blen Alen)) ;filter length (xlen (length x)) ;signal length (res (make-array xlen :initial-element 0))) ;filtered signal (if (= 1 clen) ;; simple memory-less scaling (let ((abr (/ (aref B 0) (aref A 0)))) (loop for m from 0 to (1- xlen) do (setf (aref res m) (* (aref x m) abr))) res) ;; with-memory filtering (let* ((st (make-array (1- clen) :initial-element 0)) ;filter state (B2 (make-array clen :initial-element 0)) ;zero-padded coefficients (A2 (make-array clen :initial-element 0))) ;zero-padded coefficients ;; normalize the first feedback filter coefficient to be 1.0 (if (/= (aref A 0) 1) (let ((a0 (aref A 0))) (loop for n from 0 to (1- Alen) do (setf (aref A n) (/ (aref A n) a0))) (loop for n from 0 to (1- Blen) do (setf (aref B n) (/ (aref B n) a0))))) ;; zero-pad the coefficients (loop for m from 0 to (1- Blen) do (setf (aref B2 m) (aref B m))) (loop for m from 0 to (1- Alen) do (setf (aref A2 m) (aref A m))) ;; apply filter (loop for m from 0 to (1- xlen) do (let* ((xm (aref x m)) (val (+ (* (aref B 0) xm) (aref st 0)))) ;calculate output ;; update the filter state (loop for n from 0 to (- (length st) 2) do (setf (aref st n) (+ (* val (- (aref A2 (1+ n)))) (+ (* xm (aref B2 (1+ n))) (aref st (1+ n)))))) (setf (aref st (1- (length st))) (- (* xm (aref B2 (length st))) (* (aref A2 (length st)) val))) ;; set output (setf (aref res m) val))) res))))
使い方
使い方はMatlabのfilter関数と同様で、(filt B A signal)とします。いずれの引数もベクトルで渡してあげます。出力は入力信号と同じ長さになるので、下記のふたつめの例のように必要に応じてゼロ詰めしましょう。
SBCL 2.2.8 CL-USER> (filt #(1 1) #(1) #(3 1 4 1 5 9 2)) #(3 4 5 5 6 14 11) CL-USER> (filt #(1) #(1 -0.5) #(1 0 0 0 0 0 0)) #(1 0.5 0.25 0.125 0.0625 0.03125 0.015625)
Matlab 2022a >> filter([1 1], [1], [3 1 4 1 5 9 2]) ans = 3 4 5 5 6 14 11 >> filter([1], [1 -0.5], [1 0 0 0 0 0 0]) ans = 1.0000 0.5000 0.2500 0.1250 0.0625 0.0312 0.0156
同じ結果になっているので大満足です。やったね。
ちなみに、任意の長さのベクトルを2nの長さにゼロ詰めする関数は以下のように書いてみました。
(defun nextpow2 (n) "Find the exponent of the next smallest power of two larger than n." (ceiling (log n 2))) (defun zeropad (x n) "Return the vector of length of N with the contents of vector X in the earlier part." (assert (>= n (length x))) (let ((xx (make-array n :initial-element 0))) (loop for k from 0 to (1- (length x)) do (setf (aref xx k) (aref x k))) xx)) (defun zeropad-pow2 (x) "Return the vector of length of 2^n with the contents of vector X in the earlier part." (zeropad x (expt 2 (nextpow2 (length x)))))
CL-USER> (zeropad-pow2 #(1 2 3 4 5)) #(1 2 3 4 5 0 0 0)
*1:例えばJulius O. Smith IIIの『Introduction to Digital Filters with Audio Aplications』の「Difference Equation」あたりを参照。 https://ccrma.stanford.edu/~jos/filters/Difference_Equation_I.html
