Wiener Filter: Optimal Linear Filtering Theory and Python Implementation — SNR-Based Noise Reduction Design Example

Introduction

The Wiener filter is the optimal linear filter that minimizes mean squared error (MSE). For estimating a signal from noisy observations, it provides the theoretically optimal solution within the framework of stationary stochastic processes.

Problem Formulation

Noisy observation:

where \(s(t)\) is the desired signal and \(n(t)\) is noise. We estimate \(s(t)\) by filtering:

MSE Minimization and the Full Derivation of the Wiener-Hopf Equation

Minimize the mean squared estimation error:

Equation (2) is an infinite-dimensional convolution integral, but in practice we approximate it with a finite-tap FIR filter. Here we derive, directly from the gradient of the discrete, finite-dimensional vector formulation, exactly why the optimal solution takes the form of a normal equation (the Wiener-Hopf equation).

Vector Formulation

We stack the observed signal into a tapped delay line of length \(M\) :

With filter coefficient vector \(\mathbf{w} = [w_0, w_1, \ldots, w_{M-1}]^\top\) , the estimate is

which is exactly the discrete, \(M\) -tap truncation of the integral in Equation (2). The MSE cost function is

Expansion and the Zero-Gradient Condition

Expanding the square:

Defining the autocorrelation matrix \(\mathbf{R}_{xx} = E[\mathbf{x}[n]\mathbf{x}[n]^\top]\) (the matrix version of \(R_{yy}\) in this article), the cross-correlation vector \(\mathbf{r}_{xd} = E[\mathbf{x}[n]s[n]]\) (the matrix version of \(R_{sy}\) ), and the signal variance \(\sigma_s^2 = E[s[n]^2]\) , the cost function becomes a quadratic form in \(\mathbf{w}\) :

Differentiating with respect to \(\mathbf{w}\) :

Setting the gradient to zero yields the normal equation — the discrete version of the Wiener-Hopf equation:

Generalizing this to continuous time with infinitely many taps gives the convolution equation derived from the orthogonality principle:

where \(R_{yy}\) is the autocorrelation of the observed signal and \(R_{sy}\) is the cross-correlation between desired and observed signals. Equations (6) and (7) express the same optimality condition (the orthogonality principle: the estimation error \(s[n]-\hat{s}[n]\) must be uncorrelated with every component of the observation \(\mathbf{x}[n]\) ) in finite- and infinite-dimensional form respectively.

Convexity of \(J(\mathbf{w})\) and Positive Semi-Definiteness of \(\mathbf{R}_{xx}\)

To confirm Equation (6) truly minimizes the MSE, we need to show \(J(\mathbf{w})\) is convex. Differentiating the gradient once more with respect to \(\mathbf{w}\) gives the Hessian:

For any vector \(\mathbf{a} \neq \mathbf{0}\) :

so \(\mathbf{R}_{xx}\) is always positive semi-definite (the expectation of a squared quantity is non-negative). Hence the Hessian \(2\mathbf{R}_{xx}\) is also positive semi-definite, and \(J(\mathbf{w})\) is a convex function of \(\mathbf{w}\) . Since a stationary point of a convex function is a global minimum, any \(\mathbf{w}\) satisfying Equation (6) is guaranteed to globally minimize the MSE.

Furthermore, if \(\mathbf{R}_{xx}\) is positive definite (all eigenvalues strictly positive), the minimizer is unique. Conversely, if the components of \(\mathbf{x}[n]\) are linearly dependent (a pure sinusoid, discussed below, is the canonical example), \(\mathbf{R}_{xx}\) becomes singular, and the solution is either non-unique or numerically unstable.

Frequency-Domain Solution: Derivation

Equation (7) is a time-domain convolution equation. Taking the Fourier transform of both sides, the convolution theorem turns it into a product:

By the Wiener-Khinchin theorem, the Fourier transform of an autocorrelation (or cross-correlation) function equals its power spectral density (or cross spectral density):

so Equation (7) decomposes, in the frequency domain, into an independent scalar equation at each frequency bin:

Dividing both sides by \(S_{yy}(f)\) :

The normal equation, which was an \(M\) -dimensional (or infinite-dimensional) linear system in the time domain, reduces to a per-frequency scalar division in the frequency domain — this is precisely why the Wiener filter has such a clean closed form in frequency.

When signal and noise are uncorrelated:

where \(\text{SNR}(f) = S_{ss}(f) / S_{nn}(f)\) . This follows immediately from Equation (8) since \(y=s+n\) with \(E[sn]=0\) gives \(S_{yy}(f)=S_{ss}(f)+S_{nn}(f)\) and \(S_{sy}(f)=S_{ss}(f)\) .

Intuition: Frequencies with high SNR pass through; frequencies with low SNR are attenuated.

Numerical Experiment: Solving the Wiener-Hopf Equation and Verifying MMSE Optimality

We now verify, in numpy, the claim behind Equations (5)-(6): that solving the normal equation genuinely minimizes the MSE. Using the same signal as the implementation below (\(fs=1000\) , a \(5\,\mathrm{Hz}+20\,\mathrm{Hz}\) sinusoid mixture, additive white noise with standard deviation \(0.8\) , np.random.seed(42)), we build an \(M=8\) -tap delay line \(\mathbf{x}[n]\) and desired signal \(d[n]=s[n]\) (the clean signal), compute the empirical \(\mathbf{R}_{xx}\) and \(\mathbf{r}_{xd}\) , and solve Equation (6) with np.linalg.solve.

import numpy as np

np.random.seed(42)
fs = 1000
t = np.arange(0, 1, 1/fs)
signal = np.sin(2*np.pi*5*t) + 0.5*np.sin(2*np.pi*20*t)
noise = np.random.normal(0, 0.8, len(t))
observed = signal + noise

M = 8
N = len(observed)
X = np.zeros((N - M + 1, M))
for i in range(M):
    X[:, i] = observed[M-1-i : N-i]
d = signal[M-1:]

R = (X.T @ X) / X.shape[0]
r = (X.T @ d) / X.shape[0]
w_opt = np.linalg.solve(R, r)

The results:

• Solution \(\mathbf{w}_{\text{opt}} = [0.1507,\ 0.1398,\ 0.1276,\ 0.1134,\ 0.1054,\ 0.0890,\ 0.0744,\ 0.0599]\) , minimum cost \(J(\mathbf{w}_{\text{opt}}) = 0.088442\)
• Zero-gradient check: \(\|\nabla J(\mathbf{w}_{\text{opt}})\|_\infty = 2.22\times10^{-16}\) (effectively zero, at double-precision round-off level)
• Positive-definiteness of \(\mathbf{R}_{xx}\) : minimum eigenvalue \(0.5708\) , maximum \(5.7854\) (all positive \(\Rightarrow\) positive definite \(\Rightarrow\) unique solution)
• Direct optimality check: over \(2{,}000\) trials of Gaussian perturbations with standard deviation \(0.05\) added to \(\mathbf{w}_{\text{opt}}\) , every single trial (100.00%) satisfied \(J(\mathbf{w}_{\text{opt}}+\boldsymbol{\delta}) \ge J(\mathbf{w}_{\text{opt}})\) , with the smallest cost increase still \(1.69\times10^{-3} > 0\)

The numerically zero gradient, the positive-definite \(\mathbf{R}_{xx}\) , and the fact that no perturbation direction decreases the cost together confirm, on real data, the proof above that the solution of Equation (6) is a global minimum.

Python Implementation

1D Signal Denoising

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, ifft, fftfreq

np.random.seed(42)

# --- Signal generation ---
fs = 1000
t = np.arange(0, 1, 1/fs)
signal = np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 20 * t)
noise = np.random.normal(0, 0.8, len(t))
observed = signal + noise

# --- Wiener filter (frequency domain) ---
N = len(observed)
Y = fft(observed)
freqs = fftfreq(N, 1/fs)

# Power spectrum estimation
S_yy = np.abs(Y)**2 / N
noise_power = np.var(noise) * N
S_nn = noise_power / N * np.ones(N)
S_ss = np.maximum(S_yy - S_nn, 0)

# Wiener filter
H = S_ss / (S_ss + S_nn + 1e-10)
S_filtered = ifft(Y * H).real

# --- Visualization ---
fig, axes = plt.subplots(3, 1, figsize=(12, 8))

axes[0].plot(t, signal, 'b-', alpha=0.7, label='Original')
axes[0].set_title('Original Signal')
axes[0].set_ylabel('Amplitude')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

axes[1].plot(t, observed, 'gray', alpha=0.5, label='Noisy')
axes[1].set_title('Observed (Signal + Noise)')
axes[1].set_ylabel('Amplitude')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

axes[2].plot(t, signal, 'b-', alpha=0.3, label='Original')
axes[2].plot(t, S_filtered, 'r-', alpha=0.8, label='Wiener filtered')
axes[2].set_title('Wiener Filter Output')
axes[2].set_xlabel('Time [s]')
axes[2].set_ylabel('Amplitude')
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

mse_noisy = np.mean((signal - observed)**2)
mse_filtered = np.mean((signal - S_filtered)**2)
print(f"MSE (noisy):    {mse_noisy:.4f}")
print(f"MSE (filtered): {mse_filtered:.4f}")
print(f"Improvement:    {(1 - mse_filtered/mse_noisy)*100:.1f}%")

Running this prints MSE (noisy): 0.6133, MSE (filtered): 0.1284, Improvement: 79.1%. Under the idealized assumption of a known noise variance (the true value is used), the frequency-domain Wiener filter cuts the MSE by roughly 80%.

Frequency Response of the Filter

plt.figure(figsize=(10, 4))
positive_freqs = freqs[:N//2]
plt.plot(positive_freqs, H[:N//2], 'g-', linewidth=2)
plt.xlabel('Frequency [Hz]')
plt.ylabel('|H(f)|')
plt.title('Wiener Filter Frequency Response')
plt.grid(True, alpha=0.3)
plt.xlim(0, 100)
plt.tight_layout()
plt.show()

Causality Constraints: Non-Causal vs Causal Wiener Filters, and the Paley-Wiener Condition

The \(H(f)\) obtained from Equations (8)-(9) is, in general, non-causal (it uses future observation samples too). The solution of Equation (8) only needs to satisfy an algebraic equation in the frequency domain, with no guarantee that the corresponding impulse response \(h(\tau)\) vanishes for \(\tau<0\) . Indeed, this article’s Python implementation (batch fft → multiply → ifft) uses the entire signal at once and is therefore non-causal — unusable for real-time processing.

For real-time applications, a separate causal Wiener filter (\(h(\tau)=0\) for \(\tau<0\) ) must be constructed. The standard procedure is spectral factorization:

1. Factor \(S_{yy}(f)\) using a minimum-phase (causally invertible) factor \(L(f)\) as \(S_{yy}(f) = \sigma^2\, L(f)\, L^\ast(f)\) .
2. The causal Wiener filter is then

where \([\cdot]_+\) denotes keeping only the non-negative-time (\(\tau \ge 0\) ) part of the inverse Fourier transform.

The necessary and sufficient condition for this spectral factorization to exist is the Paley-Wiener condition:

Intuitively, this holds unless \(S_{yy}(f)\) is a pathological spectrum that vanishes identically over some band (i.e., is perfectly silent over a finite interval). In practice, any additive noise gives \(S_{yy}(f) = S_{ss}(f)+S_{nn}(f) \ge S_{nn}(f) > 0\) , so the condition is usually satisfied automatically. What matters more in practice than the condition itself is the consequence: a causal filter always has MSE greater than or equal to the non-causal solution. The non-causal filter can only do better (or tie), since it has access to future information, and the causal filter is an approximation that accepts this degradation. The Kalman filter, as we show next, is precisely the recursive, state-space procedure for computing this causal Wiener filter.

Ill-Conditioned Autocorrelation Matrices and Regularization

As noted in the convexity discussion, when the components of \(\mathbf{x}[n]\) are nearly linearly dependent, \(\mathbf{R}_{xx}\) becomes singular or ill-conditioned, and the solution \(\mathbf{w}=\mathbf{R}_{xx}^{-1}\mathbf{r}_{xd}\) becomes unstable. A canonical ill-conditioned case is a deterministic sinusoidal signal. A pure sinusoid has only two degrees of freedom (amplitude and phase), so for an FIR filter with \(M>2\) taps the autocorrelation matrix theoretically has rank at most 2.

import numpy as np

fs = 1000
t = np.arange(0, 1, 1/fs)
pure_sine = np.sin(2*np.pi*10*t)  # noiseless pure sinusoid

M = 5  # tap count (> 2)
N = len(pure_sine)
X = np.zeros((N - M + 1, M))
for i in range(M):
    X[:, i] = pure_sine[M-1-i : N-i]
R = (X.T @ X) / X.shape[0]
d = pure_sine[M-1:]
r = (X.T @ d) / X.shape[0]

cond = np.linalg.cond(R)
eigvals = np.linalg.eigvalsh(R)

Running this gives \(\text{cond}(\mathbf{R}_{xx}) = 7.13\times10^{17}\) (an extreme ill-conditioning near the limit of double precision), with eigenvalues \([\,{-0.0},\ 0.0,\ 0.0,\ 0.0196,\ 2.4902\,]\) — exactly three eigenvalues near zero (with tiny negative values from floating-point error), matching the theoretical rank-2 prediction. Calling np.linalg.solve(R, r) on this matrix actually raises LinAlgError: Singular matrix. Using np.linalg.lstsq instead yields the minimum-norm solution \(\mathbf{w}=[0.5981,\ 0.3998,\ 0.2000,\ -0.0006,\ -0.2013]\) (\(\|\mathbf{w}\|=0.7734\) ), but this solution is sensitive to tiny perturbations of the observed data. Indeed, adding negligible noise (standard deviation \(10^{-6}\) ) to the observed signal before rebuilding \(\mathbf{R}_{xx}\) and \(\mathbf{r}_{xd}\) (100 trials) changes the solution by \(\|\Delta\mathbf{w}\|\) on the order of \(5.5\text{–}7.5\times10^{-5}\) — an amplification of the input perturbation by roughly a factor of 60.

The remedy is a Ridge-regression-style regularization: add a small positive constant \(\lambda\) to the diagonal of \(\mathbf{R}_{xx}\) and solve using \(\mathbf{R}_{xx}+\lambda I\) .

Strengthening the regularization to \(\lambda=10^{-2}\) improves the condition number to \(2.5\times10^{2}\) , and the solution’s sensitivity to the same \(10^{-6}\) -scale perturbation drops to \(2.8\times10^{-9}\text{–}4.4\times10^{-8}\) — even smaller than the input perturbation itself (more than a three-order-of-magnitude improvement in sensitivity compared to \(\lambda=10^{-8}\) ). The practical takeaway: check the condition number of \(\mathbf{R}_{xx}\) , and when ill-conditioning is suspected (near-periodic signals, or a tap count exceeding the signal’s effective degrees of freedom), choose a regularization parameter \(\lambda\) that balances stability against performance degradation.

Impact of Noise Statistics Estimation Error on Filter Performance

Equation (9)’s \(H(f)=S_{ss}/(S_{ss}+S_{nn})\) assumes the true power spectra, but in practice \(S_{nn}\) (noise power) must be estimated from finite data, and estimation error directly affects filter performance. Using a broadband signal (an AR(1) process \(s_k = 0.9\,s_{k-1}+w_k\) , \(w_k\sim\mathcal{N}(0,0.05)\) ) plus white noise with variance \(0.25\) , we deliberately shifted the assumed noise-variance estimate away from the true value when constructing the frequency-domain Wiener filter, and measured the resulting MSE against the true signal.

Underestimation (\(-50\%\) to \(-30\%\) ) monotonically worsens the MSE, as expected: underestimating the noise power inflates \(H(f)\) , letting residual noise pass through. What is surprising is that overestimating the noise variance (\(+30\%\) to \(+200\%\) ) actually improves MSE relative to using the true value. The cause is that the single-snapshot periodogram \(S_{yy}=|Y|^2/N\) itself is a high-variance estimator whose variance does not shrink with data length (see Window Functions and Power Spectral Density (PSD): Theory and Python Implementation for details). Mildly over-estimating the noise floor — an “over-subtraction” — suppresses the gain variability induced by this periodogram noise, a well-known heuristic in spectral-subtraction-based speech enhancement. But as the overestimation grows further (\(+300\%\) onward), bias (erroneously attenuating true signal energy) starts to dominate and MSE degrades monotonically. The experiment shows a U-shaped curve with a minimum near \(+200\%\) (\(0.076315\) ), quantitatively confirming the practical lesson that using the true noise power is not always optimal when the underlying spectral estimate itself is noisy.

Relationship with Kalman Filter

The Wiener filter is optimal for stationary processes; the Kalman filter extends this to non-stationary, recursive settings. The two are not unrelated theories — the steady-state Kalman filter is exactly the causal Wiener filter. We verify this numerically using the state-space model \(s_k = a\,s_{k-1}+w_k\) (\(w_k\sim\mathcal{N}(0,q)\) ), observation \(z_k = s_k+v_k\) (\(v_k\sim\mathcal{N}(0,r)\) , \(a=0.9,\ q=1.0,\ r=2.0\) ).

First, we solve the discrete algebraic Riccati equation with scipy.linalg.solve_discrete_are to obtain the steady-state Kalman gain.

from scipy.linalg import solve_discrete_are
import numpy as np

a, q, r = 0.9, 1.0, 2.0
P_pred = solve_discrete_are(np.array([[a]]), np.array([[1.0]]),
                            np.array([[q]]), np.array([[r]]))[0, 0]
K_inf = P_pred / (P_pred + r)
P_filt = (1 - K_inf) * P_pred

This gives \(P_{k|k-1}=1.757791\) (predicted error variance), steady-state Kalman gain \(K_\infty=0.467772\) , and updated error variance \(P_{k|k}=0.935545\) .

Next, from the true autocorrelations implied by this state-space model — \(R_{ss}(k)=\frac{q}{1-a^2}a^{|k|}\) , \(R_{zz}(k)=R_{ss}(k)+r\,\delta_{k,0}\) , \(R_{sz}(k)=R_{ss}(k)\) — we solve the causal FIR Wiener filter (Equation (6), \(\mathbf{R}_{xx}\mathbf{w}=\mathbf{r}_{xd}\) ) for increasing tap counts \(M\) , computing the theoretical MMSE \(=\sigma_s^2-\mathbf{w}^\top\mathbf{r}_{xd}\) at each step.

As the tap count grows, MMSE(\(M\) ) decreases monotonically and converges to \(0.935545\) from \(M=16\) onward — matching the steady-state Kalman filter’s updated error variance \(P_{k|k}=0.935545\) to six decimal places. Moreover, the converged FIR filter’s leading coefficient \(w_0=0.467772\) exactly matches the Kalman gain \(K_\infty=0.467772\) , and the tap ratio \(w_1/w_0=0.479005\) exactly matches the equivalent pole of the Kalman recursion \((1-K_\infty)a=0.479005\) . This is expected theoretically: expanding the steady-state Kalman recursion \(\hat{s}_{k|k}=(1-K_\infty)a\,\hat{s}_{k-1|k-1}+K_\infty z_k\) recursively gives

which is exactly an infinite-tap causal Wiener filter whose coefficients form the geometric sequence \(w_j = K_\infty\bigl[(1-K_\infty)a\bigr]^j\) .

Finally, we generated 30 independent series of \(N=5000\) samples each (np.random.default_rng(seed), \(\text{seed}=0,\ldots,29\) ) and applied both the steady-state Kalman filter and the \(M=40\) -tap FIR Wiener filter to the same data, comparing the Monte Carlo average MSE. Both filters gave identical results, \(0.934818 \pm 0.003926\) (standard error), consistent with the theoretical value \(P_{k|k}=0.935545\) within one standard error. This confirms, both analytically and numerically, that “the Kalman gain converges to the Wiener filter in steady state.”

Adaptive Wiener Filtering: LMS and RLS as Sequential Approximations

Everything above assumed \(\mathbf{R}_{xx}\) and \(\mathbf{r}_{xd}\) (or \(S_{yy}(f)\) , \(S_{sy}(f)\) ) are known. In real environments, however, signal and noise statistics are often unknown or non-stationary (time-varying), making the batch solution of Equation (6) infeasible. Adaptive Filters (LMS/RLS): Theory and Python Implementation , LMS and NLMS Adaptive Algorithms: Theory and Python Implementation , and RLS (Recursive Least Squares) Algorithm: Theory and Python Implementation detail how LMS and RLS iteratively approach the solution of Equation (6) using only observed samples — positioning them as practical means of realizing this article’s Wiener solution without explicitly estimating the underlying statistics.

• LMS updates coefficients using the instantaneous gradient \(\nabla J \approx -2e[n]\mathbf{x}[n]\) (a stochastic approximation using a single-sample product in place of the true expectation), and converges linearly, on average, to the Wiener solution \(\mathbf{w}^\ast=\mathbf{R}^{-1}\mathbf{p}\) .
• RLS updates the exact normal equation \(\mathbf{R}(n)\mathbf{w}(n)=\mathbf{r}(n)\) (the same form as Equation (6)) at each time step in \(O(M^2)\) via the matrix inversion lemma (Sherman-Morrison formula), achieving quadratic convergence — reaching the normal equation’s solution in the shortest possible time, independent of the input statistics. Measured results show the number of samples needed to reach twice the noise floor’s MSE is 57 for LMS versus 27 for RLS — roughly half.

In other words, the \(\mathbf{R}_{xx}\mathbf{w}=\mathbf{r}_{xd}\) solved in this article gives the “ground truth” answer only when the statistics are known in advance; LMS and RLS are practical algorithms that approximate this ground truth online in unknown environments. When stationarity breaks down (i.e., \(\mathbf{R}_{xx}\) or \(\mathbf{r}_{xd}\) vary over time), RLS with a forgetting factor, or NLMS, becomes effective — adaptive filters fill the gap left by the Wiener filter’s “known, stationary statistics” assumption using a data-driven approach, distinct from the model-based approach of the Kalman filter.

Related Articles

• Adaptive Filter Theory and Applications in Digital Signal Processing - The theoretical starting point that derives this article’s Wiener-Hopf equation (Eq.(6)) as an MMSE optimization problem from first principles, and numerically verifies three edge cases: an ill-conditioned \(R\) , the convergence relationship with LMS, and non-stationarity.
• Butterworth Filter Design: Theory and Python Implementation - Manual cutoff frequency vs Wiener’s automatic band design.
• Lowpass Filter Design and Comparison - Frequency-response comparison of moving-average, Butterworth, and Chebyshev — the classical alternative to Wiener’s SNR-driven band design.
• Chebyshev Filter Design: Theory and Python Implementation - Equiripple IIR filter that realizes a sharp manual cutoff — instructive contrast against Wiener’s “let the SNR shape the passband” philosophy.
• FIR vs IIR Filters: Characteristics, Design, and Python Implementation - Understand the design philosophy difference between classical filters and Wiener filter.
• Fast Fourier Transform (FFT): Theory and Python Implementation - FFT is essential for implementing the Wiener filter.
• Kalman Filter: Theory and Python Implementation - The Kalman filter as a non-stationary, sequential extension of the Wiener filter.
• Time-Series Anomaly Detection: From Statistical Methods to Kalman Filter in Python - Anomaly detection methods as an application of the Wiener filter.
• Frequency Characteristics of the EMA Filter - EMA is a simple low-pass filter, while the Wiener filter optimally designs the passband based on SNR.
• Moving Average Filters Compared: SMA, WMA, and EMA in Python - Understand the design philosophy difference between moving average and Wiener filter.
• Adaptive Filters (LMS/RLS): Theory and Python Implementation - The LMS algorithm iteratively approximates the Wiener filter’s optimal solution, and RLS converges even faster.
• Window Functions and Power Spectral Density (PSD): Theory and Python Implementation - Covers the Welch-based PSD estimation and window function choices needed to estimate the signal/noise PSDs that drive Wiener filter design.
• Complementary Filter: Theory and Python Implementation - Fusing high-frequency and low-frequency sensors via a complementary filter can be viewed as an intuitive special case of the Wiener filter built on SNR’s frequency dependence.
• LMS and NLMS Adaptive Algorithms: Theory and Python Implementation - Iteratively approximates the Wiener solution \(\mathbf{R}^{-1}\mathbf{p}\) from this article; provides a clean batch-vs-online contrast.
• RLS (Recursive Least Squares) Algorithm: Theory and Python Implementation - An adaptive algorithm that solves this article’s normal equation exactly, recursively, via the matrix inversion lemma. Confirms faster quadratic convergence than LMS and a numerically exact match with the Kalman filter.
• Short-Time Fourier Transform (STFT): Theory and Python Implementation - Foundation for applying Wiener to non-stationary signals: estimate PSD per frame and construct a time-varying Wiener gain.
• Filter Design Selection Guide: 3 axes and characteristic comparison — Use this hub to position the Wiener filter against alternatives during selection.

References

• Haykin, S. (2014). Adaptive Filter Theory (5th ed.). Pearson. Chapters 2-3.
• Wiener, N. (1949). Extrapolation, Interpolation, and Smoothing of Stationary Time Series. MIT Press.
• Vaseghi, S. V. (2008). Advanced Digital Signal Processing and Noise Reduction (4th ed.). Wiley.
• Kailath, T. (1974). A View of Three Decades of Linear Filtering Theory. IEEE Transactions on Information Theory, 20(2), 146-181.