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

clear;

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

yr=zeros(DATLEN);
yi=zeros(DATLEN);         

// read data file from c_ft_gaussian.c/c_ft_gaussian.exe
fd=mopen(FILE_NAME, 'rt');
for i=1:DATLEN   
  [n, yr(i), yi(i)]=mfscanf(fd, '%f %f');
end
mclose(fd);

ymax=max(abs([yr, yi]));      // calculate maximum value

// plot histgram (yr)
scf(21);
clf;
gcf().axes_size = [600, 600];
histplot(40, yr, 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) real part"

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

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

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