Can anyone help me in understanding of deconvolution based on toeplitz matrix?
85 views (last 30 days)
Show older comments
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!
0 Comments
Answers (1)
William Rose
on 25 Sep 2024 at 16:01
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)
Use D to compute p=convolution of d with w:
p=D*w;
disp(p')
Compute estimate of w, from p, using the Toeplitz matrix:
wEst=inv(D'*D)*D'*p;
disp(wEst')
The result shows that the estimate of w equals the original w.
2 Comments
Paul
on 25 Sep 2024 at 18:17
Hi Wiliam,
w=[-1 0 1 2 1 3]';
d=[1 2 3 4 5]';
D = convmtx(d,numel(w))
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!