""" Random process generators and characterisation tools. Every noise type is the output of a random process. This module implements generators for the major process classes encountered in DSP and a ``characterize()`` function that estimates the statistical signature of an unknown noise sequence: distribution moments, autocorrelation shape, and power-law slope. Process classes --------------- * Gaussian (white and AR-coloured): thermal noise, Johnson-Nyquist * Poisson: shot noise, photon counting, discrete arrival events * Markov / autoregressive: the workhorse for correlated noise; AR(p) driven by white noise * Lévy alpha-stable: impulse noise with heavy tails; atmospheric, power-line, ignition interference * Cyclostationary: clock-coupled noise whose *statistics* are periodic * Telegraph (random switching): burst / RTS noise in semiconductors References ---------- * Papoulis, A. & Pillai, S. U. (2002). *Probability, Random Variables, and Stochastic Processes* (4th ed.). McGraw-Hill. * Gardiner, C. W. (2004). *Handbook of Stochastic Methods* (3rd ed.). Springer. * Kasdin, N. J. (1995). Discrete simulation of colored noise and stochastic processes and 1/f^alpha power law noise generation. *Proc. IEEE*, 83(5). * Nolan, J. P. (2020). *Univariate Stable Distributions*. Springer. * Machlup, S. (1954). Noise in semiconductors: spectrum of a two-parameter random signal. *J. Appl. Phys.*, 25(3), 341-343. """ import numpy as np from scipy.signal import lfilter, welch from scipy.stats import kurtosis, skew # --------------------------------------------------------------------------- # Gaussian process generators # --------------------------------------------------------------------------- def gaussian_white(n: int, sigma: float = 1.0, seed: int | None = None) -> np.ndarray: """Generate white Gaussian noise (thermal noise model). Parameters ---------- n : int Number of samples. sigma : float Standard deviation (RMS amplitude). seed : int or None RNG seed for reproducibility. Returns ------- x : ndarray of shape (n,) Zero-mean, uncorrelated Gaussian samples. Notes ----- White Gaussian noise models thermal (Johnson-Nyquist) noise in resistors and the aggregate of many small independent noise sources via the central limit theorem. Its ACF is :math:`R_{xx}[k] = \\sigma^2 \\delta[k]` and its PSD is flat at :math:`\\sigma^2 / f_s` [W/Hz]. """ rng = np.random.default_rng(seed) return rng.normal(0.0, sigma, size=n) def gaussian_ar1(n: int, phi: float = 0.9, sigma: float = 1.0, seed: int | None = None) -> np.ndarray: """Generate AR(1) Gaussian noise (Ornstein-Uhlenbeck discrete analogue). Parameters ---------- n : int Number of samples. phi : float AR(1) coefficient. |phi| < 1 for stationarity. sigma : float Standard deviation of the driving white noise. seed : int or None RNG seed. Returns ------- x : ndarray of shape (n,) Correlated Gaussian samples with exponential ACF decay. """ rng = np.random.default_rng(seed) w = rng.normal(0.0, sigma, size=n) # lfilter uses a[0]*y[n] = b[0]*x[n] + b[1]*x[n-1] - a[1]*y[n-1] # AR(1): x[n] = w[n] + phi * x[n-1] => x[n] - phi*x[n-1] = w[n] # AR(1) driven by noise of std sigma has variance sigma^2 / (1 - phi^2). x = lfilter([1.0], [1.0, -phi], w) return x # --------------------------------------------------------------------------- # Poisson process generator # --------------------------------------------------------------------------- def poisson_process(n: int, rate: float = 10.0, seed: int | None = None) -> np.ndarray: """Generate a Poisson (shot) noise process. Models the arrival of discrete independent events at a constant average rate. Each sample is drawn from Poisson(*rate*). Shot noise in photodiodes and transistor collector current follows this process; at high rates it approximates a Gaussian with variance equal to the mean. Parameters ---------- n : int Number of samples. rate : float Expected event count per sample interval (lambda). seed : int or None RNG seed. Returns ------- x : ndarray of shape (n,), dtype int Poisson-distributed counts. """ rng = np.random.default_rng(seed) return rng.poisson(rate, size=n) # --------------------------------------------------------------------------- # AR(p) / Markov process generator # --------------------------------------------------------------------------- def ar_process(n: int, coeffs: np.ndarray | list, sigma: float = 1.0, seed: int | None = None) -> np.ndarray: """Generate an autoregressive AR(p) process. The AR(p) model is :math:`x[t] = w[t] - \\sum_{i=1}^p a_i x[t-i]` where *w* is white noise with standard deviation *sigma*. This is the workhorse for correlated noise: by choosing the coefficients you control the ACF shape and PSD slope. Parameters ---------- n : int Number of samples. coeffs : array-like of length p AR coefficients a_1 .. a_p (the denominator of the all-pole filter). sigma : float Standard deviation of the driving white noise. seed : int or None RNG seed. Returns ------- x : ndarray of shape (n,) Correlated noise samples. """ coeffs = np.asarray(coeffs, dtype=float) rng = np.random.default_rng(seed) w = rng.normal(0.0, sigma, size=n) # AR(p): x[t] + a_1 x[t-1] + ... + a_p x[t-p] = w[t] a = np.r_[1.0, coeffs] # a[0] = 1, then a_1 .. a_p x = lfilter([1.0], a, w) return x # --------------------------------------------------------------------------- # Lévy alpha-stable / impulse noise generator # --------------------------------------------------------------------------- def levy_impulse_noise(n: int, scale: float = 5.0, outlier_fraction: float = 0.02, seed: int | None = None) -> np.ndarray: """Generate impulse noise via a contaminated Gaussian mixture. True alpha-stable distributions (Lévy flights) lack closed-form PDFs except at alpha=0.5 (the one-sided Lévy distribution itself), alpha=1 (Cauchy), and alpha=2 (Gaussian). For DSP work a contaminated Gaussian mixture captures the essential behaviour (heavy tails that break the 3-sigma assumption) and is closed-form for verification. The output is x[t] = (1-p)*N(0, 1) + p*N(0, scale) where p = *outlier_fraction*. Parameters ---------- n : int Number of samples. scale : float Ratio of outlier std to nominal std (>= 1; default 5.0 = 5× larger outliers). Larger = heavier tails. outlier_fraction : float Probability of drawing from the outlier component (default 0.02). seed : int or None RNG seed. Returns ------- x : ndarray of shape (n,) Impulse-noise samples with occasional large excursions. """ rng = np.random.default_rng(seed) mask = rng.random(size=n) < outlier_fraction x = rng.normal(0.0, 1.0, size=n) x[mask] = rng.normal(0.0, scale, size=n)[mask] return x # --------------------------------------------------------------------------- # Cyclostationary process generator # --------------------------------------------------------------------------- def cyclostationary_noise(n: int, period: int = 100, modulation_depth: float = 0.5, seed: int | None = None) -> np.ndarray: """Generate cyclostationary noise: white noise with periodic variance. Models clock-coupled noise in digital systems: the noise *statistics* (here the variance) repeat with a known period rather than being constant. Common in switched-capacitor circuits, ADC clock feedthrough, and processor supply noise that follows the instruction cycle. Parameters ---------- n : int Number of samples. period : int Period of the variance modulation in samples. modulation_depth : float Fraction of the nominal std that the modulation adds. 0 = stationary white noise; 1 = std swings from 0 to 2× nominal (variance from 0 to 4× nominal). seed : int or None RNG seed. Returns ------- x : ndarray of shape (n,) Cyclostationary noise samples. """ rng = np.random.default_rng(seed) t = np.arange(n) envelope = 1.0 + modulation_depth * np.sin(2 * np.pi * t / period) return rng.normal(0.0, 1.0, size=n) * envelope # --------------------------------------------------------------------------- # Telegraph process (RTS / burst noise) generator # --------------------------------------------------------------------------- def telegraph_noise(n: int, p_switch: float = 0.005, levels: tuple = (-1.0, 1.0), seed: int | None = None) -> np.ndarray: """Generate random telegraph signal (RTS / burst noise). Models a two-state random process where the signal switches between two discrete levels with a fixed per-sample probability. In semiconductors this is the signature of a single charge trap filling and emptying (random telegraph signal noise in MOSFETs). The ACF decays exponentially with time constant 1 / (2 * p_switch). Parameters ---------- n : int Number of samples. p_switch : float Probability of switching state at each sample (0 < p_switch < 1). levels : tuple (low, high) The two amplitude levels. seed : int or None RNG seed. Returns ------- x : ndarray of shape (n,) Two-level random telegraph signal. """ rng = np.random.default_rng(seed) x = np.empty(n) state = rng.choice(levels) switches = rng.random(size=n) < p_switch for i in range(n): if switches[i]: state = levels[0] if state == levels[1] else levels[1] x[i] = state return x # --------------------------------------------------------------------------- # Characterisation # --------------------------------------------------------------------------- def characterize(x: np.ndarray, fs: float = 1.0) -> dict: """Estimate the statistical signature of a noise sequence. Computes distribution moments, normalised autocorrelation at key lags, and a PSD slope estimate via log-log linear regression on Welch's periodogram. Returns a dictionary suitable for comparing process classes or verifying a noise generator. Parameters ---------- x : ndarray Noise samples. fs : float Sample rate (used only for the frequency axis in the PSD fit). Returns ------- info : dict Keys: 'mean', 'var', 'std', 'skewness', 'excess_kurtosis', 'acf_lag1', 'acf_lag2', 'acf_lag5', 'psd_slope', 'psd_r2', 'is_zero_mean' (bool), 'is_symmetric' (bool). """ info = {} info['mean'] = float(np.mean(x)) info['var'] = float(np.var(x)) info['std'] = float(np.std(x, ddof=0)) info['skewness'] = float(skew(x)) info['excess_kurtosis'] = float(kurtosis(x)) # excess kurtosis = kurtosis - 3 info['is_zero_mean'] = bool(abs(info['mean']) < 0.05 * info['std']) # Normalised autocorrelation at key lags. xc = x - info['mean'] acf = np.correlate(xc, xc, mode='full') acf = acf / acf[len(xc) - 1] # normalise so R[0] = 1 mid = len(xc) - 1 N = len(x) for lag, key in [(1, 'acf_lag1'), (2, 'acf_lag2'), (5, 'acf_lag5')]: info[key] = float(acf[mid + lag]) if lag < N else np.nan # PSD slope via Welch + log-log regression. nperseg = min(256, len(x) // 4) if nperseg >= 16: f, psd = welch(x, fs, nperseg=nperseg) # Exclude DC and fit above f=0. mask = f > 0 if mask.sum() > 4: log_f = np.log10(f[mask]) log_p = np.log10(psd[mask]) coeffs = np.polyfit(log_f, log_p, 1) info['psd_slope'] = float(coeffs[0]) # R-squared. pred = np.polyval(coeffs, log_f) ss_res = np.sum((log_p - pred) ** 2) ss_tot = np.sum((log_p - np.mean(log_p)) ** 2) info['psd_r2'] = float(1 - ss_res / ss_tot) if ss_tot > 0 else 0.0 else: info['psd_slope'] = np.nan info['psd_r2'] = np.nan else: info['psd_slope'] = np.nan info['psd_r2'] = np.nan # Symmetry check: |skewness| < 0.5 info['is_symmetric'] = bool(abs(info['skewness']) < 0.5) return info