Can anyone help me in understanding of deconvolution based on toeplitz matrix?

First of all, d is a trace whose size is (1,2500), p is a trace whose size is (1,2500).

For d*w=p, which is Dw=p, where D is Toeplitz matrix made by d. Take d and p as an example [1 2 3 4 5].

I construct D that the first column is [1 2 3 4 5 0 0 0 0], the first row is [1 0 0 0 0], which means d lags 2 time step up and 2 time step down, the 0 time lag is place 3 when the first place is 1.

To find w, I use deconvolution, therefore, w=((D^T)D)^(-1)(D^T)p

Confusion comes in here. the number of row of D is 9, however, the p is (1,5) size, how should I concatenate p using 0? Just under the last element of p or two 0 before and two 0 behind?

What I really want to know is the exact form of deconvoution in programs. Really Really appraciate it if anyone could help. Please!

Answer 1

[William Rose]

@Junjie Zhao,

Here is an examle, using d=[1 2 3 4 5], which you suggested. Since d has length 5, then the convolution d*w = p (where * indicates convolution) will be 4 elements longer than w. So if w has length 5, p will have length 9, etc. For this example, I will assume w has length 6, in order to demonstrate that the length of w does not have to equal the length of d. I will do the forward convoltion (compute p=d*w) first. Then I will do the inverse convolution.

Compute Toeplitz matrix:

w=[-1 0 1 2 1 3]';
d=[1 2 3 4 5]';
%nr=length(d)+length(w)-1;  % rows in Toeplitz matrix
%nc=length(w);  % columns in Toeplitz matrix
c=[d;zeros(length(w)-1,1)];  % column 1 of Toeplitz matrix
r=[1,zeros(1,length(w)-1)]; % row 1 of Toeplitz matrix
D=toeplitz(c,r);
disp(D)

     1     0     0     0     0     0
     2     1     0     0     0     0
     3     2     1     0     0     0
     4     3     2     1     0     0
     5     4     3     2     1     0
     0     5     4     3     2     1
     0     0     5     4     3     2
     0     0     0     5     4     3
     0     0     0     0     5     4
     0     0     0     0     0     5

Use D to compute p=convolution of d with w:

p=D*w;
disp(p')

    -1    -2    -2     0     3    15    22    23    17    15

Compute estimate of w, from p, using the Toeplitz matrix:

wEst=inv(D'*D)*D'*p;
disp(wEst')

   -1.0000   -0.0000    1.0000    2.0000    1.0000    3.0000

The result shows that the estimate of w equals the original w.

[Paul]

Hi Wiliam,

Can also use convmtx for D and save a few lines of code.

w=[-1 0 1 2 1 3]';
d=[1 2 3 4 5]';
D = convmtx(d,numel(w))

D = 10×6

1 0 0 0 0 0 2 1 0 0 0 0 3 2 1 0 0 0 4 3 2 1 0 0 5 4 3 2 1 0 0 5 4 3 2 1 0 0 5 4 3 2 0 0 0 5 4 3 0 0 0 0 5 4 0 0 0 0 0 5

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[William Rose]

@Paul, thank you. I didn't know about convmtx().

[Junjie Zhao]

Sorry for not replying your answer timely. Thank you so much for your effort. However, what I want to know is how should I process p in the deconvolution operation to achieve the effect of deconvolution. If you could help me with this, I will be very grateful.

[William Rose]

@Junjie Zhao,

In the discussion above, p is the convolution of d and w. If you know p and d, you can estimate w by deconvolution, as follows (and as described above):

First, let's create p, by convolution:

w=[-1 0 1 2 1 3]';
d=[1 2 3 4 5]';
p=conv(d,w);

Now, let's pretend we don't know w. We want to estimate w, by deconvolution of p and d:

nw=length(p)-length(d)+1;   % length of w, based on lengths of d and p
D=convmtx(d,nw);            % matrix for deconvolution
wEst=inv(D'*D)*D'*p;        % estimate of w, by deconvolution
disp(wEst')

   -1.0000   -0.0000    1.0000    2.0000    1.0000    3.0000

Note that wEst equals w, which shows the deconvolution worked.