
%    kurtfig3.m  -- this script has been tested with GNU Octave 2.1.35 
%    test of the recursive estimation of the kurtosis
%    with a Gaussian noise with changing mean and variance
%
%    ref:  Eric Chassande-Mottin, "Testing the normality of the 
%    gravitational wave data with a lost cost recursive estimate 
%    of the kurtosis", in Proceedings of PSIP'2002, Grenoble, 2002.
%
%
%     Copyright (C) 2002, Eric Chassande-Mottin (ecm@obs-nice.fr)
% 
%     This program is free software; you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation; either version 2 of the License, or
%     (at your option) any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program; if not, write to the Free Software
%     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
%


clear

% -- set plot context

graw('reset;set data style lines;') 
gset term x11
gset nokey
gset lmargin 9;
gset rmargin 1;

debug=0; % 1: x11, 0: postscript

% -- signal parameters

tsig=162; %s
fs=50; %Hz

m1=1; % 10;
m2=-2; % 2*m1;
m3=1; 
s1=1;
s2=1;
s3=4;

% -- window parameters

tot=5; % in %
th=20; %s, duration when decreased to tot% of start amplitude

% --  make signal

N=round(tsig*fs);
t=(0:N-1)/fs;
sig=randn(N,1);
m=[m1*ones(N/3,1);m2*ones(N/3,1);m3*ones(N/3,1)];
sigma=sqrt([s1*ones(N/3,1);s2*ones(N/3,1);s3*ones(N/3,1)]);
x=sigma.*sig+m;
xmax=1.1*max(abs(x));

% -- make window

% mean
a1=(100/tot)^(-1/(th*fs));
C1=1-a1;
% variance
a2=a1;
C2=(1+a1)*(1-a2)/2;
% 4th order moment
a4=a1;
% kurtosis
k4=2*(1-a1+2*a1^2)/(1+a1)/(1+a1^2);
k2=6*(1-a1)*(1+a1+2*a1^2)/(1+a1)^2/(1+a1^2);
C4=(1-a4)/k4;
%slope=C4*k4/(1-a4)
bias=-k2/k4

% -- short term estimates

mu1=zeros(N,1);
mu2=zeros(N,1);
kappa4b=zeros(N,1);
kappa4=zeros(N,1);

% init values

mu1(1)=0;
mu2(1)=1;
kappa4b(1)=0;
kappa4(1)=0;

%other possibilities:
%m=(x(1)+x(2))/2;
%mu1(1)=m;
%mu2(1)=((x(1)-m)^2+(x(2)-m)^2);

% recursive loops

for j=2:N,
  mu1(j)=a1*mu1(j-1)+C1*x(j);
  d2=(x(j)-mu1(j-1))^2;
  mu2(j)=a2*mu2(j-1)+C2*d2;
  d2=d2/mu2(j-1);
  kappa4b(j)=(1+C1-2*C1*d2)*kappa4b(j-1)+C1*d2^2;
  kappa4(j)=kappa4b(j)+bias;
end;
kmax=1.1*max(kappa4);

% -- make plots

if debug==0
  gset term posts eps enhanced color solid "helvetica" 20
  gset output "../figure/kurtfig3.eps";
end;  

% -- plotting...
multiplot(2,1);
idx=1:10:N;

subplot(211);
grid("on");
ylabel('signal');
eval(sprintf('gset xrange [0:%f];',t(N)));
eval(sprintf('gset yrange [%f:%f];',-xmax,xmax));
gset ytics -12,4,12
plot([t(N/3) t(N/3)],[-xmax xmax],'b');
hold on
plot([t(2*N/3) t(2*N/3)],[-xmax xmax],'b');
plot(t(idx),x(idx),'r');
hold off

% relax
gset xtics autofreq
gset ytics autofreq
gset autoscale x
gset autoscale y

subplot(212);
xlabel('time (s)');
ylabel('kurtosis');
eval(sprintf('gset xrange [0:%f];',t(N)));
eval(sprintf('gset yrange [2:%f];',kmax));
grid("on");
plot(t(idx),kappa4(idx));
hold on
%eval(sprintf("plot([t(N/3) t(N/3)],[2 %f],'b--')",kmax));
%eval(sprintf("plot([t(2*N/3) t(2*N/3)],[2 %f],'b--')",kmax));
plot([t(N/3) t(N/3)],[2 kmax],'b--');
plot([t(2*N/3) t(2*N/3)],[2 kmax],'b--');
hold off


% relax
gset xtics autofreq
gset ytics autofreq
gset autoscale x
gset autoscale y

% -- relax
oneplot();
gset autoscale x
gset autoscale y
gset term x11











