Frequency-domain deconvolution

Question

I want to estimate the impulse response of an LTI system using frequency-domain deconvolution.

Suppose that the input signal $x[n]$ of length $N$, and the unknown system impulse repsonse $h[n]$ of length $L$, I have the output signal $y[n]$ of length $M=N+L-1$.

$$ y[n] = x[n] * h[n] $$

Now convert them to the frequency domain by zero-padding and doing $M$-point FFTs

$$ Y[k] = X[k] \cdot H[k] $$

The unknown system response $H[k]$ can be calculated by frequency-domain deconvolution:

$$ H[k] = \frac{Y[k]}{X[k]} $$

My question is that when taking IFFT to convert $H[k]$ back to the time domain, the IFFT point should be $M$, which gives an $M$-point time-domain sequence $h[n]$. But we know that the length of $h[n]$ should be $L<M$, which implies the last $M-L$ points of $h[n]$ are zero. However in most situations with real world noisy data, this is not guaranteed. What should I do after taking IFFT, force the last samples to be zero?

Answer 1

This depends mostly on the requirement of your specific application. If your signals are noisy, the accuracy with which you can determine the impulse response is inherently limited by the signal to noise ratio. You need to set an accuracy target which needs to be less than the SNR.

For example: Once you have a target, you can simply start truncating and calculate the relative error of the truncation as

$$E(L) = \frac{(\sum_{k=0}^{L}h[k]x[M-k])^2}{\sum y^2[n]}$$

and determine the required impulse response length, $L$ by setting a threshold for the error.

The most suitable error criteria depends on what you want to with the impulse response. You can also evaluate the truncation error over a certain frequency band or with a specific frequency weighting, etc.

Please note that it's important to choose a suitable excitation spectrum. Assuming that the noise $N[k]$ is additive after the system, we can write the transfer function as

$$\hat{H}[k] = \frac{Y[k]+N[k]}{X[k]} = \frac{Y[k]}{X[k]} + \frac{N[k]}{X[k]} = H[k] + \frac{N[k]}{X[k]} $$

So the error to the transfer function is the noise weighted with the inverse of the input spectrum. That can be a big problem if the input signal is bandlimited or has a very different spectral shape than the noise.

Answer 2

One way to mitigate the effect of noise (assuming, like Hilmar, on the output) is to use transfer function estimators. A straight-forward one is: $$\tilde{H}(\omega) = \frac{P_{yx}(\omega)}{P_{x}(\omega)}$$ where $P_{yx}(\omega)$ is the Cross Power Spectral Density of $x$ and $y$ and $P_{x}(\omega)$ is the Power Spectral Density of $x$. These can be computed in many different ways, such as the Welch method.

---

Here's why this estimator is better suited for noisy measurements:

In Frequency domain, you have: $$Y = X\cdot H$$ so you'd be tempted to simply re-arrange to get $H = \dfrac{Y}{X}$
As Hilmar pointed out though, if there's noise $N_0$ in the output, the resulting estimate is: $$\tilde{H} = H + \frac{N_0}{X} = H+\frac{1}{\texttt{SNR}}$$ You can see the problem here, the lower the $\texttt{SNR}$, the higher your estimator's bias.

Here comes the frequency averaging estimator: Let's look at the definition of the Cross Power Spectral Density: \begin{align} P_{yx} &= \mathbb{E}\left[ \overline{X}\cdot(Y+N_0) \right]\\ &= \mathbb{E}\left[ \overline{X}\cdot Y + \overline{X}\cdot N_0\right]\\ &= \mathbb{E}\left[ \overline{X}\cdot HX + \overline{X}\cdot N_0\right]\\ &= H\cdot\mathbb{E}\left[\overline{X}X\right] + \mathbb{E}\left[\overline{X}\cdot N_0\right]\\ &= HP_{x} + \mathbb{E}\left[\overline{X}\cdot N_0\right] \end{align} Assuming the noise in the output, $N_0$, is un-correlated with the input, we have $$\mathbb{E}\left[\overline{X}\cdot N_0\right] = 0$$ $$P_{yx} = HP_{x} \implies H(\omega) = \tilde{H}(\omega) = \frac{P_{yx}(\omega)}{P_{x}(\omega)}$$

Of course this is just theory, since the Expected Value Operator $\mathbb{E}$ is a statistical operator based in infinite time. In a practical measurement scenario (finite time, noisy), the cross term $\overline{X}\cdot N_0$ is never actually $0$, but by using frequency averaging (Welch's method) to compute $P_{xy}$, that cross term does go towards 0 since $X$ and $N_0$ are un-correlated, resulting in a better estimate of $H(\omega)$ than simply dividing the output spectrum by the input spectrum.