filters - Duration of unknown rectangular pulse with additive white Gaussian noise

There is a discrete signal $f[i]$ (example below).

It is known, that $f[i]$ have a form of rectangular pulse with additive white Gaussian noise.

$f[i] = s[i] + n[i]$,
$s[i] = \alpha(\theta[i - i_{1}] - \theta[i - i_{2}]) + c$,
$i_{2} > i_{1}$,
$i_{1} > N$

Where
$\theta[i]$ is a Heaviside step function,
$n[i]$ is an additive white Gaussian noise,
$\alpha$ is a height of rectangular pulse,
$i_{1}$ is an index of the first sample of rectangular pulse,
$i_{2}$ is an index of the last sample of rectangular pulse,
$c$ is constant level of signal,
$N$ is adjustable parameter.

All parameters may have large range of values.
It is required to find value of $(i_{2} - i_{1})$ (duration of rectangular pulse in samples).

At the moment, I have tried two ways to solve this problem.

Low pass filter with threshold.

As first attempt, I have used simple scheme with low pass filter and threshold.
1. Apply FIR low pass filter with cutoff frequency equal to $0.05f_{sampling}$.
2. Estimate mean $m$ and dispersion $\sigma^{2}$ of filtered noise from first $N$ samples of the signal.
3. Set threshold $t = m + 3\sigma$.
4. Estimate $i_{1} = \min_{i}(f[i] > t)$.
5. Estimate $i_{2} = \max_{i}(f[i] > t)$.

Pros:
1. This algorithm is simple.
2. It is easy to write fast implementation.

Cons:
1. It is hard to estimate efficient value of filter's cutoff frequency. On the one hand low value can corrupt the form of short pulses. On the other hand large value decreases effect from filtration.
2. Algorithm isn't using all information, we have about the signal.

Regression analysis

As the second attempt, I have tried to approximate input sequence of samples with the function $s'[i]$.

$s'[i] = \alpha(\theta'[i - i_{1}] - \theta'[i - i_{2}]) + c$,
$\theta'[i] = \frac{1}{1+e^{-\frac{i}{\beta}}}$, where $\beta$ is a small parameter.

For approximation I have used the method of least squares with gradient descent for minimizing cost function.
1. Set initial values for $\alpha$, $c$, $i_{1}$, $i_{2}$.
2. Perform gradient descent.
3. Set threshold $t = c + 0.5\alpha$.
4. Estimate $i_{1} = \min_{i}(f[i] > t)$.
5. Estimate $i_{2} = \max_{i}(f[i] > t)$.

Pros:
1. This algorithm gives results with good precision.
2. It works for wide range of durations.

Cons:
1. It is very slow.

After all, I am not satisfied with the precision of the first algorithm and with the speed of the second one. How would you solve this problem?
Is there any classical solution, that I failed to find?
Ideas, links, any feedback will be much appreciated.
Thank you.

Answer

You want a method that removes noise while preserving edges. This cannot be achieved well by linear filtering, as you noticed yourself. I know of two approaches that might work well for your problem. The first is median filtering, where samples inside a window are replaced by their median. The following plot shows the result of median filtering with a window length of 25 samples (in red):

The other, more complex, approach is Total variation denoising, which works very well for piecewise constant signals. There is a very good description of total variation denoising including Matlab code available: link.