function [out]=rbd(sig,win,M1,M2)

%

% Recursive Bayesian Autoregressive Changepoint Detector

% Usage: [out] = rbd(sig,win,M1,M2)

%

% Input: sig ............ 1-D data vector

% win ..... length of sliding window

% (default 400), min~200, max~5000

% M1, M2 .......... order of AR model (from left and right)

% (default 6)

% Output: out ................~(d_cep1).^2

% is proportional to the signal changes

%

% by Roman Cmejla, cmejla@fel.cvut.cz

% 2011Aug15, ver.3.00

if (nargin == 0), disp('Usage: [out] = rbd(sig, win, M1, M2)'); return; end

if (nargin == 1), M1=6; M2=6; win=400; end;

if (nargin == 2), M1=6; M2=6; end;

if (nargin == 3), M1=M2; end;

% Initialization for m=length(window)/2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

bay(length(sig))=0;

evid(length(sig))=0;

ME=round((M1+M2)/2);

sig=sig(:); sig=sig/max(abs(sig));

m=fix(win/2);

data=sig(1:win);

G(length(data-m),M1+M2)=0;

N=length(data);

% Initial matrix of basis function

for j=2:M1

if j <= m, G(j,:)=[data(j-1:-1:1)' zeros(1,M1-j+1) zeros(1,M2)]; end;

if j > m, G(j,:)=[zeros(1,M1) data(j-1:-1:1)' zeros(1,M2-j+1)]; end;

end;

for j=(M1+1):N

if j <= m, G(j,:)=[data(j-1:-1:j-M1)' zeros(1,M2)]; end;

if j > m, G(j,:)=[zeros(1,M1) data(j-1:-1:j-M2)']; end;

end;

D=data'*data; % Signal energy

CHI=data'*G; % Correlation vector

FI=inv(G'*G); % Inverse correlation matrix

% Initial matrix of Bayesian evidence for estimate normalization

for j=2:ME

G_E(j,:)=[data(j-1:-1:1)' zeros(1,ME-j+1)];

end;

for j=(ME+1):N

G_E(j,:)=[data(j-1:-1:j-ME)'];

end;

CHI_E = data'*G_E; % Correlation vector of BE

FI_E = inv(G_E'*G_E); % Inverse correlation matrix of BE

for mm=m+1:length(sig)-m

% Recursive changepoint estimation for all next samples

% Observe new sample

d_last=sig(mm+m); % new signal sample

G2=[zeros(1,M1) sig(mm+m-1:-1:mm+m-M2)'];

D=D+d_last*d_last;

CHI=CHI+d_last*G2;

W=FI*G2';

LAMBDA=(1+G2*W);

FI=FI-W*inv(LAMBDA)*W';

% Removing old sample

D=D-sig(mm-m)*sig(mm-m);

if mm-m > M1, Z=[sig(mm-1-m:-1:mm-M1-m)' zeros(1,M2)]; else

Z=[sig(mm-1-m:-1:1)' zeros(1,M2) zeros(1,M1-mm+m+1)]; end;

CHI=CHI - sig(mm-m)*Z;

W=FI*Z';

LAMBDA=(1-Z*W);

FI=FI+(1/LAMBDA)*W*W';

% Position update

R=[zeros(1,M1) sig(mm-1:-1:mm-M2)']; % replacing with row of zeros

CHI=CHI - sig(mm)*R;

W=FI*R';

LAMBDA=(1-R*W);

FI=FI+(1/LAMBDA)*W*W';

Q=[sig(mm-1:-1:mm-M1)' zeros(1,M2)]; % replacing with new data

CHI=CHI + sig(mm)*Q;

W=FI*Q';

LAMBDA=(1+Q*W);

FI=FI-(1/LAMBDA)*W*W';

bay(mm)=log(D-CHI*FI*CHI');

% Bayesian evidence update

% Adding new sample

G2_E=[sig(mm+m-1:-1:mm+m-ME)'];

CHI_E=CHI_E+d_last*G2_E;

W_E=FI_E*G2_E';

LAMBDA_E=(1+G2_E*W_E);

FI_E=FI_E-W_E*inv(LAMBDA_E)*W_E';

% Removing old sample

if mm-m > ME, Z_E=[sig(mm-1-m:-1:mm-ME-m)']; else

Z_E=[sig(mm-1-m:-1:1)' zeros(1,ME-mm+m+1)]; end;

CHI_E=CHI_E - sig(mm-m)*Z_E;

W_E=FI_E*Z_E';

LAMBDA_E=(1-Z_E*W_E);

FI_E=FI_E+(1/LAMBDA_E)*W_E*W_E';

FTF=CHI_E*FI_E*CHI_E';

evid(mm)=log(D-FTF);

% Return for position m+1

end;

% Calculate changepoint position probabilities

out=real(evid)-real(bay);

out=[zeros(1,m) 4*out(round(win/2)+1:end)];