信号の包絡線を計算する (Common Lisp)

信号の包絡線 (signal envelope) は「信号の頂点にシーツをかぶせたときのような概形」を持つ波形です。

MatlabやRで書いた信号包絡線の計算プログラムをCommon Lispに移植しました。信号包絡線の計算方法は振幅の絶対値に対して、(1)低域フィルタをかける、(1')Savitzky-Golayフィルターがよい、(2)ピーク検出をしてスプラインでつなぐ、(2')いやいや秋間補間がよい、(3)ヒルベルト変換を使う、(3')結局それに低域フィルタかけるじゃろ、(4)スペクトル包絡には線形予測符号も使われてるよね、などなどいくつも提案できるのですが、ここではヒルベルト変換を用いて解析信号を計算する方法でいきます。解析信号の絶対値を計算すると信号包絡線になります。

まずは解析信号の実装。フーリエ変換をして、負の周波数成分をゼロにして、逆フーリエ変換します。

(defun hilbert (x)
  "Compute analytic signal of real signal vector X using Hilbert transform."
  (let* ((xx (napa-fft:fft (zeropad-pow2 x)))
         (yy (make-array (length xx) :initial-element #C(0 0)))
         (zz nil)
         (out (make-array (length x))))
    (setf (aref yy 0) (aref xx 0))
    (loop for n from 1 below (/ (length xx) 2)
          do (setf (aref yy n) (* 2 (aref xx n))))
    (setf (aref yy (/ (length xx) 2)) (aref xx (/ (length xx) 2)))
    (loop for n from (1+ (/ (length xx) 2)) below (length xx)
          do (setf (aref yy n) 0))
    (setq zz (napa-fft:ifft yy))
    (loop for n from 0 below (length out)
          do (setf (aref out n) (aref zz n)))
    out))

Matlabのhilbert()と比較してみましょう。

>> hilbert([1 2 -3 4])
ans =
   1.0000 + 1.0000i   2.0000 + 2.0000i  -3.0000 - 1.0000i   4.0000 - 2.0000i

CL-USER> (hilbert #(1 2 -3 4))
#(#C(1.0d0 1.0d0) #C(2.0d0 2.0d0) #C(-3.0d0 -1.0d0) #C(4.0d0 -2.0d0))

同じ結果になっていますね。ヨシ。

次に信号包絡線を計算します。

(defun envelop (x)
  "Calculate amplitude envelop of a signal vector X."
  (let* ((xh (hilbert x))
         (xe (make-array (length xh))))
    (dotimes (n (length xh))
      (setf (aref xe n) (abs (aref xh n))))
    xe))

こちらもMatlabと比較してみましょう。

>> abs(hilbert([1 2 -3 4]))
ans =
    1.4142    2.8284    3.1623    4.4721

CL-USER> (envelope #(1 2 -3 4))
#(1.4142135623730951d0 2.8284271247461903d0 3.1622776601683795d0
  4.47213595499958d0)

本物の音響信号でやってみましょう。スネアドラムを一発叩いた音の録音があるので、信号包絡線を計算してみます。wavreadは以前作ったプログラムです。（→Common Lisp用のwavread - 丸井綜研）

(ql:quickload :clgplot)
(let* ((snd (wavread "snare04-short.wav"))
       (x (first snd))
       (fs (second snd))
       (xe (envelop x))
       (xel (filt (biquad-lpf 100 (sqrt 1/2)) xe))
       (tt (list->vector (mapcar #'(lambda (v) (/ v fs)) (range (length xe))))))
  (clgp:plot x :x-seq tt)
  (clgp:plot xel :x-seq tt))

結果は以下のようになりました。左側が元の信号で右側が包絡線です。横軸は時間（秒）で縦軸は振幅（プラスマイナス1の範囲に正規化）。少しスムーズにするために包絡線にカットオフ周波数100 Hzの低域フィルタをかけてあげています。