"""Multirate signal processing: decimation, interpolation, and rate conversion.""" import numpy as np from scipy.signal import firwin, lfilter def design_multirate_filter(L_or_M: int, num_taps: int = 63, cutoff_factor: float = 0.9) -> np.ndarray: """Design a lowpass FIR filter for multirate rate change. The cutoff is placed at pi/L_or_M (the Nyquist frequency of the lower-rate signal), pulled in slightly by cutoff_factor to leave room for the transition band. Parameters ---------- L_or_M : int Rate change factor (interpolation L or decimation M). num_taps : int Number of FIR filter coefficients (odd recommended). cutoff_factor : float Fraction of the ideal cutoff (1.0 = at Nyquist, 0.9 = 10 % guard band). Returns ------- h : np.ndarray FIR filter coefficients, length num_taps. """ cutoff = cutoff_factor / L_or_M return firwin(num_taps, cutoff) def decimate(x: np.ndarray, M: int, fir_order: int = 63) -> np.ndarray: """Decimate signal x by integer factor M with anti-aliasing FIR filter. Applies a lowpass filter with cutoff at pi/M, then keeps every M-th sample. Parameters ---------- x : np.ndarray Input signal. M : int Decimation factor (keep every M-th sample). fir_order : int FIR filter length. Returns ------- y : np.ndarray Decimated signal of length ceil(len(x) / M). """ h = design_multirate_filter(M, num_taps=fir_order) return lfilter(h, 1.0, x)[::M] def interpolate(x: np.ndarray, L: int, fir_order: int = 63) -> np.ndarray: """Interpolate signal x by integer factor L with anti-imaging FIR filter. Inserts L-1 zeros between each sample, then applies a lowpass filter with gain L to restore amplitude. Parameters ---------- x : np.ndarray Input signal. L : int Interpolation factor. fir_order : int FIR filter length. Returns ------- y : np.ndarray Interpolated signal of length len(x) * L. """ x_up = np.zeros(len(x) * L) x_up[::L] = x h = design_multirate_filter(L, num_taps=fir_order) * L return lfilter(h, 1.0, x_up) def resample_rational(x: np.ndarray, L: int, M: int, fir_order: int = 63) -> np.ndarray: """Rational sample rate conversion: interpolate by L, decimate by M. The output sample rate is fs_in * L / M. A single combined lowpass filter is used with cutoff at pi / max(L, M). Parameters ---------- x : np.ndarray Input signal. L : int Interpolation (upsample) factor. M : int Decimation (downsample) factor. fir_order : int FIR filter length. Returns ------- y : np.ndarray Resampled signal. """ # Upsample by L (zero-insert) x_up = np.zeros(len(x) * L) x_up[::L] = x # Single combined filter at the more restrictive cutoff h = design_multirate_filter(max(L, M), num_taps=fir_order) * L x_filtered = lfilter(h, 1.0, x_up) # Downsample by M return x_filtered[::M] def polyphase_decimate(x: np.ndarray, h: np.ndarray, M: int) -> np.ndarray: """Efficient polyphase decimation. Decomposes filter h into M polyphase sub-filters, each operating at the output (low) rate. This avoids computing filter outputs that would be discarded by the downsampler. The key efficiency insight: instead of filtering at the high rate and keeping every M-th output, we split the input into M branches (commutator), filter each branch with a short sub-filter, and sum. Parameters ---------- x : np.ndarray Input signal. h : np.ndarray FIR filter coefficients (anti-aliasing lowpass). M : int Decimation factor. Returns ------- y : np.ndarray Decimated signal. """ # Pad h to a multiple of M n_pad = (M - len(h) % M) % M h_padded = np.concatenate([h, np.zeros(n_pad)]) # Decompose into M polyphase sub-filters # Phase k contains coefficients h[k], h[k+M], h[k+2M], ... phases = [h_padded[k::M] for k in range(M)] # Branch k must filter the input stream x[mM - k] (the input delayed by k, # then downsampled), so prepend M zeros and take every M-th sample starting # at index M - k. Pairing phase k with x[k::M] instead is a common slip: it # filters a commutated sequence and does not equal the direct method # lfilter(h, 1, x)[::M]. This form reproduces it to floating-point precision. n_out = (len(x) + M - 1) // M x_padded = np.concatenate([np.zeros(M), x]) y = np.zeros(n_out) for k in range(M): x_branch = x_padded[M - k::M][:n_out] filtered = lfilter(phases[k], 1.0, x_branch) y[:len(filtered)] += filtered return y def cic_decimate(x: np.ndarray, M: int, N_stages: int = 1) -> np.ndarray: """CIC (cascaded integrator-comb) decimation, no multiplications. Implements H(z) = ((1 - z^{-M}) / (1 - z^{-1}))^N_stages. The integrator section runs at the high rate, the comb section at the low rate (the efficient Hogenauer decimator arrangement). Parameters ---------- x : np.ndarray Input signal. M : int Decimation factor. N_stages : int Number of CIC stages (higher = steeper rolloff, more droop). Returns ------- y : np.ndarray Decimated signal. """ y = x.astype(np.float64).copy() # Integrator stages (running sum): operate at high rate for _ in range(N_stages): y = np.cumsum(y) # Downsample # Downsample: keep every M-th sample y = y[::M] # Comb stages: operate at low rate for _ in range(N_stages): y_delayed = np.concatenate([[0], y[:-1]]) y = y - y_delayed return y def cic_response(w: np.ndarray, M: int, N_stages: int = 1) -> np.ndarray: """Magnitude response of a CIC decimator at the output rate. Evaluated as a function of the *output* digital frequency ``w`` in ``[0, pi]``. Normalised to unity DC gain, so this is the sinc^K passband droop the compensator has to undo: |H(w)| = |sin(w/2) / (M * sin(w/(2M)))|^N_stages Parameters ---------- w : np.ndarray Output-rate angular frequencies in [0, pi]. M : int Decimation factor. N_stages : int CIC order. Returns ------- mag : np.ndarray Normalised magnitude, 1.0 at DC. """ w = np.asarray(w, dtype=np.float64) mag = np.ones_like(w) nz = w > 0 mag[nz] = np.abs(np.sin(w[nz] / 2.0) / (M * np.sin(w[nz] / (2.0 * M)))) ** N_stages return mag def cic_compensator(M: int, N_stages: int = 1, num_taps: int = 9, fc: float = 0.5, grid: int = 512) -> np.ndarray: """Design an FIR that flattens CIC passband droop (inverse-sinc). The compensator runs at the CIC *output* rate and approximates ``1 / |H_CIC|`` across the passband, undoing the sinc^K droop so the cascade CIC -> compensator has a flat passband. It is a short, symmetric (linear-phase) FIR fitted by least squares. Parameters ---------- M : int CIC decimation factor. N_stages : int CIC order (number of integrator/comb stages). num_taps : int FIR length; must be odd (Type-I linear phase, exact symmetry). fc : float Passband edge as a fraction of the output Nyquist, 0 < fc < 1. Only frequencies up to ``fc * pi`` are compensated. grid : int Number of frequency points used in the least-squares fit. Returns ------- h : np.ndarray Symmetric FIR taps, length ``num_taps``, normalised to unit DC gain. """ if num_taps % 2 == 0: raise ValueError("num_taps must be odd for a Type-I linear-phase FIR") if not 0.0 < fc < 1.0: raise ValueError("fc must be a fraction of the output Nyquist in (0, 1)") w = np.linspace(0.0, fc * np.pi, grid) target = 1.0 / cic_response(w, M, N_stages) # inverse-sinc: the desired FIR shape # Type-I linear-phase FIR: A(w) = h[c] + 2 * sum_{n>=1} h[c-n] cos(n w). # Fit the half-length cosine coefficients by least squares. c = (num_taps - 1) // 2 basis = np.column_stack([np.ones_like(w)] + [2.0 * np.cos(n * w) for n in range(1, c + 1)]) coef, *_ = np.linalg.lstsq(basis, target, rcond=None) h = np.empty(num_taps) h[c] = coef[0] for n in range(1, c + 1): h[c - n] = h[c + n] = coef[n] return h / h.sum() # normalise DC gain to exactly 1