CODEWITHSUNDEEP
matlabconvolution
jfeng92
jfeng92

Reputation: 341

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

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!

Upvotes: 1

Views: 14496

Answers (3)

Mehrwolf
Mehrwolf

Reputation: 8537

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.

Upvotes: 2

Dan Becker
Dan Becker

Reputation: 1214

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.

Upvotes: 2

Dan
Dan

Reputation: 45752

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?

Upvotes: 1

Related Questions