accompanying code to colored noise: frequency domain filtering blog post at https://zeptoblog.com/2024/04/21/colored-noise-frequency-domain-filtering.html

pinkingWithInverseFFT.py

import numpy as np
from typing import Any
from nptyping import NDArray, Shape, Float, ComplexFloating
import matplotlib.pyplot as plt


def generateShapingFunction(n_pts: int,
                            exp: float
                            ) -> NDArray[Shape['Any'], Float]:
    f = np.ones(n_pts)
    for i in range(1, int(n_pts/2)+1):
        f[i] = 1/i**(exp/2)
        f[n_pts-i] = f[i]

    return f


def generateColoredNoise(n_pts: int,
                         exp: float
                         ) -> NDArray[Shape['Any'], ComplexFloating]:
    # generate white noise with 0 mean
    x = np.random.randn(n_pts)
    # x = x - np.mean(x)

    # generate the shaping function
    f = generateShapingFunction(n_pts, exp)

    # frequency domain filtering
    x_f = np.fft.fft(x)
    y_f = x_f*f
    y = np.fft.ifft(y_f)
    return y


n_pts = 2**14
n_runs = 1000
f = np.linspace(0, 1, n_pts)

fil_f = generateShapingFunction(n_pts, 1)
fil_f1p5 = generateShapingFunction(n_pts, 1.5)
fil_f2 = generateShapingFunction(n_pts, 2)
fil_f3 = generateShapingFunction(n_pts, 3)

scale_fil_f = np.sqrt(np.sum(fil_f**2)/n_pts)
scale_fil_f1p5 = np.sqrt(np.sum(fil_f1p5**2)/n_pts)
scale_fil_f2 = np.sqrt(np.sum(fil_f2**2)/n_pts)
scale_fil_f3 = np.sqrt(np.sum(fil_f3**2)/n_pts)

pr = range(1, int(n_pts))
fig, ax = plt.subplots()
ax.axvline(.5, color='black', linestyle='dashed')
ax.semilogx(f[pr], 20*np.log10(fil_f[pr]/scale_fil_f), label='1/f')
ax.semilogx(f[pr], 20*np.log10(fil_f1p5[pr]/scale_fil_f1p5), label='1/f^1.5')
ax.semilogx(f[pr], 20*np.log10(fil_f2[pr]/scale_fil_f2), label='1/f^2')
ax.semilogx(f[pr], 20*np.log10(fil_f3[pr]/scale_fil_f1p5), label='1/f^3')
ax.grid(True)
ax.set_xlim([1/n_pts, 1])
ax.set_xlabel('normalized frequency')
ax.set_ylabel('gain [dB]')
ax.legend()
plt.show()

y_fft_db_avg = np.zeros(n_pts)
y_fft_fil_f_db_avg = np.zeros(n_pts)
y_fft_fil_f1p5_db_avg = np.zeros(n_pts)
y_fft_fil_f2_db_avg = np.zeros(n_pts)
y_fft_fil_f3_db_avg = np.zeros(n_pts)
for i in range(n_runs):
    y = np.random.randn(n_pts)
    y = y - np.mean(y)

    y_fft = np.fft.fft(y)/n_pts

    y_fft_fil_f = y_fft*fil_f/scale_fil_f
    y_fft_fil_f1p5 = y_fft*fil_f1p5/scale_fil_f1p5
    y_fft_fil_f2 = y_fft*fil_f2/scale_fil_f2
    y_fft_fil_f3 = y_fft*fil_f3/scale_fil_f3

    y_fil_f = np.fft.ifft(y_fft_fil_f*n_pts)
    y_fil_f1p5 = np.fft.ifft(y_fft_fil_f1p5*n_pts)
    y_fil_f2 = np.fft.ifft(y_fft_fil_f2*n_pts)
    y_fil_f3 = np.fft.ifft(y_fft_fil_f3*n_pts)

    y_fft_db = 20*np.log10(np.abs(y_fft))
    y_fft_fil_f_db = 20*np.log10(np.abs(y_fft_fil_f))
    y_fft_fil_f1p5_db = 20*np.log10(np.abs(y_fft_fil_f1p5))
    y_fft_fil_f2_db = 20*np.log10(np.abs(y_fft_fil_f2))
    y_fft_fil_f3_db = 20*np.log10(np.abs(y_fft_fil_f3))

    y_fft_db_avg += y_fft_db/n_runs
    y_fft_fil_f_db_avg += y_fft_fil_f_db/n_runs
    y_fft_fil_f1p5_db_avg += y_fft_fil_f1p5_db/n_runs
    y_fft_fil_f2_db_avg += y_fft_fil_f2_db/n_runs
    y_fft_fil_f3_db_avg += y_fft_fil_f3_db/n_runs

plt.semilogx(f[pr], y_fft_db_avg[pr], label='input')
plt.semilogx(f[pr], y_fft_fil_f_db_avg[pr], label='1/f')
plt.semilogx(f[pr], y_fft_fil_f1p5_db_avg[pr], label='1/f^1.5')
plt.semilogx(f[pr], y_fft_fil_f2_db_avg[pr], label='1/f^2')
plt.semilogx(f[pr], y_fft_fil_f3_db_avg[pr], label='1/f^3')
plt.grid(True)
plt.xlabel('normalized frequency')
plt.ylabel('psd [dB]')
plt.legend()
plt.xlim([1/n_pts, 0.5])
plt.show()

plt.plot(y, label='input')
plt.plot(y_fil_f, label='1/f')
plt.plot(y_fil_f1p5, label='1/f^1.5')
plt.plot(y_fil_f2, label='1/f^2')
plt.plot(y_fil_f3, label='1/f^3')
plt.xlim([0, n_pts])
plt.grid(True)
plt.xlabel('sample')
plt.ylabel('amplitude')
plt.legend()
plt.show()