サイトロゴ

【Python】1/fゆらぎ（ピンクノイズ）のデジタルフィルタ

著者画像
Toshihiko Arai
投稿日: 2021-01-31

前回の記事では、電子回路のラグ・リードフィルタを使ってアナログ的にピンクノイズを作ってみた。今回は、プログラミングでデジタルピンクノイズを作ってみようと思う。なお、プログラミング言語はPythonで行った。

ピンクノイズフィルタのプログラミング

こちらのは、ホワイトノイズにピンクノイズフィルタを掛けてwaveファイル出力するプログラムである。IIRフィルタを使用している。正確な1/fノイズ（ピンクノイズ）フィルタのアルゴリズムは、今だに解明されていないようである。よって、こちらのサイトを参考に近似的なフィルタでピンクノイズを実現している。 https://ccrma.stanford.edu/~jos/sasp/Example_Synthesis_1_F_Noise.html

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.io.wavfile import write


def output_wav(filename, rate, data):
    write(filename, rate, data.astype(np.int16))


if __name__ == '__main__':
    rate = 44100
    T = np.arange(0, 2, 1 / rate)
    white = np.random.rand(len(T)) - 0.5
    data = white * 32767
    output_wav('white_noise.wav', rate, data.astype(np.int16))

    B = [0.049922035, -0.095993537, 0.050612699, -0.004408786]
    A = [1, -2.494956002, 2.017265875, -0.522189400]
    # a[0]*y[n] = b[0]*x[n] + b[1]*x[n-1] + ... + b[M]*x[n-M]
    #                       - a[1]*y[n-1] - ... - a[N]*y[n-N]
    pink = signal.lfilter(B, A, white) * 8

    fft_data = 10*np.log(np.abs(np.fft.rfft(pink)))  # 縦軸
    freqList = np.fft.rfftfreq(len(pink), 1.0 / rate)  # 横軸
    ax = plt.subplot()
    ax.set_xlim([10, 20000])
    ax.set_ylim([5, 100])
    plt.loglog(freqList, fft_data)
    plt.xlabel('Frequency')
    plt.ylabel('Power')
    plt.show()

    data = pink * 32767
    output_wav('pink_noise.wav', rate, data.astype(np.int16))

ピンクノイズの周波数特性

上記のプログラムを実行して出力されたWaveファイルの周波数特性を調べると図のようになった。-3dB/oct(-10dB/dec)となっていることがわかる。 https://soundcloud.com/user-999425084/pink-noise-by-iir-filter

ピンクノイズを作る方法は、これ以外にも「The Voss algorithm」「The Voss-McCartney algorithm」などが存在する。詳しくはこちらのサイトを参考に。 https://www.firstpr.com.au/dsp/pink-noise/

他にもピンクノイズのフィルタに関する海外のサイトが役立ったので紹介しておこう。 http://web.cvxr.com/cvx/examples/filter_design/html/one_over_f_filter.html https://pjb.com.au/comp/lua/digitalfilter_talk/index.html

関連記事