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

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

2021-01-30

ピンクノイズフィルタのプログラミング
ピンクノイズの周波数特性

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

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

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))

サウンドプログラミングにおすすめの書籍はこちら。