How to do an unusual convolution...

So, I have three things which I need to convolve together.

1. - A list of times of occurrences and amplitude of a signal. Each item in the list can vary in time from 0 to 10^3. These values are produced by an external program.

2. - A function describing the precise form of the signal, which is best measured in units of 10^-6.

3. - The analytical result of passing a unit impulse into filtering electronics.

I require the entire waveform being produced by the electronics, 10,000 times, quickly, accurately and without using a great deal of memory to store the output.

Currently, I perform a time simulation of the electronics as it is the only way I have of maintaining the requisite time-accuracy, and handling the irregular distribution of the list of occurrences without taking excessive time and overloading the memory burden of the program. Any ideas?

[José-Luis]

How to ask a question

[Matt J]

What does it mean to convolve with a "list", your data in item 1? And why isn't the obvious solution of CONV function or fft-based convolution using FFT applicable to you?

[Image Analyst]

If the list items (time points) are not uniformly spaced, you will have to first call interp1() to get uniformly spaced times, then call conv().

[Iain]

The list, is a list of "unit impulses". That is infinitely thin, infinitely high pulses that can be integrated to 1, but are multiplied by their amplitudes. You can think of the list as a badly formatted "sparse" matrix.

The obvious solution of conv is too memory intensive. - In the worst case, I could have a vector extending from 0 to 1000, in steps of a ten-millionth. Thats 40GB, if I use "singles". Clearly this is beyond my system's capabilities.

[Iain]

The time points are neither uniformly spaced nor guaranteed to be unique.

[Image Analyst]

Some data or a screenshot would help.

[Iain]

The data can look like:

In a benign case:

time = 500 + 3*randn(1000,1); amplitude = 10e-13 + randn(1000,1)*4e-14;

In a typical case

time = [500+3*randn(1000,1); 550]; amplitude = [10e-13+randn(1000,1)*4e-14; 1e-9];

In the worst case

time = [0 500+3*randn(1000,1); 1000]; amplitude = [1e-8; 10e-13+randn(1000,1)*4e-14; 1e-9];

The values are the result of an extensive and highly computationally expensive process being run in a different (non-matlab) application. The benign case will usually happily work with conv, but the benign case is the least valuable.

[Matt J]

The obvious solution of conv is too memory intensive. - In the worst case, I could have a vector extending from 0 to 1000, in steps of a ten-millionth. Thats 40GB, if I use "singles". Clearly this is beyond my system's capabilities.

That problem isn't related to CONV. You're going to need 40GB just to hold the input data and the result, no matter what you do. What happened to the 10^-6 sampling interval that you mentioned initially? How did it drop down to 10^-7?

[Iain]

I do have the option of storing the result as a "sparse" vector. As I attempted to illustrate in the "worst" case, the data does contain protracted periods of no activity.

By "2. - A function describing the precise form of the signal, which is best measured in units of 10^-6." I meant that the entire DURATION is best measured in units of 10^-6. The sampling frequency, therefore, would need to be more than 10 times that. (10^-7.) - Obviously though, even 10^-6 is problematic (4GB)

Answer 1

[Matt J]

I do have the option of storing the result as a "sparse" vector. As I attempted to illustrate in the "worst" case, the data does contain protracted periods of no activity.

It's not entirely clear from the discussion which data is sparse and which data is not. I understand that your "list-mode" data set is sparse, i.e., the data set (1) that is held as a time/amplitude list. However, you have not indiicated that the data sets in (2) or (3) are sparse.

If you're saying that all signals being convolved with each other are going to be sparse, you could try using the CONVN method of my ndSparse class located here.

If (2) and (3) are not sparse, then their convolutions with the list mode data are not going to be sparse either. They will inevitably consume 4-40GB.

[Iain]

Datasets 2 and 3 are not sparse, as they are relatively short period (compared to the first one), and are continuious. However, there's no real reason to not cast them as sparse to try that method.

Answer 2

[Matt J]

Well, first of all I would bin your time/amplitude "list" data to the same sampling intervals as your other signal. If the signal is "continuous" when sampled at 1e-6, it's probably not going to be necessary to have time resolution of the "list" data any better than that.

 L=length(signal);
 timeAxis=(0:L-1)*1e-6;
 [~,sub]=histc(time, timeAxis);

 impulses=accumarray(sub,amplitude(:),[L,1]);

Now you would convolve "signal" with "impulses". One way is with FFTs

 out= ifft(fft(impulses,2*L).*fft(signal,2*L),L,'symmetric');