"""Lock-in detection: synchronous demodulation of a modulated signal. A lock-in amplifier measures the amplitude (and phase) of a signal at a known reference frequency by multiplying the input with a reference waveform and low-pass filtering the product. Digitally that is nothing new: over a block of samples it is exactly the coherent single-bin DFT from the estimating-a-sinusoid topic. What the lock-in *adds* is the measurement strategy around it: modulate the physical quantity you care about up to a carrier frequency where the front-end noise is white, so that the demodulated estimate obeys the white-noise averaging law instead of drowning in 1/f noise and drift at DC. The governing relation, verified in the tests against an exact Fisher information matrix built for colored (AR(1)) noise: var(A_hat) = S(f0) / T = 2 S(f0) B_n where S is the one-sided noise PSD at the carrier f0, T the integration time, and B_n the equivalent noise bandwidth of the post-demodulation low-pass (1/(2T) for a block average, 1/(4 tau) for a one-pole filter). Everything on the topic page is a consequence of reading that formula at different carrier frequencies. Frequency convention -------------------- As in the sinusoid module: frequencies are in the units of the sample rate *fs*. Leave *fs* at its default of 1.0 to work in cycles per sample; pass a real sample rate to work in Hz. """ import numpy as np def front_end_noise(n: int, sigma_white: float, corner: float, fs: float = 1.0, rng: np.random.Generator | None = None) -> np.ndarray: """Generate white + 1/f noise with a calibrated PSD. The model is the noise a real analogue front end delivers: a white floor plus a 1/f component that overtakes it below the corner frequency. The one-sided PSD is S(f) = S_w (1 + corner / f), S_w = 2 sigma_white^2 / fs so the white floor and the 1/f component are equal at ``f = corner``. Parameters ---------- n : int Number of samples (at least 2). sigma_white : float Standard deviation of the white component per sample. corner : float 1/f corner frequency, in the units of *fs*. Zero gives pure white noise. fs : float Sample rate. rng : numpy Generator, optional Source of randomness. A fresh default generator if omitted. Returns ------- ndarray of shape (n,) The noise sequence. Notes ----- The 1/f component is built by spectral shaping, like ``generate_power_law_noise_psd`` in the noise-whitening topic, but with the gain calibrated so the PSD is the stated S(f) rather than normalised to unit variance. The DC bin is set to zero, so the lowest frequency actually present is fs/n: a 1/f process has no well-defined mean, and any finite record imposes such a cutoff. """ if n < 2: raise ValueError("n must be at least 2") if sigma_white < 0 or corner < 0: raise ValueError("sigma_white and corner must be non-negative") if rng is None: rng = np.random.default_rng() white = sigma_white * rng.standard_normal(n) if corner == 0: return white # Shape unit white noise: |H(f)|^2 * (2/fs) must equal S_w * corner/f, # so H(f) = sigma_white * sqrt(corner / f). W = np.fft.rfft(rng.standard_normal(n)) f = np.fft.rfftfreq(n, d=1.0 / fs) H = np.zeros_like(f) H[1:] = sigma_white * np.sqrt(corner / f[1:]) pink = np.fft.irfft(W * H, n) return white + pink def lockin_block(x: np.ndarray, f0: float, fs: float = 1.0, reference: str = 'sine') -> tuple[float, float]: """Demodulate a whole record at the reference frequency. Multiplies the record by quadrature references at *f0* and averages: the digital dual-phase lock-in, identical to ``coherent_dft`` in the estimating-a-sinusoid topic when the reference is a sine. The magnitude/phase (rather than single-phase) output makes the result insensitive to the phase of the reference, which is what lets two free-running clocks disagree about phase without costing accuracy. Parameters ---------- x : ndarray of shape (n,) Real input samples. f0 : float Reference frequency, in the units of *fs*. For an unbiased result the record should span a whole number of reference cycles (see the coherent-integration callout on the estimating-a-sinusoid page). fs : float Sample rate. reference : {'sine', 'square'} ``'sine'`` demodulates with sampled cosine/sine references. ``'square'`` demodulates with square waves, the classic analogue lock-in reference: only additions, but the noise variance rises by pi^2/8 (about 0.9 dB) and the output responds to input at the odd harmonics 3 f0, 5 f0, ... with weights 1/3, 1/5, ... The amplitude is calibrated in both cases, so a clean on-bin tone reads (nearly) the same either way; see the Notes for the sampled-square residual. Returns ------- amplitude : float Estimated amplitude A of ``A cos(2 pi f0 t + phase)``. phase : float Estimated phase in radians, wrapped to (-pi, pi]. Notes ----- With white noise of variance sigma^2 the sine-reference amplitude estimate attains the Cramer-Rao bound ``var(A_hat) = 2 sigma^2 / N``. In terms of the one-sided PSD ``S = 2 sigma^2 / fs`` and the integration time ``T = N / fs`` this is ``S / T``, and that form generalises: for noise that is merely *locally* white around f0 the bound is ``S(f0) / T``, verified in the tests against the exact Fisher information matrix for AR(1) noise. The 1/f mountain at DC never enters, which is the entire point of modulating. A *sampled* square reference is not quite the continuous ideal. The reference transitions are placed between samples (a half-sample offset, corrected in phase internally), which keeps the pattern an exact ±1 sequence, but sampling folds the square's harmonics near multiples of the sample rate back onto the fundamental. With M samples per reference cycle the calibration residual is about 1.7/M^2 (measured: +2.6% at M = 8, +0.65% at M = 16, +0.04% at M = 64), and an M not divisible by 4 additionally biases the phase. Prefer fs/f0 divisible by 4; the sine reference has no such residual at any M. """ x = np.asarray(x, dtype=float).ravel() n = len(x) if n == 0: raise ValueError("x must not be empty") if reference not in ('sine', 'square'): raise ValueError(f"unknown reference {reference!r}") idx = np.arange(n) if reference == 'sine': z = np.mean(x * np.exp(-2j * np.pi * f0 * idx / fs)) return 2.0 * np.abs(z), float(np.angle(z)) # Square references: the fundamental of sign(cos) is (4/pi) cos, so # each quadrature product averages to (2/pi) A cos/sin(phase) and # the calibration factor is pi/2 (on top of the usual 2). The # half-sample offset keeps samples off the transitions (an exact # ±1 pattern); the resulting phase shift of pi f0/fs is undone on z. theta_half = 2 * np.pi * f0 * (idx + 0.5) / fs zi = np.mean(x * np.sign(np.cos(theta_half))) zq = np.mean(x * np.sign(np.sin(theta_half))) z = (zi - 1j * zq) * np.exp(1j * np.pi * f0 / fs) return float(np.pi / 2 * np.abs(z)), float(np.angle(z)) def lockin_stream(x: np.ndarray, f0: float, fs: float = 1.0, tau: float = 1.0) -> tuple[np.ndarray, np.ndarray]: """Track a slowly varying amplitude and phase sample by sample. The streaming form of the lock-in: multiply by quadrature references, low-pass each product with a one-pole filter of time constant *tau*, and read amplitude and phase off the filtered pair. Where :func:`lockin_block` answers "what was the amplitude over this record", this answers "what is the amplitude *now*", at the price of the usual bandwidth trade: the output follows changes on time scales slower than tau and averages away noise in an equivalent noise bandwidth of ``B_n = 1/(4 tau)``. Parameters ---------- x : ndarray of shape (n,) Real input samples. f0 : float Reference frequency, in the units of *fs*. fs : float Sample rate. tau : float Output filter time constant, in the time units implied by *fs* (samples when fs = 1, seconds when fs is in Hz). Must span at least a few reference cycles, or the 2 f0 ripple of the demodulated product leaks into the output. Returns ------- amplitude : ndarray of shape (n,) Instantaneous amplitude estimate. phase : ndarray of shape (n,) Instantaneous phase estimate in radians. Notes ----- The in-phase and quadrature outputs each carry noise variance ``2 S(f0) B_n`` (S one-sided; both sidebands of the carrier fold onto the output, hence the 2), so the amplitude noise standard deviation is ``sqrt(S(f0) / (2 tau))``. Verified in the tests. Filtering I and Q *before* forming the magnitude matters: the magnitude of a noisy pair is biased upwards (a Rayleigh floor when no signal is present), and averaging magnitudes never removes that bias, while averaging I and Q does. See the topic page's detection section. """ x = np.asarray(x, dtype=float).ravel() if len(x) == 0: raise ValueError("x must not be empty") if tau <= 0: raise ValueError("tau must be strictly positive") if tau * fs < 3: raise ValueError("tau must span at least a few samples") theta = 2 * np.pi * f0 * np.arange(len(x)) / fs b = np.exp(-1.0 / (fs * tau)) from scipy.signal import lfilter i_out = lfilter([1 - b], [1, -b], 2 * x * np.cos(theta)) q_out = lfilter([1 - b], [1, -b], -2 * x * np.sin(theta)) return np.hypot(i_out, q_out), np.arctan2(q_out, i_out) def enbw(window: np.ndarray, fs: float = 1.0) -> float: """One-sided equivalent noise bandwidth of an averaging window. The bandwidth of the ideal brick-wall low-pass that would pass the same white-noise power as the window used as a filter, with both normalised to the same DC gain: B_n = (fs / 2) * sum(w^2) / (sum(w))^2 For a plain N-sample average this is ``fs / (2 N) = 1 / (2 T)``; a Hann window is 1.5 times wider, which is the noise price of its leakage suppression. A one-pole filter with time constant tau has ``B_n = 1 / (4 tau)`` by the same integral (this function only handles finite windows; the tests pin the one-pole result). Parameters ---------- window : ndarray of shape (n,) Window coefficients (any positive scaling; B_n is scale-free). fs : float Sample rate. Returns ------- float Equivalent noise bandwidth in the units of *fs*. """ w = np.asarray(window, dtype=float).ravel() if len(w) == 0: raise ValueError("window must not be empty") s = np.sum(w) if s == 0: raise ValueError("window must have nonzero sum") return float(fs * np.sum(w**2) / (2.0 * s**2)) def amplitude_noise_std(psd_at_f0: float, T: float) -> float: """Standard deviation of a lock-in amplitude estimate. The square root of ``var(A_hat) = S(f0) / T``: one-sided noise PSD at the carrier, divided by the integration time. This is the instrument-sizing formula: pick the accuracy you need, read off the integration time it costs, before any hardware exists. It equals the Cramer-Rao bound when the noise is locally white around f0 (tests check it against the exact Fisher information matrix), so no estimator does better and the sizing is a true floor. Parameters ---------- psd_at_f0 : float One-sided noise PSD at the carrier frequency, in units^2 per unit of *fs* (e.g. V^2/Hz). T : float Integration time (seconds when the PSD is per Hz). Returns ------- float Standard deviation of the amplitude estimate, in the input's units. """ if psd_at_f0 < 0: raise ValueError("psd_at_f0 must be non-negative") if T <= 0: raise ValueError("T must be strictly positive") return float(np.sqrt(psd_at_f0 / T))