detect to rising, stable and falling point in non-smooth rectangular wave

Question

I am working on basic signal processing problems in MATLAB. I have found a signal from the internet (i don't remember the site exactly). The data is organized in column wise. 1st and 2nd column is data and 3rd column is a rectangular wave. The rectangular wave is not very smooth. So I decided to detect four points as shown in the figure. I spend 2 days to find the exact solution but i couldn't find any solution to the detect those four points. I tried few things:
1) finding the difference between the points to get 1,2,3 and 4 points. but it doesn't work.
2) Median Filter to smooth the wave.
Anyone, please give some suggestion to detect these points in MATLAB

Answer 1

The usual approach to change detection is the CUSUM algorithm.

I've done an implementation that just addresses the level (mean) change issue. It's included (in R) below.

The black line is the noise-free data, the red line is the noisy data and the blue bars are the detected breaks (for this realization).

This just addresses the level change; to address the three zones: low level, changing level, and high level, you'll need to figure out the way to estimate the mean in the changing level (perhaps assume a rise time and fix the two levels?).

The code below is based on this document.

A slight change might get you the first and second changes easily: instead of working with the data, work with the difference of the data. If I do that, then I get:

but I've had to add:

sigma <- 0.001
data_before_diff <- noiseless_data + rnorm(length(noiseless_data),0,sigma)
data <- diff(data_before_diff)
mu2diff <- 0.04
thresh <- 1

and then later:

 s[k] = (mu2diff - mu1)/sigma*(data[k] - (mu2diff + mu1)/2)

and

    tmp <- mu1
    mu1 <- mu2diff
    mu2diff <- tmp

To find the changes at the end of the signal, you'll probably need to change mu2diff to -mu2diff and check for that versus 0.

Note I've had to greatly reduce the noise variance and thus the threshold to get this to work.

---

R Code Below

# 30039
# CUSUM 

N <- 30

mu1 <- 0
mu2 <- 1
noiseless_data <- c(rep(mu1,N), seq(mu1,mu2,1/N), rep(mu2,N*4), seq(mu2,mu1,-1/N), rep(mu1,N))


sigma <- 0.1
data <- noiseless_data + rnorm(length(noiseless_data),0,sigma)

thresh <- 100

breaks <- rep(0,N)
num_breaks <- 0

s <- rep(0,length(data))
capS <- rep(0,length(data))
G <- rep(0,length(data))
for (k in 1:length(data))
{
  s[k] = (mu2 - mu1)/sigma*(data[k] - (mu2+mu1)/2)
  if (k==1)
  {
    capS[k] = s[k]
  }
  else
  {
    capS[k] = capS[k-1] + s[k]
  }

  G[k] <- max(0,capS[k] - min(capS[1:k]))

  if (abs(G[k]) > thresh)
  {
    #capS[k] <- 0
    num_breaks <- num_breaks + 1
    breaks[num_breaks] <- which.min(capS[1:k])
    tmp <- mu1
    mu1 <- mu2
    mu2 <- tmp
  }
}

first_break <- min(breaks[1:num_breaks])
last_break <- max(breaks[1:num_breaks])

plot(data,col="red", type="l")
lines(noiseless_data)
lines(c(first_break,first_break),c(0,1),col="blue", lwd=10)
lines(c(last_break,last_break),c(0,1),col="blue",lwd=10)

Answer 2

@Peter I have trief to implement your code in MATLAB, but the output is not same as yours. Can you please suggest why is it so? Second thing I want to know is why you used difference data?

data_before_diff <- noiseless_data + rnorm(length(noiseless_data),0,sigma)
data <- diff(data_before_diff)

MATLAB CODE

N=30;

mu1=0;
mu2=1;
noiseless_data=[mu1*ones(mu1,N), mu1:1/N:mu2, mu2*ones(1,4*N), mu2:-1/N:mu1, mu1*ones(1,N)];

sigma=0.01;
data=noiseless_data+ rand(1,length(noiseless_data));

thresh=100;

breaks=zeros(1,N);

num_breaks=0;

s= zeros(1,length(data));

capS=zeros(1,length(data));

G=zeros(1,length(data));
for k =1:length(data)
    s(1,k) = (mu2 - mu1)/sigma*(data(1,k) - (mu2+mu1)/2);
    if (k==1)
      capS(1,k) = s(1,k);

    else

        capS(k) = capS(k-1) + s(k);

    end
    G(k) = max(0,capS(k) - min( capS(1:k)));

    if (abs(G(k)) > thresh)

%         capS(k)= 0
        num_breaks=num_breaks + 1;
        breaks(num_breaks) = find(min(capS(1:k)));
        tmp = mu1;
        mu1 = mu2;
        mu2 = tmp;
    end
end


first_break = min(breaks(1:num_breaks));
last_break = max(breaks(1:num_breaks));

plot(data,'r')
hold on
plot(noiseless_data)
line([first_break first_break],[0,1],'LineWidth',4,'Color','b')
line([last_break last_break],[0,1],'LineWidth',4,'Color','b')
hold off