IAAFT for Log-Normal Blue/Violet Noise: Matching Target Power Spectral Density and Marginal Distribution

Create log-normal signals with blue/violet power spectral density using IAAFT. Prevent whitening from exponentiation by enforcing PSD and marginal distribution.

Generating log-normal signals with a rising power spectral density sounds simple on paper: draw Gaussian noise, shape its spectrum, and exponentiate. In practice, the last step derails the spectral structure for high-pass flavors like blue or violet noise. This guide walks through why that happens and how to produce signals that match both a target marginal distribution (e.g., log-normal) and a target PSD using IAAFT.

Problem overview

The goal is to synthesize random sequences whose samples follow a long-tailed distribution with support on the positive reals (for example, log-normal), while the power spectral density is low at low frequencies and increases with frequency (e.g., blue noise at +10 dB/decade or violet noise at +20 dB/decade). Drawing from a log-normal distribution is straightforward. Shaping Gaussian noise to a target PSD via the FFT is also standard. The difficulty arises when these are combined by applying a nonlinear mapping (the exponential) to an already shaped process.

Repro case: shaping first, then exponentiating

The following code generates colored Gaussian noise with several canonical PSDs, then applies an exponential mapping to obtain log-normal samples. The PSD shape is preserved for white and pink noise after exponentiation, but blue and violet noise regress toward a white-like spectrum.

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import welch
# Spectral shaping via rFFT: multiply by target magnitude spectrum
# and invert.
def shape_spectrum(n_points, psd_fn=lambda f: 1.0):
    freq_side = np.fft.rfft(np.random.randn(n_points))
    spec = psd_fn(np.fft.rfftfreq(n_points))
    spec = spec / np.sqrt(np.mean(spec ** 2))
    colored = freq_side * spec
    return np.fft.irfft(colored)
# Helper to freeze a PSD function into a generator of length->signal.
def with_psd(psd_fn):
    return lambda n: shape_spectrum(n, psd_fn)
@with_psd
def make_white(f):
    return 1.0
@with_psd
def make_blue(f):
    return np.sqrt(f)
@with_psd
def make_violet(f):
    return f
@with_psd
def make_brown(f):
    return 1.0 / np.where(f == 0, float('inf'), f)
@with_psd
def make_pink(f):
    return 1.0 / np.where(f == 0, float('inf'), np.sqrt(f))
def synthesize_colored(n_points, kind="white"):
    if kind == "white":
        return make_white(n_points)
    elif kind == "pink":
        return make_pink(n_points)
    elif kind == "blue":
        return make_blue(n_points)
    elif kind == "violet":
        return make_violet(n_points)
    elif kind == "brown":
        return make_brown(n_points)
# Demo: compare Gaussian vs. log-normalized versions across PSD types.
n_points = 100_000
kinds = ["white", "pink", "blue", "violet"]
colors = ["lightgray", "magenta", "blue", "violet"]
fig, axes = plt.subplots(2, 2, sharex=False, sharey="row")
bins_gauss = np.linspace(-5, 5, 101)
for nm, col in zip(kinds, colors):
    x = synthesize_colored(n_points, nm)
    axes[0, 0].hist(x, bins=bins_gauss, color=col, label=nm, histtype="step", density=True)
    fx, px = welch(x)
    axes[1, 0].loglog(fx, px, color=col, label=nm)
bins_logn = np.linspace(0, 10, 101)
for nm, col in zip(kinds, colors):
    y = np.exp(synthesize_colored(n_points, nm))
    axes[0, 1].hist(y, bins=bins_logn, color=col, label=nm, histtype="step", density=True)
    fy, py = welch(y)
    axes[1, 1].loglog(fy, py, color=col, label=nm)
axes[0, 0].set_title("Gaussian samples")
axes[0, 1].set_title("Exp(Gaussian) 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()

Why the spectral shape collapses after exponentiation

Exponentiation is a nonlinear operation. Unlike linear filtering, it does not preserve the linear correlation structure that the PSD encodes. In the case of high-pass spectra such as blue and violet noise, applying the exponential substantially alters the spectrum, and the resulting PSD tends to resemble white noise. This matches the observed behavior.

A separate, orthogonal point is that normalizing the shaping curve to preserve average energy, as is often done when generating colored noise in the frequency domain, does not compensate for the subsequent nonlinear distortion. The whitening effect is due to the nonlinear mapping, not the normalization.

IAAFT: enforce PSD and marginal distribution together

A robust way to meet both constraints—target marginal distribution and target PSD—is to iterate between them. The Iterative Amplitude Adjusted Fourier Transform (IAAFT) algorithm does exactly that. It alternates two projections: one step imposes the desired magnitude spectrum; the other restores the target sample distribution via rank remapping. This approach has been described for stochastic signals and is well documented in the referenced literature.

An additional implementation detail that significantly improves the fidelity of the spectral target is to construct the DFT magnitude using piecewise-linear interpolation in log–log space when the PSD is specified in log–log form. This matches how such slopes are typically defined (e.g., “+10 dB/decade”) and yields much tighter PSD control.

import numpy as np
from scipy import signal as sig
from scipy.stats import rankdata
from numpy.fft import rfft, irfft
# Configuration
n_pts = 2 ** 16
fs_rate = 1.0
max_iter = 50
eps_tol = 1e-6  # relative PSD error threshold
# Target marginal: log-normal shaped by exponentiating Gaussian draws
target_seed = np.exp(np.random.randn(n_pts))
asc_template = np.sort(target_seed)
# Target magnitude spectrum: example +10 dB/decade ("blue noise")
# Build via piecewise-linear interpolation on log–log axes.
brk_freqs = np.array([1e-5 * fs_rate, fs_rate / 2])
brk_psd_db = 10 * np.log10(brk_freqs)
freq_bins = np.fft.rfftfreq(n_pts, d=1 / fs_rate)
# Insert a first point to avoid log10(0) and set DC explicitly.
brk_freqs = np.insert(brk_freqs, 0, fs_rate / n_pts)
brk_psd_db = np.insert(brk_psd_db, 0, brk_psd_db[0])
psd_interp_db = np.interp(
    np.log10(freq_bins[1:]),
    np.log10(brk_freqs[1:]),
    brk_psd_db[1:]
)
# Convert to amplitude spectrum; set DC bin to zero for high-pass shapes.
target_amp = np.empty_like(freq_bins)
target_amp[0] = 0.0
target_amp[1:] = 10 ** (psd_interp_db / 20)
# Rank-based remapping: replace values by those of the target marginal with same ranks.
def remap_by_rank(values, sorted_ref):
    ranks = rankdata(values, method='ordinal') - 1
    return sorted_ref[ranks]
# IAAFT loop: alternate spectrum enforcement and marginal reconstruction.
synth = target_seed.copy()
prev_err = np.inf
for i in range(max_iter):
    # Impose target magnitude spectrum.
    X = rfft(synth)
    X = target_amp * np.exp(1j * np.angle(X))
    shaped = irfft(X, n=n_pts)
    # Restore target marginal via rank remapping.
    synth = remap_by_rank(shaped, asc_template)
    # Check convergence periodically.
    if (i % 5 == 0) or (i == max_iter - 1):
        cur_err = np.linalg.norm(np.abs(rfft(synth)) - target_amp) / np.linalg.norm(target_amp)
        print(f"iter {i:2d}: PSD rel-err = {cur_err:.3e}")
        if cur_err > prev_err - eps_tol:
            print("converged.")
            break
        prev_err = cur_err

How this compares to filtering-based approaches

Linear filtering of log-normal samples with blue-noise shaping filters (FIR or IIR) can achieve the desired PSD slopes. The outcome remains strictly positive and right-skewed, but it does not remain log-normal due to linear mixing of samples. In the provided results, the IAAFT approach offered the best simultaneous match of PSD and marginal distribution among the tested methods, with the added benefit of using log–log interpolation to track PSD targets defined in decades.

That said, tests showed that for a log-normal target distribution and non-white PSD, the empirical marginal can deviate and exhibit a mode at zero rather than near one, even though it remains long-tailed. This is an important caveat to keep in mind when validating results.

Why this matters

Many downstream algorithms assume specific PSD features or marginal properties. If either side is off—say, the PSD collapses after a nonlinear mapping, or the distribution drifts after filtering—benchmarks, detectors, or simulation-based analyses can mislead. Enforcing both constraints together guards against such mismatches.

Practical takeaways

If you need both a prescribed PSD and a specific sample distribution such as log-normal, an iterative scheme like IAAFT is an effective tool. It alternates between spectrum enforcement and rank-based amplitude restoration, and with log–log interpolation for PSD construction it can closely track slopes such as +10 dB/decade. When validating, examine both the PSD and the marginal distribution. Be aware that for non-white targets the empirical distribution may still deviate; verify the mode and tail behavior against expectations. For tasks where only the PSD matters, direct filtering can be simpler; for tasks where both constraints are essential, the iterative approach is more appropriate.

In short, shape spectra and restore amplitudes jointly, verify both views of the data, and prefer log–log handling when your PSD targets are specified in decades.