2
$\begingroup$

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)

Fig 1: Echo Return Loss Enhancement (ERLE). Black = original, Green = after Wiener, Blue = after LMS.

Fig 4: residue (e = d - w * u). Green: Wiener filter, Red: FLMS filter. Blue: adaptive FLMS (irrelevant to this question).

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: Wiener 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: 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 Figure 5: Wiener (orange) vs. FBLMS (green) for synthetic white noise signal and FIR transfer function.

Figure 6 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:

$\endgroup$
4
  • $\begingroup$ Can you include code for how you're calculating wiener ? $\endgroup$
    – Peter K.
    Sep 26 at 20:22
  • 2
    $\begingroup$ "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). $\endgroup$
    – TimWescott
    Sep 27 at 4:16
  • $\begingroup$ I added the code in which I calculate the Wiener filter in the frequency domain. $\endgroup$ Sep 29 at 7:31
  • $\begingroup$ 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). $\endgroup$ Sep 30 at 13:12
0
$\begingroup$

The Wiener filter minimizes the MSE under a series of assumptions, which usually don’t match reality. That means that, in practice, it will not give the solution with the smallest MSE possible.

Using the Wiener filter for deconvolution is the simplest thing you can do, and there are many other, more complex deconvolution algorithms that produce better results.

$\endgroup$
4
  • $\begingroup$ But why does it fail here? The signal is a synthetic (Gaussian) white noise and the transfer function is a FIR with 150 taps, and there is no interfering signal or external noise. These are theoretical settings in which Wiener should be optimal. $\endgroup$ Oct 3 at 7:48
  • $\begingroup$ @Triceratops: One of the reasons could be that you swapped mic and ref here: error_wiener = mic - signal.lfilter(wiener_taps, 1, ref). But your whole Wiener implementation is odd, as you leave out the noise power spectral density altogether, and I miss a 1/H term (with 0 noise, as is your case, the filter in frequency domain should be equal to 1/fft_h). $\endgroup$ Oct 3 at 14:25
  • $\begingroup$ 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. $\endgroup$ Oct 4 at 8:02
  • $\begingroup$ @Triceratops Yes, I know quite well what Wiener deconvolution is. The goal is to recover the original signal from the convolved one knowing the impulse response function. If you want, you can swap the impulse response and the original signal in those equations, thereby estimating the impulse response from the original and convolved signals. But that is not what you are doing, and it is not how you are testing the output either. You are actually applying the estimated Wiener filter to the original signal, and comparing it to the convolved signal. $\endgroup$ Oct 4 at 11:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.