/* plot histgram           */
/* (c) 2025 cepstrum.co.jp */

clear;

FILE_NAME='wt_out.txt';     // input data file name
DATLEN   =2^20;             // 2^20=1048576
HISTLEN  =300;              // histgram plot parameter
XYLEN    =36000;            // X-Y plot data length

x=zeros(DATLEN);
y=zeros(DATLEN);         

// read data file from c_wt_gaussian.c/c_wt_gaussian.exe
fd=mopen(FILE_NAME, 'rt');
for i=1:DATLEN   
  [n, x(i), y(i)]=mfscanf(fd, '%e %e');
end
mclose(fd);

// plot histgram (x)
scf(11);
clf;
gcf().axes_size = [600, 600];
histplot(40, x, axesflag=1, normalization=%f, rect=[-1.1, 0, 1.1, 40000], style=2);
title fontsize 2
title("cepstrum.co.jp");
xlabel "input (uniform)"

// plot histgram (y)
scf(21);
clf;
gcf().axes_size = [600, 600];
ymax=max(abs(y));
histplot(40, y, axesflag=1, normalization=%f, rect=[-1.05*ymax, 0, 1.05*ymax, 120000], style=2);
title fontsize 2
title("cepstrum.co.jp");
xlabel "output (Gaussian)"

// plot Gaussian histgram curve
gaus=101000.0*exp(-0.00013*(-HISTLEN:HISTLEN)^2);
scf(21);
plot2d(ymax*(-HISTLEN:HISTLEN)/HISTLEN, gaus, style=13);

// X-Y plot
scf(31);
clf;
gcf().axes_size = [600, 600];
xyrange=max(abs(y(1:2*XYLEN)));
plot2d([0, 0], [-1.05*xyrange, 1.05*xyrange], style=1);
plot2d([-1.05*xyrange, 1.05*xyrange], [0, 0], style=1);
plot2d(y(1:XYLEN), y(XYLEN+1:2*XYLEN), axesflag=1, rect=[-1.05*xyrange, -1.05*xyrange, 1.05*xyrange, 1.05*xyrange], style=0);