ピンクノイズを作るPython関数

Voss法でピンクノイズを作るプログラムをMatlabからPythonに移植しました。ベクトルではなくループで書いているところなんて、Python的・numpy的にも汚い書き方だと思うので、適当に修正して使って下さい。ピンクノイズについては「DSP Generation of Pink Noise」が詳しく、下記プログラムもこのページを全面的に参考にさせていただいています。

Audacityで10秒間のピンクノイズのスペクトルを見てみたところ、周波数特性もほぼ-10 dB/decade（≈ -3 dB/octave）になってるみたいですし、まぁまぁなんじゃないでしょうか。

# 必要なライブラリ
import math
import numpy as np

# ピンクノイズの関数定義
def synth_pinknoise(duration=1.0, samprate=44100):
    """
    Generate pinknoise using Voss algorithm
    http://www.firstpr.com.au/dsp/pink-noise/
    """
    f_low = 10   # lowest frequency to keep pink (Hz)
    levels = math.ceil(np.log2(samprate/f_low))
    t = np.arange(start=0, stop=duration, step=1/samprate)
    out = np.zeros(len(t))
    x = np.random.normal(size=levels)
    for n in range(math.ceil(samprate * duration)):
        for m in range(levels):
            if n % 2**(m+1) == 0:
                x[m] = np.random.normal()
            out[n] = np.random.normal() + np.sum(x)
    return(out / max(abs(out)), t)


# 使い方（戻り値は、x：ピンクノイズ信号、t：各サンプルの時刻）
x, t = synth_pinknoise(duration=10.0, samprate=48000)

# WAVファイルに書き出し
import soundfile as sf
sf.write("pinknoise.wav", x, 48000, subtype='PCM_24')

追記：このままだと遅いのですが、numbaを使ったら上記コードを一切変更せずに爆速になりました。JITすごい。