Skip to content

Instantly share code, notes, and snippets.

Show Gist options
  • Save sdwfrost/e40be35990bfde49c665fe5d1455ed94 to your computer and use it in GitHub Desktop.
Save sdwfrost/e40be35990bfde49c665fe5d1455ed94 to your computer and use it in GitHub Desktop.
Change Detection Using Volatility Filters
#LMS change Detector
function [state,h]=ChangeDetectionLMSShiftMultipleSensor1(x,WinSlow,WinFast,WinDesired,muBase)
%Slow Response to change, accurate in steady state
[m,n]=size(x);
if m>n
x=x';
end
[NumSensors,n]=size(x);
shift=300;
xSlow=[zeros(NumSensors,WinSlow),x];
for i=1:length(xSlow)-WinSlow
OutSlow(:,i)=std(xSlow(:,i:i+(WinSlow-1))');
end
OutSlow=[zeros(NumSensors,shift),OutSlow(:,1:end-shift)];
%Fast Response to change, poor accuracy in steady state
xFast=[zeros(NumSensors,WinFast),x];
for i=1:length(xFast)-WinFast
OutFast(:,i)=std(xFast(:,i:i+(WinFast-1))');
end
%Fast Response to change, poor accuracy in steady state
xDesired=[zeros(NumSensors,WinDesired),x];
for i=1:length(xDesired)-WinDesired
OutDesired(:,i)=std(xDesired(:,i:i+(WinDesired-1))');
end
%Adaptive Filter sparse LMS convex combination of two variance estimation
%algorithms
epsilon=0.001;
omega=0;
state=zeros(1,length(x));
i=1;
future=WinSlow+shift+50;
for i=1:NumSensors
mu(i)=muBase/var(x(1:future));
end
count=1;
flag=0;
i=1;
W(1,1:NumSensors)=0;
while i<length(OutSlow)-WinDesired
WTemp=mean(W(i,:));
h(i)=WTemp;
for ns=1:NumSensors
yEstimate=WTemp*OutFast(ns,i) + (1-WTemp)*OutSlow(ns,i);
%yOut(i)=yEstimate;
variation=(randn(1)*epsilon);
e(i)=OutDesired(ns,i+WinDesired)-yEstimate;
%variation=epsilon;
W(i+1,ns)=WTemp + mu(ns)*e(i)*(abs(WTemp)-variation)*(OutFast(ns,i)-OutSlow(ns,i));
if W(i+1,ns)>1
W(i+1,ns)=1;
end
if W(i+1,ns)<0
W(i+1,ns)=0;
end
end
if h(i)>0.8 && flag==0 && i<(length(x)-future)
if state(i)==1
state(i+2:i+future)=1;
else
state(i:i+future)=1;
end
for f=1:NumSensors
mu(f)=muBase/var(x(f,i:i+future));
end
%flag=1;
flag=1;
count=1;
end
if flag==1
count=count+1;
% flag=0;
if count>future
flag=0;
end
end
i=i+1;
end
##Statistical Change Detector
function [ChangePoint,ActChange]=PeakChangeDetector(x,WinSlow,sigLevel)
xSlow=[zeros(1,WinSlow),x];
for i=1:length(xSlow)-WinSlow
OutSlow(i)=std(xSlow(i:i+(WinSlow-1)));
end
%Differencing of the volatility filter output
DiffVolTimeSeries=tsmom(OutSlow',WinSlow);
DiffVolTimeSeries(1:WinSlow)=0;
N=length(DiffVolTimeSeries);
count=1;
mu=sqrt(2)*gamma((WinSlow+1)/2)/(gamma((WinSlow)/2)); %Estimate of the first moment of the chi distribution
ChangePoint=zeros(1,N);
h=1;
ActChange=zeros(1,N);
%Loop to calculate the changePoints
while count<N+1
if count == 1
count=WinSlow*2;
VarEstimate=var(x(1:count));
StdEstFilter=2*sqrt(VarEstimate*(WinSlow-(mu^2))/WinSlow);
count=count+1;
end
if abs(DiffVolTimeSeries(count))>StdEstFilter*sigLevel && count+WinSlow<N
VarEstimate=var(x(count:count+WinSlow));
StdEstFilter=2*sqrt(VarEstimate*(WinSlow-(mu^2))/WinSlow);
%store(h)=count;
ActChange(count)=1;
h=h+1;
[~,I] = max(abs(DiffVolTimeSeries(count:count+WinSlow)));
ChangePoint(count+I-1-WinSlow)=1;
count=count+WinSlow*2;
else
count=count+1;
end
end
##Generating Random Change Point locations
function [x,timeinstants,sigma,time]=randomPieceWiseGeneratorMC(Num,c)
N=randi([10000,40000]);
x=zeros(Num,N);
xrand=x;
init=1;
flag=0;
count=1;
timeinstants=zeros(1,N);
n=Num;
C=(ones(n,n)*c)-(eye(n)*c)+eye(n);
[E,L] = eig(C);
while flag<1
initnew=randi([init+1000,init+4000]);
time(count)=initnew;
timeinstants(initnew)=1;
if count==1
sigma(count)=randi([1,10]);
else
r=round(rand(1));
rLow=(rand(1)*0.6)+0.1;
rHigh=(rand(1)*3)+1.5;
sigma(count)=sigma(count-1)*(r*rLow + (1-r)*rHigh);
end
x(:,init:initnew)=sigma(count)*E*sqrt(L)*randn(Num,initnew-init+1);
count=count+1;
if (initnew+4000)>N
flag=1;
timeinstants(time(count-1))=0;
end
init=initnew;
end
x(:,initnew:N)=E*sqrt(L)*randn(Num,N-initnew+1)*sigma(count-1);
time=time(1:end-1);
##Script
%% LMS test script
[x,timeinstants,sigma,time]=randomPieceWiseGeneratorMC(4,0.1);
WinDesired=10;
[state,h]=ChangeDetectionLMSShiftMultipleSensor1(x,250,20,10,0.3);
plot(x')
hold
plot(state*std(x(1,:)))
%% Statistical Change Detector
[x,timeinstants,sigma,time]=randomPieceWiseGeneratorMC(1,0);
[ChangePoint,ActChange]=PeakChangeDetector(x,200,3);
plot(x)
hold
plot(ActChange*std(x),'r') %Detection of change point
plot(ChangePoint*std(x),'g')%Estimate of change point location
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment