Rectangular or trapezoidal signal filtering. Beginner

Question

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.

Answer 1

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.