Comparison of pink noise pinking filters

pinking.py

#
# Compares the RBJ, PK3, and PKE pinking filters from:
# https://www.firstpr.com.au/dsp/pink-noise/
#

from scipy.signal import zpk2tf
import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as signal
from scipy.signal import butter, lfilter, freqz

rbj_gain = 0.2
pke_gain = 0.12
pk3_gain = 0.12


def apply_rbj(x):
    zeros = ( 0.98443604, 0.83392334, 0.07568359)
    poles = ( 0.99572754, 0.94790649, 0.53567505)

    b, a = zpk2tf(zeros, poles, rbj_gain)
    
    print("RBJ coefficients:")
    print(f"b = {b}")
    print(f"a = {a}")
    print("")
    
    return lfilter(b, a, x)


def apply_pke(x):
    y = np.zeros(len(x))
    
    b0 = 0
    b1 = 0
    b2 = 0
    
    for i in range(0, len(x)):
        white = x[i];

        w0 = white * 0.0990460 * pke_gain
        w1 = white * 0.2965164 * pke_gain
        w2 = white * 1.0526913 * pke_gain
        
        b0 = 0.99765 * b0 + w0
        b1 = 0.96300 * b1 + w1
        b2 = 0.57000 * b2 + w2

        pink = b0 + b1 + b2 + (white * 0.1848 * pke_gain);

        y[i] = pink;

    return y;


def apply_pk3(x):
    y = np.zeros(len(x))
    
    b0 = 0
    b1 = 0
    b2 = 0
    b3 = 0
    b4 = 0
    b5 = 0
    b6 = 0

    for i in range(0, len(x)):
        white = x[i];

        w0 = white * 0.0555179 * pk3_gain
        w1 = white * 0.0750759 * pk3_gain
        w2 = white * 0.1538520 * pk3_gain
        w3 = white * 0.3104856 * pk3_gain
        w4 = white * 0.5329522 * pk3_gain
        w5 = white * 0.0168980 * pk3_gain
        w6 = white * 0.115926  * pk3_gain
        
        b0 = 0.99886 * b0 + w0;
        b1 = 0.99332 * b1 + w1;
        b2 = 0.96900 * b2 + w2;
        b3 = 0.86650 * b3 + w3;
        b4 = 0.55000 * b4 + w4;
        b5 = -0.7616 * b5 - w5;

        pink = b0 + b1 + b2 + b3 + b4 + b5 + b6 + (white * 0.5362 * pk3_gain);
        b6 = w6;

        y[i] = pink

    return y


def plot(name, x):
    fs = 44100
    
    max_sample = max(np.max(x), -np.min(x))
    max_sample_db = 20 * np.log10(max_sample)
    
    print(f"{name} max sample at {max_sample_db:.1f} dBFS")
    
    f, Pxx = signal.welch(
        x,
        fs=fs,
        window='hann',
        nperseg=1024,
        noverlap=512,
        scaling='density'
    )

    plt.plot(f, 10 * np.log10(Pxx))


x = np.random.random_sample(44100 * 10)
x *= 2
x -= 1

rbj_y = apply_rbj(x)

plt.figure(figsize=(8, 4))

plot("White", x)
plot("RBJ", rbj_y)
plot("PKE", apply_pke(x))
plot("PK3", apply_pk3(x))

plt.xlabel("Frequency Hz")
plt.ylabel("Magnitude dB")
plt.grid(True)
plt.tight_layout()
plt.show()