Moving mean square error between 2 arrays, 'valid', where they fully overlap

user_8324480 0
Viewed 161

I have a noisy square signal which looks like this:

square signal

The amplitude is known. To match the complete square, I can create a pattern of the square and apply np.correlate to find where the signal and the pattern maximally overlap. I wanted to apply a similar approach to find the edge, try to correlate with the 2 patterns below:

pattern

As the correlation is nothing more than a convolution, this doesn't work. Half the pattern is equal to 0, and the convolution of this half will return 0 no matter the position on the signal; while the other half is equal to -X with X the amplitude. This second half convoluted with the signal will be maximal when the signal amplitude is maximal. On the signal plot, you can observe that the square is not perfect and that the beginning has a slightly larger amplitude. Basically, both correlation leads to a match on the beginning of the square, where the convolution is maximal. The ramp up (end of the square) is not detected.

To avoid this problem, I would like to use a different operation. As I do know the amplitude of the square signal, I can generate a pattern with the correct amplitude, in this case about -0.3. Thus, I would like to take the pattern and slide it across the signal. At each step, I would compute the mean square error and my pattern would match with the signal at the position where the mean square error is minimized. Moreover, I would like to use the same type of setting as for a convolution, 'valid', where the operation is performed only when the 2 arrays fully overlap.

Do you know of an other method; and/or which function, methods I should use? I couldn't find a all-in-one function line np.convolve or np.correlate.

EDIT: Since I couldn't find a pre-coded function in a library, I've coded mine with a while loop. It's pretty inefficient... It's up here on codereview for upgrades.

convolution correlation numpy python
2 Answers

I think that convolving/correlating your signal with a step function is still pretty close to the optimal solution, since this is similar to matched filtering, which can be proven to be optimal (under certain conditions, noise likely needs to be Gaussian).

The only issue you have is that your template (the step function) contains a DC part. Removing this, will give you the result you want:

import numpy as np
import matplotlib.pyplot as plt


# simulate the signal
signal = np.zeros(4000)
signal[200:-400] = -0.3
signal += 0.005 * np.random.randn(*signal.shape)

plt.plot(signal)
plt.title('Simulated signal')
plt.show()


# convolve with template with non-zero DC
templ = np.zeros(200)
templ[100:] = 1  # step from 0 to 1
plt.plot(np.convolve(signal, templ))
plt.title('Convolution with template with DC component')
plt.show()

# convolve with template without DC
templ_ac = templ - templ.mean()  # step from -0.5 to +0.5
plt.plot(np.convolve(signal, templ_ac))
plt.title('Convolution with template without DC component')
plt.show()

Results:

enter image description here enter image description here enter image description here

The way to understand this is that convolve(signal, template) = convolve(signal, template_DC) + convolve(signal, template_AC), where template_DC = mean(template) and template_AC = template - template_DC. The first part is the convoltion of the signal with a flat template, which is just a smoothed version of your signal. The second part is the 'edge detection' signal you want. If you do not subtract the AC part of the template, the uninteresting first part dominates the interesting part.

Note that the scaling of the template is not important, the step in the template doesn't have to be 0.3. This will just cause a scale factor in the end result. Also note that this method does not depend on the exact value of the step, so a larger step in your signal will cause a large effect in the edge detection.

If you know that the step is always exactly 0.3, and you want to be insensitive to steps of different amplitude, you could do some sort of least square fitting of the signal with the template, for every possible shift of the templatate, and only trigger a detection of an edge if the residual is small enough. This will be slow, but might give better rejection of steps with the wrong amplitude.

user_26472790

Since you have very little noise you can calculte where the signal drastically changes with a loop, for example:

for i in range(begin+10,end):
    if(abs(data[i-10]-data[i])>0.1):
        foundChange()
user_134271550
Related