[Python] Generating blue noise with values sampled from a log normal distribution

Aim


I am trying to generate random signals with the following two properties:

1. 
   The values should be approximately log-normally distributed (any long-tailed distribution bounded form below with non-zero mode would do).
   
   
2. 
   The power spectral density (PSD) of the signal should have low power in the lower frequencies and increase in power with increasing frequency.
   

Current approach


Generating samples from a log-normal distribution is straightforward: sample from a normal distribution and apply the exponential function 1.

Python code:

import numpy as np
import matplotlib.pyplot as plt

size = 100_000
samples = np.exp(np.random.randn(size))
plt.hist(samples, bins=np.linspace(0, 10, 101), histtype="step", density=True)
plt.show()




Generating signals with normally distributed values and a defined power spectra is covered well in this SO answer.

The Problem


Combining both approaches works well for white noise and pink noise. However, when applying the exponential to violet or blue noise, the PSD shape is not preserved. Instead, the PSDs closely resemble the PSD of white noise.

Any insights or suggestions for alternative approaches would be much appreciated.



import numpy as np
import matplotlib.pyplot as plt

from scipy.signal import welch

# from https://stackoverflow.com/a/67127726/2912349
def noise_psd(N, psd = lambda f: 1):
X_white = np.fft.rfft(np.random.randn(N));
S = psd(np.fft.rfftfreq(N))
# Normalize S
S = S / np.sqrt(np.mean(S**2))
X_shaped = X_white * S;
return np.fft.irfft(X_shaped);

def PSDGenerator(f):
return lambda N: noise_psd(N, f)

@PSDGenerator
def get_white_noise(f):
return 1;

@PSDGenerator
def get_blue_noise(f):
return np.sqrt(f);

@PSDGenerator
def get_violet_noise(f):
return f;

@PSDGenerator
def get_brownian_noise(f):
return 1/np.where(f == 0, float('inf'), f)

@PSDGenerator
def get_pink_noise(f):
return 1/np.where(f == 0, float('inf'), np.sqrt(f))

def get_colored_noise(size, color="white"):
if color == "white":
return get_white_noise(size)
elif color == "blue":
return get_blue_noise(size)
elif color == "violet":
return get_violet_noise(size)
elif color == "brownian":
return get_brownian_noise(size)
elif color == "pink":
return get_pink_noise(size)

# --------------------------------------------------------------------------------
# Sample the different noise distributions, and determine the spectral density of each series.

size = 100_000

fig, axes = plt.subplots(2, 2, sharex=False, sharey="row")

bins = np.linspace(-5, 5, 101)

noise_types = ["white", "pink", "blue", "violet"]
colors = ["lightgray", "magenta", "blue", "violet"]

for noise_type, color in zip(noise_types, colors):
samples = get_colored_noise(size, noise_type)
axes[0, 0].hist(samples, bins=bins, color=color, label=noise_type, histtype="step", density=True)

f, p = welch(samples)
axes[1, 0].loglog(f, p, color=color, label=noise_type)

# --------------------------------------------------------------------------------
# The samples above are normally distributed. Applying the exponential function yields lognormally distributed samples.
# However, the shape of the power spectrum is not preserved in every case.
# Specifically, noise with a low power in the lower frequencies and high power in high frequencies (blue, violet) appears indistinguishable from white noise.
#
# [1] https://en.wikipedia.org/wiki/Log-normal_distribution

bins = np.linspace(0, 10, 101)

for noise_type, color in zip(noise_types, colors):
samples = np.exp(get_colored_noise(size, noise_type))
axes[0, 1].hist(samples, bins=bins, color=color, label=noise_type, histtype="step", density=True)

f, p = welch(samples)
axes[1, 1].loglog(f, p, color=color, label=noise_type)

# --------------------------------------------------------------------------------
# decorate axes

axes[0, 0].set_title("Normally distributed samples")
axes[0, 1].set_title("Lognormally distributed samples")

axes[0, 0].set_ylabel("Density")
axes[0, 0].set_xlabel("Value")
axes[0, 1].set_xlabel("Value")

axes[1, 0].set_ylabel("Power")
axes[1, 0].set_xlabel("Frequency [Hz]")
axes[1, 1].set_xlabel("Frequency [Hz]")

axes[0, 0].legend(loc="upper left")

fig.tight_layout()
plt.show()

Continue reading...