ピンクノイズを生成するLISP関数

以前PythonやJSFX用のものも作ったのですが、ピンクノイズはいろいろなところで活躍するので、Common Lisp用にも作っておきます。

Voss-McCartney法

DSP generation of Pink (1/f) Noiseには、(1)ホワイトノイズにフィルタをかけてピンクにする方法、(2)Voss-McCartney法、の二つがあると紹介されています。Voss-McCartney法はホワイトノイズを低域に重みをつけて重ね合わせていく方法です。

具体的には、いくつかの乱数生成器を用意して、それらで1, 2, 4, 8, 16, ...サンプルごとに新しい乱数を作らせ、すべてを加算して出力するというものです。たとえば4サンプルごとの生成をする乱数生成器は、4サンプルに一度だけ乱数を更新しますが、それ以外の3サンプルではその前に生成した乱数を返します。より低い頻度（ここでは1/4の標本化周波数）での乱数生成をしていることになりますね。つまり、1, 2, 4, 8, 16, ...サンプルごと生成して加算するということは、オクターブ低いホワイトノイズを次々と作って足し合わせてピンクノイズを作っていくことになります。2ⁿが大きくなればなるほど低い周波数まで伸びるピンクノイズになります。

Common Lispでの実装

とりあえず標本化周波数を決めておきます。

(defparameter +samprate+ 48000)

Voss-McCartney法の本体です。指定した秒数ぶんのピンクノイズを生成します。返り値ですが、多値を使ってピンクノイズ振幅値の列と各サンプルの時刻の列の二つを返しています。どちらも1次元配列（simple-vector）です。

(defun synth-pinknoise (dur &optional (fs +samprate+))
  "Generates pinknoise (wideband noise having -3dB/oct slope) of specified duration (seconds). Uses Voss-McCartney algorithm. http://www.firstpr.com.au/dsp/pink-noise/

Example:
(wavwrite \"pinknoise.wav\" (synth-pinknoise 2.0))"
  (let* ((lowfreq 20)
         (levels (ceiling (log (/ fs lowfreq) 2)))
         (tt (sample-time dur fs))
         (out (make-array (length tt)))
         (x (make-randn-array levels)))
    (dotimes (n (length tt))
      (dotimes (m levels)
        (if (= (1- (expt 2 m)) (mod n (expt 2 (1+ m))))
            (setf (aref x m) (randn)))
        (setf (aref out n) (+ (randn) (vecsum x)))))
    (values (normalize out) tt)))

sample-timeは、指定した秒数の各サンプルの時刻を計算する関数です。横軸に時間・縦軸に振幅をとってグラフを書く時などに使います。標本化周波数が48000 Hzのときに1秒間の時刻列は、#(0 1/48000 2/48000 ... 47999/48000)というようになります。以下のように定義しています。rangeはPythonの、linspaceはMatlabの関数をそれぞれ真似しています。

(defun sample-time (dur &optional (fs +samprate+))
  "Make a vector of time values of each sample in seconds for specified duration."
  (coerce (mapcar #'(lambda (x) (/ x fs))
                  (range (round (* dur fs)))) 'vector))

(defun range (n)
  "Generates a list of integers from 0 to N-1."
  (linspace 0 (1- n) n))

(defun linspace (x1 x2 &optional (n 100))
  "Generates a list of N (100 by default) equally spaced numbers between X1 and X2."
  (let ((x (loop for x from 0 to (1- n)
                 collect (/ x (1- n)))))
    (mapcar #'(lambda (v) (+ v x1))
            (mapcar #'(lambda (v) (* v (- x2 x1)))
                    x))))

randnとmake-randn-arrayは標準正規分布に従う乱数を作る関数で、以下のエントリに書いたものです。

marui.hatenablog.com

vecsumはベクトル内の数の総和を計算します。

(defun vecsum (vec)
  (loop as x across vec
        sum x))

ピンクノイズは音波形として使用したいので、出力時には±1の範囲に収めたいです。これをノーマライズといい、以下のnormalizeで実行しています。±1ビチビチに入れてしまうと（0 dBFSにすると）、オーディオ編集ソフトによっては「レベルオーバーしてるよ」という表示になることがあるので-1 dBFSをデフォルトにしていますが、-3 dBFSとかのほうが安全かもしれません。

(defun normalize (x &optional (dB -1.0))
  "Makes a sampled signal normalized to specified level (in dB)."
  (let ((res (make-array (array-dimensions x)
                         :initial-element 0.0))
        (gain (db->amp dB))
        (tmpmax (loop for n from 0 to (1- (length x))
                      maximize (abs (aref x n)))))
    (dotimes (n (length x))
      (setf (aref res n) (* gain (/ (aref x n) tmpmax))))
    res))

(defun db->amp (dB &optional (p0 1.0))
  (* (expt 10 (/ dB 20)) p0))

できあがったピンクノイズをオーディオ編集ソフトに取り込んで30分にまとめたのが以下の動画です。YouTubeにアップロードするときに聴覚符号化で圧縮されてしまうので、本来の音とは少し違いますが……。

youtu.be