4
$\begingroup$

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?

$\endgroup$

1 Answer 1

4
$\begingroup$

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)))
$\endgroup$
7
  • $\begingroup$ See more: dsp.stackexchange.com/questions/71811/… $\endgroup$
    – ZR Han
    Apr 21, 2021 at 2:28
  • $\begingroup$ Thanks ZR: I do understand the theory reasonably well. However in practice for real data I dont see any difference between the two forms of spectral division, (unless I need to add stabilisation to the denominators). $\endgroup$ Apr 21, 2021 at 2:42
  • $\begingroup$ Thanks again ZR. Does the theory of estimators cover the cases where there is noise on both input and output. If the input and output noises are correlated would H1 and H2 give similar estimates. $\endgroup$ Apr 26, 2021 at 1:13
  • 1
    $\begingroup$ @telemeister If there is noise on both input and output, you can use H3 estimator which equals to the average of H1 and H2. But all of them cannot handle correlated noise. $\endgroup$
    – ZR Han
    Apr 26, 2021 at 1:31
  • $\begingroup$ Still not really able to demonstrate significant differences between the different estimators using real data. I don't have matlab but I'd like to simulate your code in Python. I can google what the matlab functions do. Can you tell me what is the form of the input signal "mls"? $\endgroup$ Apr 30, 2021 at 9:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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