How to detect decrease of a digital signal?

Question

I have below given digital signal (sampled with period $T = 1\,\mathrm{ms}$) for which I need to evaluate its decrese. Namely I am interested in regions where the signal decreases faster than $-1000\,\mathrm{V}\cdot s^{-1}$.

My first idea how to do that was to evaluate the first derivative of the signal and compare it against zero. Due to the fact that the signal is noisy I have been filtering the signal via the first order low pass filter with time constant $\tau = 10\,\mathrm{ms}$ at first

$$ y_1(k) = \frac{T}{T + \tau}\cdot x(k) + \frac{\tau}{T + \tau}\cdot y(k-1) $$

The processing method can be basically described via following block diagram

The first derivative is evaluated via below given difference equation

$$ y_2(k) = \frac{1}{2\cdot T}\cdot\left[y_1(k) - y_1(k-2)\right] $$

corresponding to the average of the two consecutive first differences. The proposed processing method gives below given first derivative in following form.

I have doubts regarding the proposed method because the derivative is also little bit noisy despite the fact it is evaluated from the filtered signal.

Here I can see a possible solution consisting of further decrease the pass band of the low pass filter or apply some hysteresis band around zero in the comparator. Nevertheless both of these approaches seems to me to be fragile solution. My question is whether you know any better method suitable for signal decrease detection. Thanks in advance for any ideas.

Answer 1

I agree with Ben: use the CUSUM algorithm first up and see if that meets your needs.

If I do a rudimentary attempt at simulating your signal and implementing CUSUM, then I get the output shown below. The orange line is the threshold.

You'll need to play with $N$, the length over which the statistic is calculated and $\kappa$ the alarm / threshold parameter.

---

Code Below

import numpy as np
import matplotlib.pyplot as plt 
import scipy

# Let's try a CUSUM implementation for change in 
# mean following Basseville Example 2.1.1

mu0 = 800.0
mu1 = 0
Sigma = 10
NumSamples = 1024 # Number of samples
DropSample = int(NumSamples/3)
DropDuration = int(NumSamples/5)

x = np.concatenate((mu0*np.ones(DropSample),   
    DropLevel*np.ones(DropDuration), 
    StartLevel*np.ones(NumSamples-DropSample-DropDuration)))

xn = x + np.random.normal(0,NoiseStdDev,len(x))


def sufficientStatistic(yi, mu0, mu1, sigma):
    return (mu1-mu0)/sigma/sigma*(yi-(mu0+mu1)/2)

def S1N(yN, mu0, mu1, sigma):
    N = len(yN)
    return np.sum(sufficientStatistic(yN, mu0, mu1, sigma))

def decisionFunction(y, N, mu0, mu1, sigma):
    K = int(len(y)/N)
    values = np.zeros(K)
    
    for k in range(K):
        values[k] = S1N(y[range(k*N, (k+1)*N)], mu0, mu1, sigma)
        
    return values

N = 8
decisionValues = decisionFunction(xn, N, mu0, mu1, Sigma)    
kappa = 100
threshold = mu0 - kappa*Sigma/np.sqrt(N)

plt.figure(3)
plt.plot(decisionValues)
plt.plot([0,len(decisionValues)], [threshold, threshold])