/* Gaussian random number sequence generator by FFT        */
/* pseudorandom numbers ---> Guassian pseudorandom numbers */
/*               xr, xi ---> yr, yi                        */
/*                                                         */
/* (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;                            // data length = FFT length (2^20=1048576)
HIST_DIV=100;                             // histragem parameter

// generate random signal without DC offest
xr=grand(1, DATLEN, 'lgi');              // uniform distribution
xi=grand(1, DATLEN, 'lgi');              // uniform distribution
xr=xr-sum(xr)/DATLEN;                    // remove DC offset
xi=xi-sum(xi)/DATLEN;                    // remove DC offset
xr=xr/max(abs(xr));                      // scaling
xi=xi/max(abs(xi));                      // scaling
x=xr+%i*xi;                              // make complex sequence

// FFT 
y=fft(x);

// split real part and imaginary part
yr=real(y);
yi=imag(y);

// save as wave file
//ymax=max(abs([yr, yi]));
ymax=3300;
savewave('outr.wav', 0.9*yr/ymax, 8000);
savewave('outi.wav', 0.9*yi/ymax, 8000);

// plot waveform 
scf(0);
clf;
plot2d(yr(1:1000), axesflag=1, style=2, rect=[0, -1.2*ymax, 1000, ymax]);
plot2d(yi(1:1000), axesflag=1, style=5, rect=[0, -1.2*ymax, 1000, ymax]);

//ymax=1.05*max(abs([yr, yi]));
ymax=3500;

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

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

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

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