Why Wiener filter is not optimal?

I am investigating Wiener filter in the frequency domain, given by $$W(\omega) = \frac{S_{ud}(\omega)}{S_{uu}(\omega)}$$ where $$S_{xy}(\omega)$$ is the cross-power spectral density (CPSD), see scipy.signal.csd or MATLAB cpsd. I used the Python's scipy.signal.csd for estimating $$S_{xy}$$.

Here $$d[n]$$ is the desired signal ("mic") and $$u[n]$$ ("ref") is the reference signal.

Now, here I describe how I generate the signal (in Python). I generate a 10 seconds white noise, and also generate a FIR transfer function true_h (there are 2 ways to generate it: one is designing a FIR bandpass filter and using it, the second is measured from a real room experiment).

from scipy import signal

# Generate filter
true_L = 150 # filter size for simulation
true_h = signal.firwin(true_L, cutoff=[1000, 3000], pass_zero="bandpass", fs=FS)
fft_h = np.fft.rfft(true_h)

# Generate 10 seconds of white noise
np.random.seed(0)
ref = np.random.randn(10*FS)/100
mic = signal.lfilter(true_h, 1, ref)

The goal is to estimate $$w = \mathcal{F}^{-1}(W)$$ where $$W$$ is the estimated Wiener filter in frequency and $$\mathcal{F}^{-1}$$ is the inverse discrete Fourier transform.

The Wiener in the frequency domain is calculated as:

from scipy import signal
# ref = reference signal
# mic = desired signal
freqs, Pxy = signal.csd(ref, mic, fs=FS, nperseg=nfft)
freqs, Pxx = signal.csd(ref, ref, fs=FS, nperseg=nfft)

wiener = Pxy / Pxx

The application of the estimated filter $$w$$ is done as follows:

# Inverse FFT for Wiener filter
wiener_taps = np.fft.irfft(wiener)

# Apply Wiener
error_wiener = mic - signal.lfilter(wiener_taps, 1, ref)

The desired outcome is that $$e[n] = d[n] - w[n] \ast u[n]$$ will be almost $$0$$. But this is not the case. I get a residual which is only around -10 to -20 dB.

When I use frequency/fast block LMS (FBLMS or FLMS), I obtain a residual $$e[n]$$ which is about -40 dB, i.e. the estimation of true_h is almost perfect.

I also plot both Wiener and FLMS estimates of true_h compared to the true transfer function (both in the taps and frequency domain) and get a very good match, but the match is better in FLMS than in Wiener).

This is odd, since Wiener is considered as the optimal filter in the MSE (mean square error) sense. But in this example, Wiener is performing worse than FLMS. See the results below:

The results: Fig 1: Echo Return Loss Enhancement (ERLE). Black = original, Green = after Wiener, Blue = after LMS. Fig 2: residue (e = d - w * u). Green: Wiener filter, Red: FLMS filter. Blue: adaptive FLMS (irrelevant to this question).

Added on edit: Here I further investigate Wiener. I generated a 150-taps FIR, and apply it on a white noise signal (generated by normal distrubtion):

# Generate filter
true_L = 150 # filter size for simulation
true_h = signal.firwin(true_L, cutoff=[1000, 3000], pass_zero="bandpass", fs=FS)
# Generate 10 seconds of white noise
np.random.seed(0)
ref = np.random.randn(10*FS)/100
mic = signal.lfilter(true_h, 1, ref)

Then calculated Wiener by:

1. The default CSD parameters of scipy.
f1, Pxy = signal.csd(ref, mic, fs=FS)
f1, Pxx = signal.csd(ref, ref, fs=FS)
wiener_freq = Pxy/Pxx Fig 3: Misalignment and residual after Wiener (default framing)

1. Taking one long FFT on the entire signal (one "frame"):
f1, Pxy = signal.csd(ref, mic, fs=FS, nperseg=len(ref))
f1, Pxx = signal.csd(ref, ref, fs=FS, nperseg=len(ref))
wiener_freq = Pxy/Pxx Fig 4: Misalignment and residual after Wiener (one big frame)

As we can see, the second method yielded far better result, but still not a perfect 0 residual (i.e. residual < -100 dB).

Added in Edit: It seem that this is the case when the signals are synthetic: i.e. the reference is a white noise generated by numpy.random.randn and the transfer function is generated by scipy.signal.firwin with about 50-250 taps. For real data, it seems that Wiener performs better.

Added in another edit: Figure 5: Wiener (orange) vs. FBLMS (green) for synthetic white noise signal and FIR transfer function. Figure 6: Wiener (orange) vs. FBLMS (green) for real signal (recorded using an audio device)

As it can be seen, for the synthetic case FBLMS is better than Wiener while for the real case it is vice versa, Wiener is better and FBLMS converges to Wiener.

References:

• Can you include code for how you're calculating wiener ?
– Peter K.
Sep 26 '21 at 20:22
• "But in this example, Wiener is performing worse" No. your attempt at a Wiener filter is performing worse. The Wiener filter is optimal in the MSE sense assuming that you have constructed the filter correctly. This means knowing $S_{ud}$ and $S_{uu}$ exactly, for starters, which you don't -- you're starting by determining $S_{ud}$ experimentally. I suspect, but don't know, that even if you're doing the Wiener calculations correctly, the FLMS algorithm is doing a better job because it doesn't have to calculate $S_{ud}$ explicitly (also assuming I remember what the FLMS does). Sep 27 '21 at 4:16
• I added the code in which I calculate the Wiener filter in the frequency domain. Sep 29 '21 at 7:31
• Another note: there is an improvement in the performance of Wiener as I increase the number of taps or FFT length. For example: if the signal is $N_{TOT}$ samples long I take a single FFT of $N_{TOT}$-points instead of the usual overlap-add of the CSD function in Python (that is: framing and windowing the entire signal to small frame, each is $N_{FFT}$ samples long). Sep 30 '21 at 13:12

• You are probably confusing with the noise reduction Wiener filter, used in image processing ($H(k) = S_{xx}(k)/(S_{xx}(k) + S_{nn}(k))$) which its goal is to estimate the target signal from a noisy image. I used the more general Wiener filter, which its goal is to estimate the transfer function from two signals: reference $u$ and desired $d$ (to compare against). The problem is to find $\hat{h}$ such that $y = \hat{h} \ast u$ where the real process is $d = h \ast u$ and $h$ is unknown. Oct 4 '21 at 8:02