变化检测的常用方法是CUSUM 算法。

我已经完成了一个只解决水平（平均）变化问题的实现。它包含在下面（在 R 中）。

黑线是无噪声数据，红线是噪声数据，蓝条是检测到的中断（用于此实现）。

这只是解决了级别的变化；要解决三个区域：低电平、变化电平和高电平，您需要找出估计变化电平平均值的方法（也许假设一个上升时间并固定这两个电平？）。

以下代码基于此文档。

轻微的变化可能会让您轻松完成第一次和第二次更改：而不是处理数据，而是处理数据的差异。如果我这样做，那么我得到：

但我不得不补充：

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

然后：

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

和

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

要找到信号末尾的变化，您可能需要更改mu2diff为-mu2diff并检查它与 0 的对比。

请注意，我必须大大降低噪声方差，从而降低阈值才能使其正常工作。

下面的R代码

# 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)

@Peter 我想在 MATLAB 中实现您的代码，但输出与您的不同。你能建议为什么会这样吗？我想知道的第二件事是您为什么使用差异数据？

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

MATLAB 代码

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