Deconvolution with noisy measurement of impulse response function

Question

I observe a noisy complex-valued signal that has passed through some linear time-invariant filter: $$y(t) = (h * x)(t) + n(t)$$ where $h(t)$ is a causal and finite impulse response, and $n(t)$ is additive Gaussian white noise, $n \sim \mathcal{CN}(0, \sigma^2)$.

I also have a noisy approximation of the filter: $\hat{h}(t) = h(t) + \epsilon$. You can assume $\epsilon$ is also AGWN.

The goal is to recover $x(t)$ as best as possible given only $y(t)$ and $\hat{h}(t)$.

What is the best deconvolution scheme to use in this context? Additionally, how would you avoid overfitting (i.e., is there a nice way to regularize the solution for cases when the signal-to-noise ratio is quite low)?

More deconvolution literature talks specifically about 2D images or stochastic processes. I can't find much literature on deconvolution in a signal-processing context.

---

To add some context to answer some questions in the comments:

I encountered this problem in the context of radio astronomy. I have a radio source that emits from signal $x(t)$. As the signal propagates through the interstellar medium, it can take many different paths as it propagates towards us - this is the filter, $h(t)$. The noise $n(t)$ primarily comes from terrestrial sources and the receiver at the radio telescope (hence, why no noise actually goes through the filter).

I only have noisy approximations of the filter because the radio source I am observing happens to randomly emit incredibly bright "impulses" (few nanoseconds wide) every now and then, which lets me directly but noisily measure $h(t)$ (literally, the impulse response).

In practice, the length of $h[n]$ is $\sim\!10000$, and the length of $x[n]$ is $\sim\!30,000,000$.

Answer 1

Since we're dealing AWGN noise, we can try formulate the problem as some kind of a Least Squares problem.

The way to model the problem is as we model the Total Least Squares Problem:

$$ \boldsymbol{y} \approx H \boldsymbol{x} \implies \boldsymbol{y} = \left( H + E \right) \boldsymbol{x} + \boldsymbol{n} $$

Where in the model $E$ and $n$ are composed by AWGN noise.
Since we're dealing with zeros mean noise which is Gaussian, minimizing their variance is equivalent to the Maximum Likelihood estimator, hence we can solve the following problem:

$$ \begin{align} \arg \min_{\boldsymbol{x}, E, \boldsymbol{n}} \quad & {\left\| \begin{bmatrix} E & \boldsymbol{n} \end{bmatrix} \right\|}_{F}^{2} \\ \text{subject to} \quad & \begin{aligned} \boldsymbol{y} & = \left( H + E \right) \boldsymbol{x} - \boldsymbol{n} \end{aligned} \end{align} $$

The solution is given by $ \boldsymbol{x} = -{V}_{H \boldsymbol{y}} {V}_{\boldsymbol{y} \boldsymbol{y}}^{-1} $. Pay attention that the negation of $ \boldsymbol{n} $ has no effect, I just wanted to match the derivation of Wikipedia. So if you want to go deeper into the derivation, things will match.

Where we define the matrices by the SVD:

$$ \begin{bmatrix} H & \boldsymbol{y} \end{bmatrix} = \begin{bmatrix} {U}_{H} & {U}_{\boldsymbol{y}} \end{bmatrix} \begin{bmatrix} {\Sigma}_{H} & 0 \\ 0 & {\Sigma}_{\boldsymbol{y}} \end{bmatrix} {\begin{bmatrix} {V}_{H H} & {V}_{H \boldsymbol{y}} \\ {V}_{\boldsymbol{y} H} & {V}_{\boldsymbol{y} \boldsymbol{y}} \end{bmatrix}}^{*} $$

In the above we use $ H $ is the matrix formulation of the convolution operator. You may build it according to Generate the Matrix Form of 1D Convolution Kernel.

So the recipe is quite simple in your case.

How to Solve

Given we know how to generate $H$ matrix (See Generate the Matrix Form of 1D Convolution Kernel), the solution for the case above, where $ \boldsymbol{y} $ is a vector, is pretty straight forward:

function [ vX ] = TlsRegression( mH, vY )
numCols = size(mH, 2);
[~, ~, mV] = svd([mH, vY]);
vX = -mV(1:numCols, numCols + 1) / mV(numCols + 1, numCols + 1);
end

For instance, for Linear Regression of a quadratic polynomial we can see a comparison between the Least Squares and the Total Least Squares solution:

The difference is about the error to be minimized:

The Total Least Squares (TLS) minimizes the perpendicular error, the Least Squares (LS) minimizes the vertical distance, namely it assumes the model is perfect.

Answer 2

I used $h[n]$, $\hat{h}[n]$, $x[n]$ in my answer below to represent samples of $h(t)$, $\hat{h}[t]$, $x[t]$ since I am proposing a discrete time solution using FFTs.

Proposed Approach

If $h[n]$ is significantly shorter than $x[n]$ and $y[n]$, (as would be typically expected) my solution would be to zero pad $\hat{h}[n]$ and $y[n]$ to the same length and at least double the length of the longest of the two (to avoid time domain aliasing), and then use the FFT of each as $\hat{H}[k]$ and $Y[k]$. $x[n]$ is then estimated by taking the inverse FFT of the ratio as:

$$x[n] = \text{ifft}(Y[k] / \hat{H}[k])$$

This approach is effective as long as there are not deep nulls in the channel within the occupied bandwidth of the signal (both of which can be confirmed via the FFT's performed). Deep nulls in the channel can occur when the delay spread (length of the dominant portions of the channel impulse response) significantly exceeds the inverse of the signal bandwidth. Such nulls will lead to noise enhancement degrading the result if not compensated for. By deep nulls, this would typically be channel nulls of -25 dB or more. I show an example simulation result at the end of this post with an in-band channel null of -10 dB, which as demonstrated is not an issue. This is a motivation for decision feedback equalization in communication systems, but for this application, the nulls can be compensated for by limiting the minimum value in $\hat{H}(k)$ for the portions of the channel response that are within the occupied bandwidth of the signal (I suggest -25 dB, but confirm in simulation what would be best based on the characteristics of the waveform and channel). Further, if the energy for $h[n]$ is known to be in a constrained portion of the spectrum, further benefit can be reached by averaging or filtering $\hat{h}[n]$ to reduce the noise improving the estimate for $\hat{h}[n]$ to be used in the processing above.

As Royi has clarified in the comments, the above approach is a least squares solution, in a very straightforward and efficient computation!

A simulation confirming the proposed approach is shown in the plot below, where $x[n]$ was estimated from $\hat{h}[n]$ and $y[n]$ alone:

Answer 3

One practical method would be to use the matched filter. Matched filters are commonly thought of as useful for finding the scale and delay of a known signal buried in noise, but if you think about it, the same principle allows the matched filter to perform channel estimation (which is essentially what your problem is).

This is why matched filters are used in remote sensing/image formation. For example, synthetic aperture radar attempts to form an image of an unknown scene by convolving a known "training sequence" (i.e., the transmitted radar waveform, $h(t)$ in your case) with the unknown scene/channel impulse response ($x(t)$ in your case). Of course, $x(t)$ has some noise added to it, and perhaps you don't know the waveform the radar actually transmitted due to hardware distortion, so you have a noisy estimate of the waveform, $\hat{h}(t)$, rather than $h(t)$.

This problem is solved by cross-correlating the received signal, $y(t)$, with transmitted signal, $\hat{h}(t)$ (matched filtering). Or equivalently, convolving $y(t)$ with the conjugated and time reversed version of $\hat{h}(t)$, $\hat{h}^*(-t)$. For Gaussian noise, the matched filter maximizes SNR.

Another Method

One could also use regularized deconvolution. In the frequency domain, deconvolution just means dividing the observed spectrum $Y(f)$ by the filter spectrum $H(f)$. But this can be subject to dividing by zero at points, so typically $H(f)$ is regularized to have some moderate minimum value first.

A similar solution is Wiener deconvolution.