We had certain discussions about the issue that you observe.
The behaviour that you observe is not a bug, but it is caused by limits of numerical precision. The deconv function performs a polynomial division and for long inputs, this can run into problems because the generated polynomial coefficients that are necessary to perform the division can get very large. Namely, since the denominator coefficient vector in the call to “filter” within deconv is very long, things can “blow up” and the returned coefficients can become Inf and NaN when the polynomial division becomes numerically unstable.
I assume that by “overlapping convolutions”, you mean that the non-zero entries of the Events array are spaced far enough apart that an entire copy of the Template signal can exist between events. In this case, the effect of the denominator coefficients has time to “wind back down” in the polynomial division. Put another way, events spaced closely together make the polynomial division very numerical unstable. This is an unfortunate weakness of computing a deconvolution via pure polynomial division. The same is true for a repeated sequence of events, though it depends significantly on how the frequency at which the sequence is repeated.
There are a few options for alternative approaches that you can use. Since the signal is calculated using the full linear convolution of Events with Template, the deconvolution can be achieved with the FFT:
RestoredEvents = ifft(fft(Signal) ./ fft(Template, length(Signal)), 'symmetric');
RestoredEvents = RestoredEvents(1:length(Events));
This is, in general, a numerically robust approach unless the template signal is very spectrally concentrated (i.e., if it has energy in only a small band of frequencies). If anything ever “blows up”, it is because the transform of the Template signal contains coefficients that are very small in value. This can be countered by a regularization of the denominator coefficients.
It’s important to note, however, that this only gives a good answer because for this case we are certain that the Signal that we have is exactly the linear convolution of the template signal with the series of events (and the template signal has good spectral properties). However, this may not always be the case. The suspicion is that the larger goal here is to take an observed signal and attempt to determine the locations of sparse events that trigger the template response, possibly in the presence of noise. More importantly, I suspect that using circulant convolution may not be an appropriate way to model this problem if it is more general than the example that was given.
Another alternative approach to attempt to discover the events is to set this problem up as a linear convolution system:
T = toeplitz([Template(:); zeros(numel(Signal)-numel(Template), 1)], [Template(1), zeros(1, length(Events)-1)]);
RestoredEvents = T \ Signal(:);
This is a reasonable approach, as it obtains the deconvolution by solving a linear system to return the least-squares deconvolution of the two inputs. However, it’s an inefficient tactic because it requires the explicit construction of a (large) Toeplitz matrix. If the Template is known to be short relative to the length of Signal, then there is possibly some advantage to constructing this as a sparse matrix.
Hope this helps!
