Estimate the Filter Coefficients of 1D Filtration (Convolution)

Question

I have an output signal $y$ which is an input signal $x$ convolved $\star$ with an impulse response function $h$ with some added noise $n$ :

$$y(t) = h(t) \star x(t) + n(t)$$

I know the input signal $x$ and output signal $y$ and would like to calculate $h$ the impulse response function. I found that deconvolution is not as straightforward as convolution because the input signal contains zeros and then division in the frequency domain would not be defined. Looking around the internet for ways to "deconvolve" if found two methods: Wiener deconvolution and regularized deconvolution. The Wiener deconvolution seemed easier to understand so I wanted to try and implement it in Matlab (the Matlab function deconv gives me errors about the input signal having a zero at the first entry and if I read the help file it only seems to work correctly for polynomials?).

So per the Wikipedia - Wiener Deconvolution explanation you want to find $g$ so that:

$$\hat{x}(t) = g(t) \star y(t)$$

But then in the definition of how to calculate $G$ they use all variables in the original equation. Also they only show how to find $x$ but probably $x$ and $h$ can be exchanged because convolution is commutative, but I'm unsure about the correct length of both vectors. Currently they are the same length.

$$G(f) = \frac{H^*(f) S(f)}{|H(f)|^2 S(f) + N(f)}$$

where:

• $H = {\tt fft}(h)$
• $G = {\tt fft}(g)$
• $S =$ power spectral density of $x$ ? is this ${\tt fft}(x)$?
• $N = $mean power spectral density of $n$, don't really understand what this is

My question is how do I get the impulse response function $h$ without already knowing it (as it is both in the definition of $G$ and in the original equation)? Since I know both input and output it should not be very different from finding the input with a known impulse response function. I want to do this in Matlab.

Answer 1

This is a nice question.
I will try solving it using 2 approaches (Which are basically the same).

The solution is the Least Squares Solution:

$$ \hat{h} = \arg \min_{h} \frac{1}{2} \left\| h \ast x - y \right\|_{2}^{2} $$

We assume data is given in finite discrete form (As it is in practice).
The convolution is done in valid mode (Like in MATLAB valid property).

Direct Solution

One could write the above as:

$$ \hat{h} = \arg \min_{h} \frac{1}{2} \left\| X h - y \right\|_{2}^{2} $$

Basically using the commutative property of the convolution one could the above as $ x $ was the filter and hence build a Convolution Matrix from $ x $.

Now, the solution is the usual Least Squares:

$$ \hat{h} = \arg \min_{h} \frac{1}{2} \left\| X h - y \right\|_{2}^{2} = \left( {X}^{T} X \right)^{-1} {X}^{T} y $$

Pay attention that for long signals the matrix is huge.
Hence this method is practical only for small signals.

Iterative Solution

Well, one could use the Gradient Descent Method to minimize the cost function:

$$ f \left( h \right) = \frac{1}{2} \left\| h \ast x - y \right\|_{2}^{2} $$

The Gradient is given by (Correlation is the Adjoint Operator of Convolution):

$$ \frac{d}{d h} f \left( h \right) = \left( h \ast x - y \right) \star h $$

Where $ \star $ stands for Cross Correlation.

In MATLAB Code it is given by:

hObjFun = @(vH) 0.5 * mean( (conv2(vX, vH, 'valid') - vY) .^ 2 );
vG      = conv2(flip(vX, 1), (conv2(vX, vHEst, 'valid') - vY), 'valid');

Where all data is in Column Vector (The conv2() is used as it is faster since conv() is just a wrapper around it).

This method is fast and not limited by the size (Or less, to be accurate).

The full MATLAB Code is in my StackExchange Signal Processing Q18993 GitHub Repository (Look at the SignalProcessing\Q18993 folder).