Frequency-domain deconvolution: "Direct" filtering vs "Wiener" filtering

Question

Can someone help with clarifying the difference between two approaches to frequency domain "deconvolution:

For the frequency domain problem:
We want to find a filter $F(\omega)$ which will convert some observed signal $B(\omega)$ into a desired signal $D(\omega)$, which is known i.e. $D(\omega) = B(\omega)F(\omega)$

1. The intuitive / direct approach is to just calculate $F(\omega)=D(\omega)/B(\omega)$

2. Application of the least-squares criterion leads to $F(\omega)=R_{DB}(\omega)/R_{BB}(\omega)$ i.e. $D(\omega)B^*(\omega)/B(\omega)B^*(\omega)$

Q1: If the direct division can be done successfully (i.e. no zeros in the denominator), the two approaches seem to always give identical answers in practice. Is that because the direct division is actually handled in code (e.g. Python) using the conjugates.

Q2: Does the distinction between the two approaches only arise when the division requires the addition of some stabilizing noise in the denominator?

Answer 1

The direct approach is noise sensitive and the second one (which is also known as $H_1$ estimator) is somehow noise resistant.

$H_1$ estimator assumes that there is no noise at the input signal, $B(\omega)$ in this case, and the noise exists only at output $D(\omega)$ and is uncorrelated to input $B(\omega)$. Therefore, the observed output signal is

$$ D'(\omega) = D(\omega)+N(\omega) $$

where $N(\omega)$ is the uncorrelated noise. Under this assumption, the $H_1$ estimator gives a better result

$$ F'(\omega) = \frac{D'(\omega)B^*(\omega)}{B(\omega)B^*(\omega)} = \frac{D(\omega)B^*(\omega)+N(\omega)B^*(\omega)}{B(\omega)B^*(\omega)} = F(\omega) + \frac{N(\omega)B^*(\omega)}{B(\omega)B^*(\omega)} $$

Since the noise is uncorrelated to the input signal, the second term can be eliminated by time averaging.

---

Here is some code and you can try some other SNR values. In practice for real data, you cannot obtain the clean signal without noise. All signal you measure is kind of noisy.

clear

% Input signal
x = mls;

% Comb filter (unknown system)
b = [1, zeros(1, 50), 1];
y = filter(b, 1, x);
fvtool(b, 1)

% Add white noise to the output
SNR = 10; % try some lower value and see what's happening
yn = awgn(y, SNR);

% H1 estimator
pxy = cpsd(x, yn);
pxx = cpsd(x, x);
H1 = pxy./pxx;
figure; subplot(211); plot(mag2db(abs(H1)))

% direct division
X = fft(x);
Y = fft(yn);
X = X(1:ceil(length(X)/2));
Y = Y(1:ceil(length(Y)/2));
H0 = Y./X;
subplot(212); plot(mag2db(abs(H0)))