Faster Algorithm to convolve/correlate two sparse 1-D signals in python (or any language)

Question

I have two signals which I need to correlate or convolve. Each signal is sampled non-uniformly and the values of the signal I have with me are the timestamp and the magnitude of the signal at that timestamp. The values of the signal at all the other times can be assumed to be zero. The timestamps of the signal have a resolution running in microseconds.

An example of how the signal looks is shown below:

As can be seen, the resolution of the signal is in microseconds and signal is mostly sparse.

If I were to the convolve two signals of this type, I would first have to pad the signals with zeros (since I would have to discretise the signal). While the padding can be done with resolution of microseconds, The number of values to be multiplied becomes too big and the operation becomes increasingly slow. Most of the multiplications in this convolution would be multiplication of zeros(which are pretty much useless). I have therefore chosen a round off value of 2 places (0.xxxxxx becomes 0.xx),since I have to perform 40,000 similar convolutions. I have written my resampling function as shown below.

import numpy as np
import math

def resampled_signal_lists_with_zeros(signal_dict, xlimits):
  '''
      resamples the given signal with precision determined by the round function. 
      signal_dict is a dictionary with timestamp as key and signal magnitude is the value of the key. 
      xlimits is an array containing the start and stop time of the signal.
  '''
  t_stamps_list = list(signal_dict.keys())
  t_list = list(np.arange(int(math.floor(xlimits[0])), int(math.ceil(xlimits[1])), 0.005))
  t_list = [round(t, 2) for t in t_list]
  s_list = list()
  time_keys = [round(t, 2) for t in t_stamps_list]

  i = 0
  for t in t_list:
    if i < len(t_stamps_list):
      if t==time_keys[i]:
        s_list.append(signal_dict[t_stamps_list[i]])
        i+=1

      else:
        s_list.append(0)
    else:
      s_list.append(0)
  return t_list, s_list

The correlation of two signals padded in the above manner is done using scipy as follows:

from scipy.signal import correlate
output = correlate(s_1, s_2, mode='same')

The output calculated in the above manner is pretty slow .Since the signal is pretty sparse and most of the multiplications in the signal are multiplications of zero, I think there should be a better way to do the same operations. Is there a way to get the result of the convolutions of the two sparse signals faster?

Answer 1

Interpolating discontinuous waveforms is usually not a good idea.

The way I would approach your problem would be to recognize your signals as pulse trains. So I would assume that the convolution a single pulse from one waveform with a single pulse from the other waveform was also a pulse.

so signals A and B could be represented as a train of continuous Dirac functions. $$ s_A(t)=\sum_{i=1}^{N-1} a_i \delta(t-\tau_A(i)) $$ $$ s_B(t)=\sum_{i=1}^{M-1} a_i \delta(t-\tau_B(i)) $$ where $\tau_A(i)$ and $\tau_B(i)$ are your non uniform arrival times.

There is an identity for delta functions $$ f(t) \star \delta(t-a)= f(t-a) $$ where $\star$ denotes convolution, so $$ s_A(t)\star c_1\delta(t-a)= c_1\sum_{i=1}^{N-1} a_i \delta(t-(\tau_A(i)+a)) $$ and $$ s_A(t)\star \left[c_1\delta(t-a) + c_2\delta(t-b)\right]= c_1\sum_{i=1}^{N-1} a_i \delta(t-(\tau_A(i)+a))+c_2\sum_{i=1}^{N-1} a_i \delta(t-(\tau_A(i)+b)) $$ which can be generalized to $s_A(t)\star s_B(t)$

you really just need calculate the time shifts and amplitude products.

I will leave it up to you how to handle the amplitudes but the code below assumes all the impulses are magnitude one.

clear all
close all
clc
%%
M=2 ; % number of spikes in waveform A
N=5; % number of spikes in waveform B
T=10; %  arbitrary duration
%%
waveform_a_times=rand(1,M)*T;
waveform_b_times=rand(1,N)*T;

%%
for i=1:N
    time_diff(i,:)=waveform_a_times+waveform_b_times(i);
end
conv_times=sort(reshape(time_diff,1,M*N));
figure(1)
subplot(3,1,1),stem(waveform_a_times,ones(1,M));
title('waveform A')
subplot(3,1,2),stem(waveform_b_times,ones(1,N));
title('waveform B')
subplot(3,1,3),stem(conv_times,ones(1,N*M));
title('convolution of waveform A with waveform B')

for correlation

$$ f(t)\star \delta(t+a)=f(t+a)$$

you’re also more likely to have coincident times for correlation which you should test for.