Compare KLT and kurtosis

global N = 2^20;  # Number of samples
R = 1/30;  # Signal-to-noise ratio
P = 2^10;  # Signal period
Q = 5;     # Number of runs
X = 2^13;  # Block size

function klt( d )
  disp("Calculating autocorrelation...");
  fflush(stdout);
  y = autocor(d, 511);
  disp("Calculating Toeplitz matrix...");
  fflush(stdout);
  z = toeplitz(y);
  disp("Calculating eigenvalues...");
  fflush(stdout);
  [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

printf("Sample count: %d\n", N);
printf("Block size: %d\n", X);
printf("SNR: %f\n", R);
printf("Signal period: %d\n", P);
disp("");
title(sprintf("%d samples, SNR=%f", N, R));

fd = fopen("random.bin");
n = 0;
s = sinewave(P,P)*R;
x = 0;
for k = 0:Q-1
  printf("Run: %d/%d\n", k+1, Q);
  fflush(stdout);
  for i = 1:N
   ii = mod(i-1,X)+1;
   if (ii == 1)
      b = fread(fd, X, "int32");
      msg = sprintf("Loading block: %d/%d...", (i-1)/X+1, N/X);
      printf(strcat(msg,strrep(blanks(length(msg))," ","\b"))); #Print message and back up cursor
      fflush(stdout);
    endif
    n(i) = b(ii)/2^31; #Convert b to value between 0 and 1
    x(i) = s(mod(i-1,P)+1) + n(i);
  endfor

  disp("\nProcessing 'Noise'...");
  fflush(stdout);
  subplot(Q,2,2*k+1);
  klt(n);
  disp("Calculating kurtosis...");
  title(sprintf("Noise (Kurtosis: %f)",kurtosis(n)))

  disp("Processing 'Signal plus noise'...");
  fflush(stdout);
  subplot(Q,2,2*k+2);
  klt(x);
  disp("Calculating kurtosis...");
  title(sprintf("Signal plus noise (Kurtosis: %f)",kurtosis(x)))
endfor