
%
%    kurtfig1.m  -- this script has been tested with GNU Octave 2.1.35 
%    computes the PDF of the kurtosis estimator and the rate of
%    triggered frames for a Gaussian noise
%
%    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=600; %s
fs=50; %Hz
N=round(tsig*fs);
Nreal=50;

m=0;
s=1;

% -- window parameters

tot=5; % unit is "%"
th=20; % s, duration when decreased to tot% of start amplitude
Nh=fs*th;

% -- make window

% mean
a1=(100/tot)^(-1/(th*fs));
C1=1-a1;
% variance
a2=a1;
C2=(1+a1)*(1-a2)/2;
% 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);
bias=-k2/k4

% -- simus parameters
  
% progress index
step=10;
nstep=floor(N/step);

% histogram and threshold
Nt=25;
k_threshold=linspace(2.4,4,Nt);
h_real=zeros(Nreal,Nt);
pfa_real=zeros(Nreal,Nt);

for n=1:Nreal
  fprintf('noise realization #%d/%d\n',n,Nreal);

% --  make signal

  sig=randn(N,1);
  x=s*sig+m;
  
% -- short term estimates

% init values

  mu1p=0;
  mu2p=1;
  kappa4bp=0;
  kappa4=zeros(N,1);
  
  %other possibilities:
  %m=(x(1)+x(2))/2;
  %mu1(1)=m;
  %mu2(1)=((x(1)-m)^2+(x(2)-m)^2);
  
% recursive loops

  fprintf('recursive kurtosis: ');
  for j=2:N,
    if (rem(j,nstep)==0),
      fprintf('%g ',step*fix(j/nstep)); fflush(1);
    end;
    mu1=a1*mu1p+C1*x(j);
    d2=(x(j)-mu1p)^2;
    mu2=a2*mu2p+C2*d2;
    d2=d2/mu2p;
    kappa4b=(1+C1-2*C1*d2)*kappa4bp+C1*d2^2;
    kappa4(j)=kappa4b+bias;
    mu1p=mu1;
    mu2p=mu2;
    kappa4bp=kappa4b;
  end;
  fprintf('100 %% complete.\n'); fflush(1);
  kmax=1.1*max(kappa4);

  h=hist(kappa4(Nh:N),k_threshold,1/(k_threshold(2)-k_threshold(1)));
  h_real(n,:)=h;

  pfa=zeros(Nt,1);
  kappa4=reshape(kappa4,fs,tsig);
  kappa4max=max(kappa4);
  for j=1:Nt
    pfa(j)=sum(kappa4max(th:tsig)>k_threshold(j))/(tsig-th+1);
  end;
  pfa_real(n,:)=pfa';

end;

% -- compute mean and error bars

h_mean=mean(h_real);
h_std=std(h_real);
hmax=max(max(h_real));
pfa_mean=mean(pfa_real);
pfa_std=std(pfa_real);

% -- make plots

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

multiplot(2,1);

subplot(2,1,1);
gset xrange [2.4:4];
eval(sprintf('gset yrange [0:%f];',hmax));
errorbar(k_threshold,h_mean,h_std);
grid("on");
xlabel('kurtosis');
ylabel('PDF counts');
hold on;
plot(k_threshold,h_mean);
hold off;

subplot(2,1,2);
gset xrange [2.4:4];
gset yrange [0:100];
gset autoscale y
errorbar(k_threshold,100*pfa_mean,100*pfa_std);
xlabel('threshold eta');
ylabel('% of trig. frames');
grid("on");
hold on;
semilogy(k_threshold,100*pfa_mean);
hold off

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




