この記事はLisp Advent Calendar 2022の25日目です。Lispで音プログラミングに向けた、さらなる一歩。
目次
はじめに
先日フィルタをかける関数を作りました（IIRフィルタをかける関数（Lisp Advent Calendar 2022） - 丸井綜研）。となると、今度はフィルタが欲しくなります。Hann窓のような窓関数を作って時間信号に畳み込んでもいいですが、せっかくIIRフィルタが使えるようになったので、IIRフィルタの係数を計算していきます。
IIRフィルタと言っても、今回はシンプルに双二次フィルタ（Digital biquad filter - Wikipedia）にします。しかし自分で設計するのは大変なので、Robert Bristow-Johnson氏による「Cookbook Formulae for Audio Equalizer Biquad Filter Coefficients」という便利なフィルター設計事例集をそのまま使ってしまいましょう。（ただしフィルタを適用する周波数がナイキスト周波数に近いときには、理想的なアナログフィルタの挙動とは異なってしまいます。例えば、Orfanidis “Digital Parametric Equalizer Design with Prescribed Nyquist-Frequency Gain” (J. Audio Engineering Soc., 1997)などがこれを解決する方法を提示してくれています。）
RBJ Cookbookによる低域フィルタ
双二次フィルタは次のような伝達関数で、b0〜b2とa0〜a2の6つの係数を決めてあげればフィルタ形状が決まります。
b0 + b1*z^-1 + b2*z^-2
H(z) = ------------------------
a0 + a1*z^-1 + a2*z^-2
設計事例集には低域フィルタはこんな形で掲載されています。
b0 = (1 - cos(w0))/2
b1 = 1 - cos(w0)
b2 = (1 - cos(w0))/2
a0 = 1 + alpha
a1 = -2*cos(w0)
a2 = 1 - alpha
ここでw0とalphaがパラメータになっていますが、設計事例集の前半にそれぞれの計算式が載っています。
w0 = 2*pi*f0/Fs
alpha = sin(w0)/(2*Q) (case: Q)
= sin(w0)*sinh( ln(2)/2 * BW * w0/sin(w0) ) (case: BW)
= sin(w0)/2 * sqrt( (A + 1/A)*(1/S - 1) + 2 ) (case: S)
f0はカットオフ周波数、Fsは標本化周波数です。alphaを決める方法が3種類提示されていますが、今回はQ値を使いたいのでsin(w0)/(2*Q)を採用します。
それをCommon Lispの関数にするのは楽々。
(defun biquad-lpf (fc q fs) (let* ((omega (/ (* 2 pi fc) fs)) (alpha (/ (sin omega) (* 2 q))) (b0 (/ (- 1 (cos omega)) 2)) (b1 (- 1 (cos omega))) (b2 (/ (- 1 (cos omega)) 2)) (a0 (+ 1 alpha)) (a1 (* -2 (cos omega))) (a2 (- 1 alpha))) (list `#(,b0 ,b1 ,b2) `#(,a0 ,a1 ,a2))))
カットオフ周波数を1000 Hz、Q値を1/√2、標本化周波数を44100 Hzとして計算してみましょう。
CL-USER> (setq fltr (biquad-lpf 1000 (sqrt 1/2) 44100)) ;;=> (#(0.005066263610029209d0 0.010132527220058418d0 0.005066263610029209d0) ;; #(1.1004051468361575d0 -1.9797349455598832d0 0.8995948531638425d0))
フィルタを使ってみる
双二次フィルタの係数の計算ができました。(first fltr)と(second fltr)でそれぞれ分子と分母の係数を取ってくることができます。テスト信号として、44100サンプルのゼロのうち、8番目だけを1にした信号を作り、低域フィルタをかけてみます。つまり、フィルタのインパルス応答を見ようというわけです。
(setq fltrd-sig (let* ((sig (make-array 44100)) (fc 1000) (q (sqrt 1/2)) (fs 44100)) (setf (aref sig 8) 1) (setq fcoef (biquad-lpf fc q fs)) (filt (first fcoef) (second fcoef) sig)))
masatoiさんによるGnuPlotライブラリがmacOSでもそのまま使えたので、グラフも作れるようになりました（Common LispからGnuplotでグラフ表示するライブラリを作った： clgplot - masatoi’s blog）。clgplotはRoswellを使ってros install masatoi/clgplotとしてインストールしてあります。
(ql:quickload :clgplot) (clgplot:plots (list sig (filt (first fcoef) (second fcoef) sig)) :style '(impulse line) :title-list '("original impulse" "filtered signal") :x-label "Time (samples)" :y-label "Amplitude" :x-range '(0 127) :output-format :png :output "clgp-output.png")
紫色の線が入力（インパルス）で、緑色の線がフィルタからの出力です。波形が滑らかになっていて高域がなくなっていそうなのがわかります。
以前作成したwavreadとwavwriteを使って、音ファイルにもフィルタをかけてみましょう。
(multiple-value-bind (x fs) (wavread "finn.wav") (setq fcoef (biquad-lpf 1000 (sqrt 1/2) fs)) (wavwrite "biquad-output.wav" (filt (first fcoef) (second fcoef) x) fs))
できあがったのがこちらです。前半が原音、後半が低域フィルタからの出力音です。後半は高い音が聞こえにくくなり、こもったような音質になりましたね。
とはいえ、(filt (first fcoef) (second fcoef) x)と呼び出すのは面倒なので、filt関数の引数を(filt <filt-coef> <signal>)みたいにした方がいい気がしてきました。そうすれば(filt (biquad-lpf 1000 2 48000) sig)みたいに使えますし。（そろそろCLOSのお勉強しようかな……）
他の形状のフィルタ
その後cosを呼び出す回数を減らしたりコメントを入れたりしたので、以下にCookbookで紹介されていたすべてのフィルタ係数を計算する関数を置いておきます。
(defun biquad-lpf (fc q fs) "Biquad Low-Pass Filter. Returns a list of two filter coefficient vectors B and A. The implementation follows ``Cookbook formulae for audio equalizer biquad filter coefficients'' by Robert Bristow-Johnson. Example: (biquad-lpf 1000 (sqrt 1/2) 44100) => (#(0.005066263610029209d0 0.010132527220058418d0 0.005066263610029209d0) #(1.1004051468361575d0 -1.9797349455598832d0 0.8995948531638425d0))" (assert (> fc 0) nil "Cutoff frequency must be positive.") (assert (> fs 0) nil "Sampling frequency must be positive.") (assert (> (/ fs 2) fc) nil "Cutoff frequency must be less than the Nyquist frequency.") (assert (> q 0) nil "Q-value must be positive.") (let* ((omega (/ (* 2 pi fc) fs)) (cosw (cos omega)) (sinw (sin omega)) (alpha (/ sinw (* 2 q))) (b0 (/ (- 1 cosw) 2)) (b1 (- 1 cosw)) (b2 (/ (- 1 cosw) 2)) (a0 (+ 1 alpha)) (a1 (* -2 cosw)) (a2 (- 1 alpha))) (list `#(,b0 ,b1 ,b2) `#(,a0 ,a1 ,a2)))) (defun biquad-hpf (fc q fs) "Biquad High-Pass Filter. Returns a list of two filter coefficient vectors B and A. The implementation follows ``Cookbook formulae for audio equalizer biquad filter coefficients'' by Robert Bristow-Johnson. Example: (biquad-hpf 1000 (sqrt 1/2) 44100) => (#(0.9949337363899708d0 -1.9898674727799417d0 0.9949337363899708d0) #(1.1004051468361575d0 -1.9797349455598832d0 0.8995948531638425d0))" (assert (> fc 0) nil "Cutoff frequency must be positive.") (assert (> fs 0) nil "Sampling frequency must be positive.") (assert (> (/ fs 2) fc) nil "Cutoff frequency must be less than the Nyquist frequency.") (assert (> q 0) nil "Q-value must be positive.") (let* ((omega (/ (* 2 pi fc) fs)) (cosw (cos omega)) (sinw (sin omega)) (alpha (/ sinw (* 2 q))) (b0 (/ (+ 1 cosw) 2)) (b1 (- -1 cosw)) (b2 (/ (+ 1 cosw) 2)) (a0 (+ 1 alpha)) (a1 (* -2 cosw)) (a2 (- 1 alpha))) (list `#(,b0 ,b1 ,b2) `#(,a0 ,a1 ,a2)))) (defun biquad-bpf (fc q fs) "Biquad Band-Pass Filter (constant 0 dB peak gain). Returns a list of two filter coefficient vectors B and A. The implementation follows ``Cookbook formulae for audio equalizer biquad filter coefficients'' by Robert Bristow-Johnson. Example: (biquad-bpf 1000 (sqrt 1/2) 44100) => (#(0.10040514683615755d0 0 -0.10040514683615755d0) #(1.1004051468361575d0 -1.9797349455598832d0 0.8995948531638425d0))" (assert (> fc 0) nil "Cutoff frequency must be positive.") (assert (> fs 0) nil "Sampling frequency must be positive.") (assert (> (/ fs 2) fc) nil "Cutoff frequency must be less than the Nyquist frequency.") (assert (> q 0) nil "Q-value must be positive.") (let* ((omega (/ (* 2 pi fc) fs)) (cosw (cos omega)) (sinw (sin omega)) (alpha (/ sinw (* 2 q))) (b0 alpha) (b1 0) (b2 (- alpha)) (a0 (+ 1 alpha)) (a1 (* -2 cosw)) (a2 (- 1 alpha))) (list `#(,b0 ,b1 ,b2) `#(,a0 ,a1 ,a2)))) (defun biquad-notch (fc q fs) "Biquad Notch Filter. Returns a list of two filter coefficient vectors B and A. The implementation follows ``Cookbook formulae for audio equalizer biquad filter coefficients'' by Robert Bristow-Johnson. Example: (biquad-notch 1000 (sqrt 1/2) 44100) => (#(1 -1.9797349455598832d0 1) #(1.1004051468361575d0 -1.9797349455598832d0 0.8995948531638425d0))" (assert (> fc 0) nil "Cutoff frequency must be positive.") (assert (> fs 0) nil "Sampling frequency must be positive.") (assert (> (/ fs 2) fc) nil "Cutoff frequency must be less than the Nyquist frequency.") (assert (> q 0) nil "Q-value must be positive.") (let* ((omega (/ (* 2 pi fc) fs)) (cosw (cos omega)) (sinw (sin omega)) (alpha (/ sinw (* 2 q))) (b0 1) (b1 (* -2 cosw)) (b2 1) (a0 (+ 1 alpha)) (a1 (* -2 cosw)) (a2 (- 1 alpha))) (list `#(,b0 ,b1 ,b2) `#(,a0 ,a1 ,a2)))) (defun biquad-apf (fc q fs) "Biquad All-pass Filter. Returns a list of two filter coefficient vectors B and A. The implementation follows ``Cookbook formulae for audio equalizer biquad filter coefficients'' by Robert Bristow-Johnson. Example: (biquad-apf 1000 (sqrt 1/2) 44100) => (#(0.8995948531638425d0 -1.9797349455598832d0 1.1004051468361575d0) #(1.1004051468361575d0 -1.9797349455598832d0 0.8995948531638425d0))" (assert (> fc 0) nil "Cutoff frequency must be positive.") (assert (> fs 0) nil "Sampling frequency must be positive.") (assert (> (/ fs 2) fc) nil "Cutoff frequency must be less than the Nyquist frequency.") (assert (> q 0) nil "Q-value must be positive.") (let* ((omega (/ (* 2 pi fc) fs)) (cosw (cos omega)) (sinw (sin omega)) (alpha (/ sinw (* 2 q))) (b0 (- 1 alpha)) (b1 (* -2 cosw)) (b2 (+ 1 alpha)) (a0 (+ 1 alpha)) (a1 (* -2 cosw)) (a2 (- 1 alpha))) (list `#(,b0 ,b1 ,b2) `#(,a0 ,a1 ,a2)))) (defun biquad-peak (fc gain q fs) "Biquad Peaking EQ Filter. Returns a list of two filter coefficient vectors B and A. The implementation follows ``Cookbook formulae for audio equalizer biquad filter coefficients'' by Robert Bristow-Johnson. Example: (biquad-peak 1000 12 (sqrt 1/2) 44100) => (#(1.3997200784175376d0 -1.9797349455598832d0 0.6002799215824622d0) #(1.1004051491462712d0 -1.9797349455598832d0 0.8995948508537289d0))" (assert (> fc 0) nil "Cutoff frequency must be positive.") (assert (> fs 0) nil "Sampling frequency must be positive.") (assert (> (/ fs 2) fc) nil "Cutoff frequency must be less than the Nyquist frequency.") (assert (> q 0) nil "Q-value must be positive.") (let* ((omega (/ (* 2 pi fc) fs)) (cosw (cos omega)) (sinw (sin omega)) (A (expt 10 (/ gain 40))) (q (/ q A)) ; adjustment suggested by RBJ (alpha (/ sinw (* 2 q))) (b0 (+ 1 (* alpha A))) (b1 (* -2 cosw)) (b2 (- 1 (* alpha A))) (a0 (+ 1 (/ alpha A))) (a1 (* -2 cosw)) (a2 (- 1 (/ alpha A)))) (list `#(,b0 ,b1 ,b2) `#(,a0 ,a1 ,a2)))) (defun biquad-lowshelf (fc gain q fs) "Biquad Low-shelf Filter. Returns a list of two filter coefficient vectors B and A. The implementation follows ``Cookbook formulae for audio equalizer biquad filter coefficients'' by Robert Bristow-Johnson. Example: (biquad-lowshelf 1000 12 (sqrt 1/2) 44100) => (#(4.576605860823849d0 -7.859937860688971d0 3.4446852566817814d0) #(4.264091974420542d0 -7.920349671951922d0 3.6967878039958317d0))" (assert (> fc 0) nil "Cutoff frequency must be positive.") (assert (> fs 0) nil "Sampling frequency must be positive.") (assert (> (/ fs 2) fc) nil "Cutoff frequency must be less than the Nyquist frequency.") (assert (> q 0) nil "Q-value must be positive.") (let* ((omega (/ (* 2 pi fc) fs)) (cosw (cos omega)) (sinw (sin omega)) (A (expt 10 (/ gain 40))) (alpha (/ sinw (* 2 q))) (b0 (* A (+ (+ A 1) (* -1 (- A 1) cosw) (* 2 (sqrt A) alpha)))) (b1 (* 2 A (- (- A 1) (* (+ A 1) cosw)))) (b2 (* A (- (+ A 1) (* (- A 1) cosw) (* 2 (sqrt A) alpha)))) (a0 (+ (+ A 1) (* (- A 1) cosw) (* 2 (sqrt A) alpha))) (a1 (* -2 (+ (- A 1) (* (+ A 1) cosw)))) (a2 (+ (+ A 1) (* (- A 1) cosw) (* -2 (sqrt A) alpha)))) (list `#(,b0 ,b1 ,b2) `#(,a0 ,a1 ,a2)))) (defun biquad-highshelf (fc gain q fs) "Biquad High-shelf Filter. Returns a list of two filter coefficient vectors B and A. The implementation follows ``Cookbook formulae for audio equalizer biquad filter coefficients'' by Robert Bristow-Johnson. Example: (biquad-lowshelf 1000 12 (sqrt 1/2) 44100) => (#(4.576605860823849d0 -7.859937860688971d0 3.4446852566817814d0) #(4.264091974420542d0 -7.920349671951922d0 3.6967878039958317d0))" (assert (> fc 0) nil "Cutoff frequency must be positive.") (assert (> fs 0) nil "Sampling frequency must be positive.") (assert (> (/ fs 2) fc) nil "Cutoff frequency must be less than the Nyquist frequency.") (assert (> q 0) nil "Q-value must be positive.") (let* ((omega (/ (* 2 pi fc) fs)) (cosw (cos omega)) (sinw (sin omega)) (A (expt 10 (/ gain 40))) (alpha (/ sinw (* 2 q))) (b0 (* A (+ (+ A 1) (* (- A 1) cosw) (* 2 (sqrt A) alpha)))) (b1 (* -2 A (+ (- A 1) (* (+ A 1) cosw)))) (b2 (* A (+ (+ A 1) (* (- A 1) cosw) (* -2 (sqrt A) alpha)))) (a0 (+ (+ A 1) (* -1 (- A 1) cosw) (* 2 (sqrt A) alpha))) (a1 (* 2 (- (- A 1) (* (+ A 1) cosw)))) (a2 (- (+ A 1) (* (- A 1) cosw) (* 2 (sqrt A) alpha)))) (list `#(,b0 ,b1 ,b2) `#(,a0 ,a1 ,a2))))
