Struggling to implement this simple algorithm for ECG wave delineation I'm trying my best but my maths isn't good enough to implement the algorithm as outlined in this paper in python. It for detecting the onset and offset of a wave on an ECG, and it's using a well validated method.

My data is in a numpy array.

The steps it uses are:

A. Computation of the envelope of the ECG

B. Computation of the auxiliary signal

C. Windowing

I think I've got step A. done:

def calculate_qrs_envelope(self):
    self.hilbert = np.imag(hilbert(self.lead_data['y']))
    self.envelope = np.sqrt(np.add(self.lead_data['y'],self.hilbert))
    self.plot(x=self.lead_data['x'], y=self.envelope, pen='b')


Where my data is stored in lead_data['y'].

If I plot it, it looks correct (I can't show the picture here as I don't have enough karma).

And I'm confident that I won't struggle with C.

However, B is difficult. It says in the paper:


I can calculate AS using a simple derivative, as so:

def calculate_auxiliary_signal(self):
    self.aux_sig = np.append([0],np.multiply(2,pow(np.diff(self.envelope),2)))
    self.plot(x=self.lead_data['x'], y=self.aux_sig, pen='g')


But I don't understand how to do it with regards to what it says about the parabolic fit.

My sampling frequency is 100 Hz. 

Would anyone be able to give me a hand with this maths? It would be a great help.
I cannot find any sense on the formula exposed in the paper, but as is, the expression is just a derivative approximation, with $r=r_0$. 

$$
x'_k=\frac 1{10}(2(x_{k+2r}-x_{k-2r})+(x_{k+r}-x_{k-r}))\\
y_k=2(x'_k)^2
$$

That is what the paper say. 

That is, you have a signal $x$, and you calculate $x'$ as a moving average, taking samples above and below the current index.

For any given signal $x:x_1...x_n$, you will straightforwardly be able to calculate $x'$ and $y$ without trouble.