フィルタの可視化（freqz相当の関数を作る）

信号にフィルタをかけられるようになり、そして、レシピ通りのフィルタを作れるようになりました。

次に、フィルタの周波数特性（f特）を見られるようにしたいと思います。何年か前に、Matlabのfreqzに似たものをJuliaに作ったことがあります。その時は2次のフィルタまでしか対応していなかったので、今回Common Lisp用に作るにあたって何次のフィルタでも使えるものにしたいと思います。とは言っても、やっていることはそんなに大きく違いはありません。周波数ごとに伝達関数の値を計算するだけです。（linspaceは最小値〜最大値を指定した数に均等割した値の入ったリストを返します）

(defun freqz (fltr &optional (fs 48000) (len 1024))
  (let* ((b (first fltr))
         (a (second fltr))
         (f (linspace 0 (/ fs 2) len))
         (z (mapcar #'(lambda (x) (/ (* 2 pi x) fs)) f))
         (s (mapcar #'(lambda (x) (exp (* #C(0 -1) x))) z))
         (H (make-array len :initial-element #C(0 0))))
    (loop for i from 0 to (1- len) do
      (setf (aref H i)
            (/ (loop for ib from 0 to (1- (length b))
                     sum (* (aref b ib) (expt (nth i s) ib)))
               (loop for ia from 0 to (1- (length a))
                     sum (* (aref a ia) (expt (nth i s) ia))))))
    (values H (coerce f 'vector))))

(freqz '(#(b0 b1 b2...) #(a0 a1 a2...)))とすれば、多値で周波数ごとの応答と周波数の値が返ってきます。ここではデフォルトで標本化周波数48000 Hz、計算する点数を1024点としてありますが、次のように指定すれば、任意の標本化周波数と点数の計算ができます。（これは遅延素子ひとつだけのときの櫛形フィルタで、DCでは振幅2倍、ナイキストで振幅ゼロになっています。）

(freqz '(#(1 1) #(1)) 100.0 11)
;;=> #(#C(2.0 0.0) #C(1.951 -0.309) #C(1.809 -0.588) #C(1.588 -0.809)
;;     #C(1.309 -0.951 #C(1.0 -1.0) #C(0.691 -0.951) #C(0.412 -0.809)
;;     #C(0.191 -0.588) #C(0.049 -0.309) #C(0.0 0.0))
;;   #(0.0 5.0 10.0 15.0 20.0 25.0 30.0 35.0 40.0 45.0 50.0)

先日作った双二次フィルタから、ピークフィルタ（1000 Hz、+12 dB、Q=2）の形を見てみましょう。今回もmasatoiさんのGnuPlotライブラリを使わせてもらっています。freqzは周波数特性を複素数で返しますので、振幅特性と位相特性を計算してグラフに保存してみます。（array-mapは配列用のmap関数、amp->dbは振幅比をデシベル値に変換する関数です）

(ql:quickload :clgplot)

(setq fltr (biquad-peak 1000 12 2 48000))
(multiple-value-bind (amp freq) (freqz fltr)
  (clgp:plot (array-map #'amp->db (array-map #'abs amp))
             :x-seq freq :x-logscale t :key nil
             :x-label "Frequency (Hz)"
             :y-label "Amplitude (dB)"
             :output "freqz-amplitude.png" :output-format :png)
  (clgp:plot (array-map #'phase amp)
             :x-seq freq :x-logscale t :key nil
             :x-label "Frequency (Hz)"
             :y-label "Phase (rad)"
             :output "freqz-phase.png" :output-format :png))