The Wiener filter

A slow random signal is buried in white noise, and the noise grows. The best filter passes each frequency in proportion to the signal’s share of the power there. Watch the gain: as the noise rises, the point where it falls to ½ slides down.

The best guess of a signal in noise

In Simple smoothing filters (18.2), “A longer average: less noise, more delay” cleaned up a noisy scale with a moving average. A longer average removed more noise, but it also blurred the parcel’s landing more. So how long should the average be, and should every reading get the same weight?

This page answers both questions for a signal and a noise whose statistics we know. The answer is the Wiener filter, worked out by Norbert Wiener in the 1940s.

Here is the setting. The signal I want is $d_\text{des}[n]$, the desired signal. I never see it alone. I see it with white noise $v[n]$ of power $\sigma_v^2$ added, the white noise of “White noise forgets, coloured noise remembers” in Random processes (24.2):

Both have mean 0 and are stationary, and they are uncorrelated: $\mathbb{E}\{d_\text{des}[n]\,v[m]\}=0$ for every $n$ and $m$. The noise does not know what the signal is doing.

A filter turns $x[n]$ into an output $y[n]$, my guess at the signal. The error is what the guess misses, $e[n]=d_\text{des}[n]-y[n]$. On this page $e[n]$ is this error, not the rounding error of Quantization & noise (11.1).

One error sample can be large or small by chance, so I judge a filter by the error’s average size. The mean-square error is the error’s power:

Squaring makes a miss of 2 count four times as much as a miss of 1, and it leads to equations we can solve. That was the reason for least squares in “Least squares or least worst” of Optimal FIR design (19.3) too. The Wiener filter is the filter with the smallest $J$, and that smallest value is the minimum mean-square error $J_\text{min}$.

Two lazy filters set the bar. Passing everything, $y[n]=x[n]$, leaves the error $-v[n]$, so $J=\sigma_v^2$, the noise’s power. Passing nothing, $y[n]=0$, leaves the error $d_\text{des}[n]$, so $J$ is the signal’s power $\sigma_\text{des}^2$. A good filter beats both.

The best filter, frequency by frequency

First, let the filter use the whole record, the samples after $n$ as well as those before. That is a non-causal filter, fine for cleaning a recording afterwards. Then the best filter can be found one frequency at a time.

A filter with impulse response $h[n]$ and frequency response $H$ gives the output $y=h*d_\text{des}+h*v$. So the signal reaches the error through the system $1-H$, and the noise through $H$ with a minus sign. A sign does not change a power.

In Power spectral density (24.4), “White, pink, brown” ended on the rule $S_y=\lvert H\rvert^2S_x$: a system multiplies a PSD by its gain squared. The two parts of the error are uncorrelated, so their autocorrelations add, and by Wiener–Khinchin so do their PSDs. Using the rule twice, the error’s PSD is

where $S_\text{des}$ and $S_v$ are the PSDs of the signal and of the noise. I leave out the argument $e^{j\Omega}$ to keep it short. The noise is white, so $S_v=\sigma_v^2$, the same at every $\Omega$ (24.4).

$J$ is the error’s power, the area under that PSD divided by $2\pi$ (24.4). Each frequency adds its own share, and the filter can choose its gain at each frequency separately. So I make the error’s PSD as small as possible at every $\Omega$.

At one frequency this is a tug of war. A gain near 1 lets the noise through; a gain near 0 throws the signal away. An imaginary part of $H$ adds to both terms, so the best $H$ is real. For real $H$ the error’s PSD is a parabola in $H$, and its lowest point is where its slope, $-2(1-H)S_\text{des}+2HS_v$, is 0:

In words: at each frequency, the gain is the signal’s share of the power there. Where the signal is ten times the noise, the gain is $10/11=0.91$. Where they are equal it is ½, and where the noise is ten times the signal it is $1/11=0.091$.

Put $H_\text{opt}$ back in, and the error’s PSD becomes $S_\text{des}S_v/(S_\text{des}+S_v)$. That is $H_\text{opt}S_v$, and also $(1-H_\text{opt})S_\text{des}$, so it is smaller than both $S_v$ and $S_\text{des}$. Its area gives the least error:

So the Wiener filter beats both lazy filters at every frequency, and therefore in total. $H_\text{opt}$ is real and positive, so it adds no delay, like the zero-phase response of Linear-phase systems (17.3). It is also even in $\Omega$, so its impulse response is symmetric about $n=0$: it leans on the samples after $n$ as much as on those before. That is why it needs the whole record.

One signal for the whole page

The signal on this page is the AR(1) noise of 24.2 with $a=0.9$. Each sample is 0.9 times the one before plus fresh white noise of power 0.19, so its power is $\sigma_\text{des}^2=1$ and its autocorrelation is $0.9^{\lvert\ell\rvert}$. In “Memory in time, narrowness in frequency”, 24.4 found its PSD:

It is 12.8 dB at $\Omega=0$ and $-12.8$ dB at $\pi$, and its area divided by $2\pi$, the power, is 1.000. It is a slow signal, with its power gathered at low frequencies. The white noise spreads its power evenly, so the Wiener filter will be a low-pass filter, and the noise level decides how far up it reaches.

Pass what is mostly signal

The picture at the top of the page draws the two PSDs and the gain $H_\text{opt}$ while the noise grows. Think of listening to a friend in a noisy room: you attend to the pitch range where their voice is louder than the room.

Its noise power is in dB, $10\log_{10}\sigma_v^2$, so −10 dB is a noise power of 0.1 and 10 dB is 10. Its error power is $J_\text{min}$.

Notice where the gain is ½. As the noise rises a hundredfold, from −10 dB to 10 dB, that point moves from 0.516π down to 0.032π. More noise pulls the gain down everywhere, and only the frequencies where the signal is strongest keep a gain above ½.

At 0 dB the gain is ½ at 0.144π. That is where $S_\text{des}=1=S_v$, the frequency where “Memory in time, narrowness in frequency” in 24.4 saw the AR(1) PSD cross the white level, at $\cos\Omega=0.9$.

Now compare the error with the two lazy filters. At −10 dB the error, 0.059, is 59 % of the noise’s 0.1. At 10 dB it is 0.586, under 6 % of the noise’s 10, and also below the signal’s power, 1, which passing nothing would leave.

When the clip has finished, drag the slider to 5 dB: the gain is ½ at 0.075π, and the error power is 0.375. At −20 dB the gain stays above ½ everywhere, and at 20 dB below it.

A causal filter with a few taps

$H_\text{opt}$ needs the future, and its impulse response never quite ends. A filter running live must be causal, and I would like a short one. So let’s look for the best FIR filter of $N_h$ taps $h_0$ to $h_{N_h-1}$:

The error is linear in the taps, so $J$ is a quadratic in them, a bowl, like the mean square of 19.3’s least squares. Its lowest point comes from $N_h$ linear equations.

Parametric models and linear prediction (25.3) showed the quickest way to them. At the best taps, the error must be uncorrelated with every input sample the filter uses. If it were still correlated with ${x[n-i]}$, adding a little more of ${x[n-i]}$ would shrink it. So for $i=0,\dots,N_h-1$,

In calculus terms, this sets the slope of $J$ along each tap to 0: that slope is $-2\,\mathbb{E}\{e[n]\,x[n-i]\}$.

Put in $e[n]=d_\text{des}[n]-\sum_kh_k\,x[n-k]$. Averages add, and each $\mathbb{E}\{x[n-k]\,x[n-i]\}$ is the autocorrelation $R_x[i-k]$ of 24.2. So

The right side is a cross-correlation, as in Correlation (24.3), but averaged over the process: what I want against what I see, $i$ samples earlier. These are the Wiener–Hopf equations.

In matrix form they read $\mathbf{R}\mathbf{h}_\text{opt}=\mathbf{r}$. Here $\mathbf{h}_\text{opt}$ holds the best taps, and $\mathbf{R}$ is the $N_h\times N_h$ matrix with $R_x[i-k]$ in row $i$, column $k$: the Toeplitz matrix of 25.3. The vector $\mathbf{r}$ holds the right sides. In 25.3 it held autocorrelations; here its entries are the cross-correlations $\mathbb{E}\{d_\text{des}[n]\,x[n-k]\}$.

From linear algebra this page needs one fact. When $\mathbf{R}$ is positive definite, as 25.3 defined it, and it is on this page, the system has exactly one solution. SciPy’s solve_toeplitz finds it with a Levinson-type recursion that uses the Toeplitz pattern, as 25.3 did.

For this page’s signal both sides are easy. The signal and the noise are uncorrelated, so their autocorrelations add. The noise’s is $\sigma_v^2$ at lag 0 and 0 elsewhere (24.2), so

On the right, $\mathbb{E}\{d_\text{des}[n]\,x[n-k]\}$ is $\mathbb{E}\{d_\text{des}[n]\,d_\text{des}[n-k]\}$ plus $\mathbb{E}\{d_\text{des}[n]\,v[n-k]\}=0$, so the entries of $\mathbf{r}$ are $0.9^k$.

With noise power 1 and three taps, the Wiener–Hopf equations are

and their solution is $h_0=0.330$, $h_1=0.225$ and $h_2=0.170$.

How small is the error at the best taps? There the error is uncorrelated with every ${x[n-k]}$, and so with $y[n]$, which is made from them. So $J=\mathbb{E}\{e[n]\,d_\text{des}[n]\}$, and multiplying out gives

where $\mathbf{r}^\top\mathbf{h}_\text{opt}$ is the sum of the products of matching entries. For the three taps above it is $1-(0.330+0.9\cdot0.225+0.81\cdot0.170)=0.330$.

The best FIR filter, from correlations

The second instrument solves the Wiener–Hopf equations on the page with noise power 1, and runs the filter over a seeded record. Picture a weighted average of recent readings, whose weights come from how alike neighbouring readings are.

The record is 4000 samples of this page’s signal. I took 4500 draws of the site’s generator (seed 262), scaled them by $\sqrt{0.19}$ and passed them through the AR(1) recursion, then dropped the first 500 outputs. The noise is 4000 more draws (seed 2620). Measured over the record, the signal’s power is 0.875 and the noise’s 0.986.

Watch the solid output follow the dashed signal through the dots, and the error power level off as taps are added. The measured value is the mean of $e^2[n]$ over samples 100 to 3999.

Notice how the error levels off. One, two, four and eight taps give 0.500, 0.373, 0.314 and 0.304, which is −3.01, −4.28, −5.03 and −5.17 dB. From eight taps on, the error stays at 0.304, and even 400 taps reach no lower.

Look at the taps in the end frame. Apart from the last few, each is 0.63 times the one before. The weights fade with distance, like the signal’s own memory, but faster than its 0.9 per step. An older sample brings as much noise as a newer one, but tells less about the present. From the ninth tap on, every tap is below 0.008, so more taps add almost nothing.

There is a neat check in the captions: the first tap equals the error power, 0.500 with 0.500 and 0.314 with 0.314. That holds whenever the noise is white. The error is uncorrelated with $x[n]$, so $J=\mathbb{E}\{e[n]\,d_\text{des}[n]\}=-\mathbb{E}\{e[n]\,v[n]\}$, and only the tap $h_0$ ever meets $v[n]$, giving $J_\text{min}=\sigma_v^2h_0$.

The measured errors sit a little below the theory: 0.464 against 0.500, 0.303 against 0.314 and 0.299 against 0.304. The theory is an average over every possible record, and this one is a single record. One reason is visible in the subtitle: its signal is weaker than nominal, with power 0.875 instead of 1.

When the clip has finished, use the slider to try 2 taps: the error power drops from one tap’s 0.500 to 0.373, and the measured error to 0.351.

Looking ahead, and what it assumes

This FIR filter only looks back. Allowed to look ahead too, the error at noise power 1 falls to the first instrument’s 0.218. That is below 0.304, the least error any causal filter can reach here.

The best causal filter of unlimited length, the causal IIR Wiener filter, comes from splitting $S_x$ into a causal and an anti-causal factor, a step called spectral factorisation. I leave it out here; for this signal its error is the same 0.304 that long FIR filters reach.

The Wiener filter makes three assumptions. The signal and the noise are stationary, they are uncorrelated with each other, and their correlations are known. When the correlations are not known, Adaptive filters: LMS (26.3) learns the filter from the data instead.

The derivation never used the noise being white, only uncorrelated with the signal. With coloured noise, $S_v$ and $R_x$ change, and the method does not. If the signal also passes through a known blur first, the same equations give the filter that undoes the blur as far as the noise allows: Wiener deconvolution.

Worked example

1. The gain by hand. At noise power 1 and $\Omega=0$, $S_\text{des}=0.19/(1.81-1.8)=19$, so $H_\text{opt}=19/(19+1)=0.950$. At $\Omega=\pi$, $S_\text{des}=0.19/3.61=0.0526$, so $H_\text{opt}=0.0526/1.0526=0.050$. These are the gains in the first instrument’s middle caption.

2. One tap by hand. With one tap the Wiener–Hopf equations shrink to $h_0R_x[0]=\mathbb{E}\{d_\text{des}[n]\,x[n]\}$, that is $h_0(1+1)=1$. So $h_0=0.500$, and $J_\text{min}=1-1\cdot0.5=0.500$, half the noise power. On the record the measured error is 0.464.

3. Two taps by hand. The equations are $2h_0+0.9h_1=1$ and $0.9h_0+2h_1=0.9$. Taking 0.45 times the first from the second gives $1.595h_1=0.45$, so $h_1=0.282$ and $h_0=(1-0.9h_1)/2=0.373$. Then $J_\text{min}=1-(0.373+0.9\cdot0.282)=0.373$, which equals $h_0$, as the noise is white.

4. The non-causal limit. At noise power 1 the non-causal filter reaches $J_\text{min}=0.218$, which is −6.6 dB. The best causal filter reaches 0.304, or −5.18 dB.

Where you’ll meet this

Many speech enhancers apply a Wiener gain to each bin of the short-time spectrum of Spectrograms & the STFT (15.5). They measure the noise’s PSD in the pauses between words. Images blurred by a known lens or motion are restored by Wiener deconvolution.

In a receiver, the channel smears each symbol into its neighbours. An equaliser built from these same equations undoes the smear, in Channels and equalisation (33.2). The Wiener filter is also the target that the adaptive filters of Adaptive filters: LMS (26.3) learn. The Kalman filter of RLS and the Kalman filter (26.4) extends it to signals whose statistics change.

Compare it with Matched filters and detection (26.1). The matched filter asks whether a known pulse is there, and its output need not look like the pulse. The Wiener filter wants the waveform itself, as close as it can get.

For more, see N. Wiener, Extrapolation, Interpolation, and Smoothing of Stationary Time Series (1949); S. Haykin, Adaptive Filter Theory (5th ed., 2014), ch. 2; and J. G. Proakis and D. G. Manolakis, Digital Signal Processing, ch. 12. In SciPy, solve_toeplitz(c, r) gives the FIR Wiener taps, with c holding $R_x[0]$ to $R_x[N_h-1]$ and r the cross-correlations.

Reference card