The deconv() function in MATLAB does not invert the conv() function

Question

I would like to convolve a time-series containing two spikes (call it Spike) with an exponential kernel (k) in MATLAB. Call the convolved response "calcium1". I would like to recover the original spike ("reconSpike") data using deconvolution with the kernel. I am using the following code.

k1=zeros(1,5000);
k1(1:1000)=(1.1.^((1:1000)/100)-(1.1^0.01))/((1.1^10)-1.1^0.01);
k1(1001:5000)=exp(-((1001:5000)-1001)/1000);
k1(1)=k1(2);

spike = zeros(100000,1);
spike(1000)=1;
spike(1100)=1;

calcium1=conv(k1, spike);
reconSpike1=deconv(calcium1, k1);

The problem is that at the end of reconSpike, I get a chunk of very large, high amplitude waves that was not in the original data. Anyone know why and how to fix it?

Thanks!

score 1 · Answer 1 · answered Aug 22 '12 at 14:28

1

It works for me if you keep the spike vector the same length as the k1 vector. i.e.:

    k1=zeros(1,5000);
    k1(1:1000)=(1.1.^((1:1000)/100)-(1.1^0.01))/((1.1^10)-1.1^0.01);
    k1(1001:5000)=exp(-((1001:5000)-1001)/1000);
    k1(1)=k1(2);

    spike = zeros(5000, 1);
    spike(1000)=1;
    spike(1100)=1;

    calcium1=conv(k1, spike);
    reconSpike1=deconv(calcium1, k1);

Any reason you made them different?

answered Aug 22 '12 at 14:28

Dan

45,079
17
88
157

I made it different because I want the freedom to have spikes occur outside the time [0:5000]. For example, if I wanted at spike at t=6000. But since all values greater than 5000 are zero, I am not sure why I have this problem.Perhaps if I increase the length of k1 and just set it all to zero, I could solve this. – jfeng92 Aug 22 '12 at 14:49
Unfortunately, increase the length of k1 to 100000 does not work and I still get these waves. Perhaps this is a flaw inherent in the deconv function. Could you think of any workaround? – jfeng92 Aug 22 '12 at 14:56
A higher sampling frequency maybe? So instead of 1:1000 try 1:0.01:1000... otherwise I'm not really sure I'm afraid. – Dan Aug 22 '12 at 15:11

score 1 · Answer 2 · answered Aug 22 '12 at 21:42

You are running into either a problem with MATLAB's deconvolution algorithm, or floating point precision problems (or maybe both). I suspect it's floating point precision due to all the divisions and subtractions that take place during the deconvolution, but it might be worth contacting MathWorks directly to ask what they think.

Per MATLAB documentation, if [q,r] = deconv(v,u), then v = conv(u,q)+r must also hold (i.e., the output of deconv should always satisfy this). In your case this is violently violated. Put the following at the end of your script:

[reconSpike1 rem]=deconv(calcium1, k1);
max(conv(k1, reconSpike1) + rem - calcium1)

I get 6.75e227, which is not zero ;-) Next try changing the length of spike to 6000; you will get a small number (~1e-15). Gradually increase the length of spike; the error will get larger and larger. Note that if you put only one non-zero element into your spike, this behavior doesn't happen: the error is always zero. It makes sense; all MATLAB needs to do is divide everything by the same number.

Here's a simple demonstration using random vectors:

v = random('uniform', 1,2,100,1);
u = random('uniform', 1,2,100,1);
[q r] = deconv(v,u);
fprintf('maximum error for length(v) = 100 is %f\n', max(conv(u, q) + r - v))
v = random('uniform', 1,2,1000,1);
[q r] = deconv(v,u);
fprintf('maximum error for length(v) = 1000 is %f\n', max(conv(u, q) + r - v))

The output is:

maximum error for length(v) = 100 is 0.000000
maximum error for length(v) = 1000 is 14.910770

I don't know what you are really trying to accomplish, so it's hard to give further advice. But I'll just point out that if you have a problem where pulses are piling up and you want to extract information about each pulse, this can be a tricky problem. I know some people who work on things like this, so if you want some references let me know and I will ask them.

score 1 · Answer 3 · answered Aug 23 '12 at 16:58

You should never expect that a deconvolution can simply undo a convolution. This is because the deconvolution is an ill-posed problem.

The problem comes from the fact that the convolution is an integral operator (in the continuous case you write down an integral int f(x) g(x-t) dx or something similar). Now, the inverse of computing an integral (the de-convolution) is to apply a differentiation. Unfortunately, the differential amplifies noise in the input. Thus, if your integral only has slight errors on it (and floating-point inaccuarcies might already be enough), you end up with a total different outcome after differentiation.

There are some possibilities how this amplification can be mitigated but these have to be tried on a per-application basis.

The deconv() function in MATLAB does not invert the conv() function

3 Answers3