Matched filter, convolution with signals of various patterns. Explanation of results

Question

I am trying to use matched filters for signal detection and advanced processing - to get additional information from the signal after the convolution with various templates.

The task I have is pretty simple - to find patterns with known forms in the signal I have as an input - but I need to be able to distinguish the results - e.g. detect signals with various width so the signals of similar forms but with another width will not be detected or result of filtering will have much lower amplitude comparing it to the result where we have the exact match of signal & pattern we are trying to detect.

The problem is that after the convolution of the signal & template, I have similar amplitudes of the resulting convolution for different shapes of input signals - so it's pretty hard to interpret - which signal is more similar to the template using convolution results.

The code I use in Python to do the matched filtering:

def matched_filter(template, signal):

    # removing the DC from signal and template
    # template = template - template.mean()
    # signal = signal - signal.mean()

    # Compute the conjugate of the template
    conjugate_template = np.conjugate(template)
    
    # Flip the template for zero-phase delay
    flip_template = np.flip(conjugate_template)
    
    # Compute the convolution of the signal and flipped template
    convolution = np.convolve(signal, flip_template, mode='same')
    
    return convolution

The code for toy-demonstration of a problem.

# Define the clean Gaussian template

# template 1 - Gaussian width (1,1)
template_size=10

template_mean1 = 5
template_std1 = 1
template_A1 = 1

template = template_A1 * np.exp(-(np.arange(template_size) - template_mean1)**2 / (2 * template_std1**2))


# template 2 - Gaussian width (2, 1) - high amplitude the same width
template_mean2 = 5
template_std2 = 1
template_A2 = 2

template2 = template_A2 * np.exp(-(np.arange(template_size) - template_mean2)**2 / (2 * template_std2**2))

# template 3 - Gaussian width (1, 2) - low amplitude, wider
template_mean3 = 5
template_std3 = 2
template_A3 = 3

template3 = template_A3 * np.exp(-(np.arange(template_size) - template_mean3)**2 / (2 * template_std3**2))

# Generate a sum of Gaussians with various amplitudes as the input signal

# first
input_mean1 = 10
input_std1 = 1
A1 = 1

# second
input_mean2 = 20
input_std2 = 1
A2 = 2

# third
input_mean3 = 30
input_std3 = 2
A3 = 3

# Generate a not-noised sum of Gaussians signal
not_noised_signal = A1 * np.exp(-(np.arange(40) - input_mean1)**2 / (2 * input_std1**2)) + \
               A2 * np.exp(-(np.arange(40) - input_mean2)**2 / (2 * input_std2**2)) + \
               A3 * np.exp(-(np.arange(40) - input_mean3)**2 / (2 * input_std3**2))

# add some noise
#input_signal = not_noised_signal + np.random.normal(0, 0.5, 40)
input_signal = not_noised_signal 

# Perform matched filtering with zero-phase delay
filtered_signal = matched_filter(template, input_signal)
filtered_signal2 = matched_filter(template2, input_signal)
filtered_signal3 = matched_filter(template3, input_signal)

# Plot the results

#signal 1 & template 1
plt.cla()
plt.plot(input_signal, label='Input Signal')
plt.plot(template, label='template A='+str(template_A1)+', s='+str(template_std1), linestyle='--')

plt.plot(filtered_signal, label='Filter result')
#plt.plot(not_noised_signal, label='Not Noised Signal')
plt.legend()
plt.show()

Convolving the modelled gaussians with the gaussians of various templates (various width and amplitudes) - I have very similar results - the resulting amplitude after the convolution is larger for parts of a signal with the largest amplitude, but not for the parts of a signal where you have exact match of the template form with the signal form. See the figures for explanation:

Why is that?

I suppose I am wrong in understanding the results of convolution. Does the convolution result (matched filtering) show how similar is current part of the signal under test with the template used? Or it's not like that?

In my case - I need to find the best matches - to show the positions where the signal has the most similarity with the template used. So the amplitude and form of a Gaussian in this case must be taken into account - the max. of filtering must be on the positions where the pattern really is. I need to get the positions of possible peaks, along with the information about some properties of a signal - like the parameters of the detected signal using the convolution/filtering result.

Gaussian is taken only as an example. A real application is very noised and the normal distribution of a noise is assumed,

Signals, that can be detected have the form of Lorentzians with various parameters - so in real application various amplitudes/widths are present in the signal - that is why I am trying to build dome logic to distinguish one peak type from another.

P.S. Another thing I need to understand - is how to deal with the nonuniform grid in this case? I have non-equidistant samples, dx <> const

x(i+1) - x(i) <> x(i) - x(i-1)  

While for the template I have only a model - how can I apply matched filtering in this case?

UPDATE 1. Using proposed approach,

def matched_filter(template, signal):

    # removing the DC from signal and template
    template = template - template.mean()
    signal = signal - signal.mean()

    temp_std = np.std(template)
    signal_std = np.std(signal)

    # Compute the conjugate of the template
    conjugate_template = np.conjugate(template)
    
    # Flip the template for zero-phase delay
    flip_template = np.flip(conjugate_template)
    
    # Compute the convolution of the signal and flipped template
    convolution = np.convolve(signal, flip_template, mode='same')

    convolution = convolution / (temp_std * signal_std)
    
    return convolution

But in this case I get "rescaled" version of what I had in the input signal, not the magnitude where we have the most similarity.

Answer 1

The purpose of matched filtering is to optimize the signal to noise ratio of the result under the condition of independent identically distributed noise in each sample (such as AWGN). However if you want to use it to compute a comparative correlation coefficient, then you could also do the following processing to make it equivalent to a normalized correlation (within +/-1 where 1 is an exact match independent of amplitude scaling (and -1 is an exact match with a change in sign, and 0 would be orthogonal or uncorrelated):

Subtract any DC offset (mean) of the signal that is within the length of the template, as well as the template prior to processing.

Compute the standard deviation of that portion of the signal and the standard deviation of the template.

Divide the result by the product of the two standard deviations.

This is what occurs within the function np.corrcoef() which returns a 2x2 result as the autocorrelation of the first sequence, cross-correlation of the first sequence with the second sequence, the cross-correlation of the second sequence with the first, and the autocorrelation of the second sequence. The matched filter output when normalized as suggested above would be the cross-corelation for each index of the input. I have created the function that would do this below, using deque as a FIFO buffer to move the signal through and for each index compute the normalized cross correlation using corrcoef:

def matched_filter(template, signal):

    # flip template and conjugate as coefficients of matched filter
    # (does nothing if real symmetric, but for general case)
    coeff = np.conj(template[::-1])
    
    # set up filter as a deque
    values = deque(np.zeros_like(coeff), len(coeff))
    result = []   # store results
    for sample in signal:
        # load next sample into matched filter
        values.appendleft(sample)
        # compute normalized cross correlation of signal within view of filter
        # and conjugated reversed template
        result.append(np.corrcoef(coeff, values)[0,1])
        
    return np.array(result)    

Testing with the OP's waveform, I get the following result consistent with my expectation: The first and second pulse have a resulting normalized cross-correlation equal to 1 when aligned, which makes sense since the pulses are identical other than a gain scaling. However the third pulse which is slightly wider results in a cross-correlation close to 0.9.