Relationship between discrete deconvolution and Toeplitz matrices

Question

I have 2 vectors, $a$ & $c$, both of length $M$. I know they are related by $a*b=c$. My goal is to recover $b$. Obviously $b = \mbox{deconv} (c,a)$. I am only interested in the first $M$ elements of the convolution (i.e., discarding the last $M-1$ elements). I have two methods of doing this deconvolution in MATLAB and I get different answers. I want to know why:

The first method would be using MATLAB's default algorithm:

b=deconv([c;zeros(M-1,1)],a);

The second method would be using Toeplitz matrices:

A=toeplitz(a,[a(1);zeros(M-1,1)]);     
b=A^-1*c

where A is the Toeplitz matrix form of the vector a with 0s on upper triangle.

I believe in theory these methods should give me identical answers, but the answers are somewhat different (but still close). Why is this?

Answer 1

I think you are constructing the convolution matrix incorrectly. Let's do an easy example with some numbers to see how it is constructed

Let's imagine that youhave the following vectors of dimensions $M$ and $N$, And let's for the sake of it, say that $M = 5$ and $N = 4$

$\mathbf{a} = [a_{1},a_{2},a_{3},a_{4}]^{T}$

$\mathbf{b} = [b_{1},b_{2},b_{3},b_{4},b_{5}]^{T}$

$ \mathbf{c} = \mathbf{a}*\mathbf{b} $

In matrix form, if we assume $\mathbf{a}$ is the system response, the convolution is given by

\begin{equation} c = \mathbf{A}\mathbf{b} \end{equation}

with

\begin{equation} \mathbf{A} = \left[ \begin{array}{ccccc} a_{1} & 0 & 0 & 0 & 0 \\ a_{2} & a_{1} & 0 & 0 & 0 \\ a_{3} & a_{2} & a_{1} & 0 & 0 \\ a_{4} & a_{3} & a_{2} & a_{1} & 0 \\ 0 & a_{4} & a_{3} & a_{2} & a_{1} \\ 0 & 0 & a_{4} & a_{3} & a_{2} \\ 0 & 0 & 0 & a_{4} & a_{3} \\ 0 & 0 & 0 & 0 & a_{4} \\ \end{array}\right] \end{equation}

By doing the convolution with this matrix, you get the full operation including the transients. So the correct matlab code would be

A = toeplitz([a,zeros(1,M-1)],[a(1),zeros(1,M-1)])

(Check help toepliz to verify how it constructs the matrix)

Use this matrix and the result of the convolution is exactly as with the matrix multiplication.

Now, a word of wisdom, as you can clearly see, the full convolution matrix with the transients is not square, so in order to "invert" it and do deconvolution you will need to use the Moore-Penrose Pseudoinverse. And you will get some differences in the beginning and end of the resulting vector, but it is the best you can do, as using the pseudoinverse would be the same as performing least squares, since the deconvolution is given by

$\mathbf{b} = \left(\mathbf{A}^{T}\mathbf{A}\right)^{-1}\mathbf{A}^{T}\mathbf{c}$

This result will always hold, because the dimensions of $A$ are always $(M+N-1)\times M$, so the matrix is always "skinny" and the equation holds.

Try this, it should give you the result you were looking for.