/* Central Limit Theorem (CLT) method */
/* (c) 2025 cepstrum.co.jp            */

clear;

//SEED='a5a5a5a5';
//SEED='5a5a5a5a';
SEED='0f0f0f0f';
//SEED='f0f0f0f0';

grand('setgen', 'mt');              // use Mersenne-Twister 
grand('setsd',  hex2dec(SEED));     // set seed

DATLEN =2^20;                       // random data length (2^20=1048576)
CLTLEN =12;                         // Central Limit Theorem sum length
HIST_DIV=100;                       // histragem parameter

y=zeros(1:DATLEN);

// Central Limite Theorem (CLT) method
mprintf("           |");
for i=1:20
  mprintf("_");
end
mprintf("|\nworking... ");
for i=1:DATLEN
  if modulo(i-1, floor(DATLEN/20))==0
    printf("*");
  end
  y(i)=sum(grand(1, CLTLEN, 'def')-0.5);       
end
printf("*");

// calculate histgram
//ymax=1.05*max(abs(y));
ymax=5.5;
hist_range=2.0*ymax/(2*HIST_DIV+1);
hist=zeros(1:2*HIST_DIV+1);
for i=1:DATLEN
  ix=floor((y(i)+ymax)/hist_range)+1;
  hist(ix)=hist(ix)+1;
end

// plot histgram
scf(1);
//clf;            // clear screen
gcf().axes_size = [800, 600];
plot2d(((1:2*HIST_DIV+1)-HIST_DIV-1)*hist_range, hist, axesflag=1, style=2, rect=[-ymax, 0, ymax, 24000]);
plot2d([0, 0], [0, 24000],style=35);
title fontsize 2
title("cepstrum.co.jp");
xlabel "output histgram (Gaussian)"