2
$\begingroup$

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:

result of a convolution of a signal with the gaussian (1,1)

result of a convolution of a signal with the gaussian (2,1)

result of a convolution of a signal with the gaussian (3,2)

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.

after demean and rescale after the convolution

$\endgroup$

1 Answer 1

1
$\begingroup$

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.

normalized correlation result

$\endgroup$
6
  • $\begingroup$ Dear Dan, what I did - I have changed my procedure of matched filtering using your proposed approach - I have subtracted DC offsets from template and signals, convolved them and then divided the result on the product of SDs. But in this case I am getting the same result rescaled - so the amplitude of the result is greater for the greater amplitude (or power) of the input signal, but not where we have more similarity. - see updated version $\endgroup$
    – twistfire
    Feb 17, 2023 at 13:11
  • $\begingroup$ @twistfire see my update. The issue was it is just the portion of "signal" that is within view of the matched filter that would need to be processed as I described, resulting in a true normalized cross-correlation coefficient for each output. Without doing that, the output of the matched filter will have the optimum SNR (assuming independent identically distributed noise on each sample) but will also be arbitrarily scaled--- so your other option is to compute and compare the SNR on the matched filter output to confirm the "matched" one is maximum for SNR. $\endgroup$ Feb 19, 2023 at 15:32
  • $\begingroup$ But the way I personally compute SNR for such a condition IS to use the normalized cross correlation coefficient to the "noise free' signal, perhaps with time interpolation first to remove the effects of a small time offset in the sampled signal. So that would basically be back to doing exactly what I did here. $\endgroup$ Feb 19, 2023 at 15:32
  • $\begingroup$ Dear Dan Boschen, thank you for your great explanations. Can you clarify 2 moments? Am I correct that to align the result of the proposed code-snippet and initial signal I can do the back-shift of the result by the time_shift = int(coeff.size/2-1) ? Another thing I wanted to ask - in case I use np.correlate(a=template, v=signal, mode='full')? Why am I not getting the same result as your code gives? I mean I got only positive numbers and the result is actually much broader than the pattern is - so it actually works as a smoothing filter for a signal. $\endgroup$
    – twistfire
    Feb 22, 2023 at 4:45
  • $\begingroup$ @twistfire If (and only if) the template is symmetric, then the delay is (coeff.size-1)/2. So for an even length size you can have a delay that has an offset of 1/2 a sample. np.correlate is not the normalized correlation, however if you subtract the mean and divide by the standard deviation (over just the length of the template), you will get the same answer. $\endgroup$ Feb 22, 2023 at 5:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.