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

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

• @Peter i have deleted the comment which you though as off-topic – Aadnan Farooq A Apr 11 '16 at 17:20

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:

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)

• Thanks for the suggestion. Actually, I want to detect the 4 points 1) point where it start rising 2) point where it start stable 3) point where falling start 4) point where faling ends – Aadnan Farooq A Apr 11 '16 at 16:29
• @AadnanFarooqA Hence 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 in my answer. – Peter K. Apr 11 '16 at 16:56
• I am sorry i have deleted my comment – Aadnan Farooq A Apr 11 '16 at 17:02
• $c(rep(mu1,N), seq(mu1,mu2,1/N), rep(mu2,N*4), seq(mu2,mu1,-1/N), rep(mu1,N))$ 'c' was not defined before so what is that? – Aadnan Farooq A Apr 11 '16 at 17:07
• @AadnanFarooqA In matlab that'd just be creation of a vector: [mu1*ones(mu1,N), mu1:1/N:mu2, mu2*ones(1,4*N), mu2:-1/N:mu1, mu1*ones(1,N)]; (i.e. c is like [] in matlab). – Peter K. Apr 11 '16 at 17:47

@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

• I use "diff" data because that means what we are looking for are the edges in a signal that looks like : ________| |________| |__ rather than a linearly increasing change. The constant pieces become zero and the linearly changing pieces become non-zero sections. – Peter K. Apr 13 '16 at 21:25
• @PeterK. How can i generate the noiseless data from the noise step data. As i have just noisy step data.. – Aadnan Farooq A Apr 14 '16 at 5:48
• 2) When I random data then in some case i get the output as you showed. But when I fix the data than I cannot get the output. I have updated the code with your second suggestion – Aadnan Farooq A Apr 14 '16 at 6:10