%Exa = 1;

switch Exa
case 1
%% read and display raw data
Readm4
%% have a very brief look at source

case 2  %% determine the average over a short intervall
  y2 = abs(y-mean(y));  % subtract average value
  FilterLength = 100+1; % average over some points
  t0 = cputime();
  y3 = zeros(Nmax-FilterLength,1);
  for k = 1:Nmax-FilterLength
    y3(k) = mean(y2(k:k+FilterLength));
  end%for
  timing = cputime-t0
  figure(2); plot(y3); 
  xlabel('time t'); ylabel('signal'); grid on

case 3  %% use a convolution command for faster computations
  y2 = abs(y-mean(y));  % subtract average value
  FilterLength = 100+1; % average over some points
  Filter = ones(FilterLength,1)/FilterLength;
  y3 = conv(Filter,y2(1:end-FilterLength+1))';
  y3 = y3(FilterLength:end);
  N = length(y3);
  t = linspace(0,(N-1)*dt,N)';
  figure(2); plot(t,y3);
  xlabel('time t'); ylabel('signal filtered'); grid on

case 4  %% cut above and below, apply logarithm
  amplitude = mean(y3(1:200));
  topcut = 0.8; lowcut = 0.4;  %% choose the cut levels on top and bottom
  Nlow  = find(y3<topcut*amplitude,1);
  Nhigh = find(y3<lowcut*amplitude,1);

  y4 = log(y3(Nlow:Nhigh));
  N  = length(y4);
  t  = linspace(0,(N-1)*dt,N)';
  plot(t,y4)
  xlabel('time t'); ylabel('log(amplitude)'); grid on
  grid on

case 5 %% apply linear regression
  %% run the linear regression
  F = ones(N,2); F(:,2) = t;
  [p,y_var,r,p_var] = LinearRegression(F,y4');
  yFit = F*p;

  figure(3);
  plot(t,y4,t,yFit)  % display the real data and the regression line
  xlabel('time t'); xlabel('log(amplitude)'); grid on

  fprintf('Slope: %g, StdDev of slope: %g,\n   relative error:%g\n',...
          p(2), sqrt(p_var(2)), -sqrt(p_var(2))/p(2));
  factor = 1/p(2);

  figure(4)
  hist(y4'-yFit)
  xlabel('deviation of amplitude'); ylabel('number of events'); grid on 
end%switch
