# Rectangular or trapezoidal signal filtering. Beginner

I have the following signal:

I know, that there are 2 states ('High' and 'Low') in my real signal. I need to convert my signal to the sequence of these states (Low from 2 to 5 second, High from 5.5 to 7 and so on).

Could you give me a piece of advice about classic and well established method to extract such signal from mine?

Median filter was as an example. Is it the most classic way to do it? Or you may provide some information for a better approach?

I don't like the median filter approach, because I should rather voluntary define kernel and 2 thresholds. I am searching for a method, with more natural parameters.

# UPDATE:

Some knowledge about signal:

Amplitude and mean duration of the state may vary from approx. 1 to 3 fold. But I have some knowledge about CV and skewness/CV. CV~0.5-0.9, skewness/CV~2. Distribution of States duration is gamma-like.

Depending on your noise model, a median filter might be quite good or not. Have a look at this example:

import scipy.signal

# from http://stackoverflow.com/questions/23289976/how-to-find-zero-crossings-with-hysteresis
def hyst(x, th_lo, th_hi, initial = False):
hi = x >= th_hi
lo_or_hi = (x <= th_lo) | hi
ind = np.nonzero(lo_or_hi)[0]
if not ind.size: # prevent index error if ind is empty
return np.zeros_like(x, dtype=bool) | initial
cnt = np.cumsum(lo_or_hi) # from 0 to len(x)
return np.where(cnt, hi[ind[cnt-1]], initial)

Fs = 300
T = 9
a = np.zeros(T*Fs)
# setup the state changes
a[0*Fs] = 1
a[0.5*Fs] = -2
a[2*Fs] = 2
a[3*Fs] = -2
a[4.5*Fs] = 2
a[7*Fs] = -2
a[8*Fs] = 2
a[8.5*Fs] = -2
# create the rect-waveform
data = np.cumsum(a)

# add some noise to the signal
rx = data + 2*np.random.randn(len(data))

plt.subplot(121)
plt.plot(rx)
plt.plot(data, 'r', lw=2)

plt.subplot(122)
# plot the filtered signal (if you have access to the purple signal in your plot, just use this signal)
filtered = scipy.signal.medfilt(rx, Fs//2+1)
plt.plot(filtered, lw=3)

# Perform thresholding with Hysteresis
detected = 2*hyst(filtered, -0.3, 0.3)-1
plt.plot(detected)

# From the thresholds, get the state changes via the derivative
diffed = np.diff(detected)
print ("estimated times of state switch:", ["%.2f" % (x/Fs) for x in np.nonzero(diffed)[0]])
print ("actual times of state switch:", ["%.2f" % (x/Fs) for x in np.nonzero(a)[0]])


estimated times of state switch: ['0.18', '0.60', '2.05', '3.10', '4.56', '7.07', '8.07', '8.59']
actual times of state switch: ['0.00', '0.50', '2.00', '3.00', '4.50', '7.00', '8.00', '8.50']


It creates an artificial signal, adds some noise, does a recovery of the original signal via median filtering. Eventually, it detects the state changes in the signal.

• Thanks, but I was mainly interested is there is something smarter, than median filter? – zlon Mar 29 '17 at 6:52
• What's the problem with median filter, if it works? – Maximilian Matthé Mar 29 '17 at 6:53
• I should define the size of the kernel and thresholds for each signal. I want some algorithm with a more natural parameters definition. – zlon Mar 29 '17 at 6:55
• You can fit a statistical model. Do you have access to the minimum/maximum time of a state, including the probability of each length (e.g. is it equally likely that a state remains for 5 seconds or 1 second)? Dou you have knowledge of the amplitude of the state? – Maximilian Matthé Mar 29 '17 at 7:08
• Amplitude and mean duration of the state may vary from approx. 1 to 3 fold. But I have some knowlege about CV and skewness/CV. CV~0.5-0.9, skewnes/CV~2. Distribution of States duration is gamma-like. – zlon Mar 29 '17 at 7:13