Compare KLT and kurtosis

global N = 2^21;  # Number of samples
R = 1/30;  # Signal-to-noise ratio
P = 5;     # Number of runs

function klt( d )
  y = autocor(d, 511);
  z = toeplitz(y);
  [vctr,val] = eig(z);
  plot(val(find(val))(512:-1:480))
endfunction

function r = kurt( d )
  global N;
  m = d-mean(d);
  for k = 1:N
    a(k) = m(k)^2;
    b(k) = a(k)^2;
  endfor
  r = ((sum(b)/N)/((sum(a)/N)^2)) - 3;
endfunction

s = sinewave(N,10)*R;
for k = 0:P-1
  n = randn(1, N);
  x = s + n;

  subplot(P,2,2*k+1);
  klt(n);
  title(sprintf("Noise (Kurtosis: %f)",kurt(n)))

  subplot(P,2,2*k+2);
  klt(x);
  title(sprintf("Signal plus noise (Kurtosis: %f)",kurt(x)))
endfor