{$V-}
PROGRAM dft ;

{This program performs digital Fourier transform spectral analysis
 with option of moving window.

 Ask for :
 Input data file (first column is time record, second one is start of data)
 Input file title
 Value of missing data
 Even or Uneven data spacing
 Column of data to analyze (1 is the start column of data (2nd in the input))

 Remark :
 1. Moving window only applies to EVEN spacing data.
 2. Format of input data file
    1882 1.2 3.45 6.78 .........
    1883 3.5 6.89 34.6 .........
    1884 .......................
    .    .......................
    .
    .
    1985 .......................             }


{$I UTILITY.LIB}

type
    variable = array[1..1500] of real;
var
    out_name,f_in_name,prt_name,
    prt_name1,out_name1             : string[60];
    f_in, f_out,f_prt               : text ;
    f_id,ext                        : string[75];
    xmiss,number,count,nonm,noyear  : integer ;
    ch,eorue                        : char;
    month,column,year,nodata        : integer;
    temdata,sinwork,coswork,
    data,time                       : variable;
    avg,sd,c                        : array[1..3] of real;
    initial,nomiss,odd,ok,refe      : boolean;
    wr_print,wr_file,rangeall       : boolean;
    mask                            : packed array[1..1500] of boolean;
    ibdint,ibdbeg,ibdend,nbands,
    ifreq,iband,ibd,idata           : integer;
    twopi,pi,arg,freq,fs1,fs2       : real;
    amp,pha,rsq,fratio              : real;
    f99a,f99b,f95a,f95b,f90a,
    f90b,factor                     : real;
    leng,lag,outfile,xxi,xxf,i      : integer;


PROCEDURE setph;
begin
     factor := ((column-0.5) + 12 * (xxi - 1900))/12 ;
end;

PROCEDURE window;
{ Moving window }
begin
    write('Window length (integer): ');
    readln(leng);
    write('Time lag (integer): ');
    readln(lag);
end;

PROCEDURE harmonic;
{Harmonic freqs or specify the freq spacing between integral band }
begin
    writeln;
    write('Do you wish to use the Harmonic freqs ...(y/n): ');
    if yorno then
    ibdint := 1       {number of steps between integral bands}
    else
    begin
        write('Then specify # of bands between integral bands (integer) : ');
        readln(ibdint);
    end;
end;

PROCEDURE range;
{-------- entire spectra or some cut-off freqs ---------------}
begin
    writeln;
    write('Do you wish to analyze the entire spectra ...(y/n): ');
    if yorno then
    begin
        ibdend := nbands;
        ibdbeg := 0;
        rangeall := true;
    end
    else
    begin
        write('Specify the low band #  (integer) : ');
        readln(ibdbeg);
        write('Specify the high band # (integer) : ');
        readln(ibdend);
        if ibdbeg < 0 then ibdbeg := 0;
        if ibdend > nbands then ibdend := nbands ;
        rangeall := false;
    end;
end;

PROCEDURE refeph;
begin
     writeln;
     writeln('Reference time for initial phase is first data point.');
     write('Do you want the phase referenced to midnight 1 Jan 1900 (y/n)? ');
     refe := readchar(['Y','N']) = 'Y';
end;

PROCEDURE getmsk ;
{ put data in a logical array, true=non-missing,false=missing }
var
    i                           : integer;
begin
    {initialize mask}
    for i := 1 to nodata do
    mask[i] := true;

    {set missing data = false}
    nonm := 0;  {number of non-missing data}
    for i := 1 to nodata do
    begin
        if (round(data[i]) = xmiss) then
        mask[i] := false
        else
        nonm := nonm +1;
    end;

    if (nonm = nodata) then nomiss := true
    else nomiss := false;

end;

PROCEDURE outprt1;
{-------- output device = printer -------------}
begin
    writeln('Name of printer file(s)');
    writeln('Enter any of DRIVE:\DIRECTORY\FILE NAME (NO EXTENSION)');
    write('Press <ENTER> for no printer file ===> ');
repeat
    readln(prt_name1);
    if length(prt_name1) = 0 then
    begin
         wr_print := false;
         exit;
    end;
    if pos('.',prt_name1) > 0 then
    delete(prt_name1,pos('.',prt_name1),4);
    prt_name := concat(prt_name1,'.',ext);
    ok := not checkfname(prt_name);  {true if file is found }
    if not ok then   { file is found }
    begin
         write('File exists. Do you want to overwrite it (y/n) : ');
         ok := yorno;
    end;
    if ok then   { Want to overwrite it/ file does not exist }
    begin
         assign(f_prt,prt_name);
         {$I-} rewrite(f_prt) {$I+};
         ok := (IOresult = 0);
         if not ok then       { file name is not right }
         writeln('Bad file name, try again.');
         close(f_prt);
    end;
    if not ok then { Bad file name/ does not want to overwrite it}
         write('File name again or <Enter> to quit : ');
until ok;
    wr_print := ok;
end;

PROCEDURE outprt2 ;
{--------- open printer file -------------}
begin
    if wr_print then
    begin
       prt_name := concat(prt_name1,'.',ext);
       assign(f_prt,prt_name);
       rewrite(f_prt);
    end;
end;


PROCEDURE outtrc1;
{--- output device = output file serves as a input file for TSR
        Time Series Reconstruction -------------------------------}
var
   okk    : boolean;
begin
    writeln;
    writeln('Name of output file');
    writeln('Enter any of DRIVE:\DIRECTORY\FILE NAME (NO EXTENSION)');
    write('Press <ENTER> for no output file ===> ');
repeat
 repeat
    readln(out_name1);
    if length(out_name1) = 0 then
    begin
         wr_file := false;
         exit;
    end;
    if pos('.',out_name1) > 0 then
    delete(out_name1,pos('.',out_name1),4);
    out_name := concat(out_name1,'.',ext);
    okk := out_name <> prt_name; { false if name is same as printer output}
    if not okk then
      write('Same file name as printer file, type again or <Enter> to quit : ');
 until okk ;
    ok := not checkfname(out_name);  {true if file is found }
    if not ok then   { file is found }
    begin
         write('File exists. Do you want to overwrite it (y/n) : ');
         ok := yorno;
    end;
    if ok then   { Want to overwrite it/ file does not exist }
    begin
         assign(f_out,out_name);
         {$I-} rewrite(f_out) {$I+};
         ok := (IOresult = 0);
         if not ok then       { file name is not right }
         writeln('Bad file name, try again.');
         close(f_out);
    end;
    if not ok then { Bad file name }
         write('File name again or <Enter> to quit : ');
until ok;
    wr_file := ok;
end;

PROCEDURE outtrc2;
{-------- open output file ----------------}
begin
    if wr_file then
    begin
       out_name := concat(out_name1,'.',ext);
       assign(f_out,out_name);
       rewrite(f_out);
    end;
end;

PROCEDURE getmn (var mn:real; var x:variable);
{ calculate the mean of a data set }
var
    i,num                     : integer;
begin
    mn := 0.0;
    if (not nomiss) then  {there are missing data in data set}
    begin
        for i := 1 to nodata do
        begin
             if mask[i] then
             mn := mn + x[i];
         end;
    num := nonm;
    end;
    if nomiss then  { no missing data in data set }
    begin
        for i := 1 to nodata do
        mn := mn + x[i];
        num := nodata;
    end;
    mn := mn/num;
end;

PROCEDURE getsd(var sdev:real; var x:variable; mean:real);
{ calculate the standard deviation of a data set }
var
    i,num                   : integer;
begin
    sdev := 0.0;
    if (not nomiss) then {there are missing data in data set}
    begin
        for i := 1 to nodata do
        begin
             if (mask[i]) then
             sdev := sdev + sqr((x[i] - mean))
        end;
    num := nonm;
    end;
    if nomiss then   {there are no missing data in data set}
    begin
        for i := 1 to nodata do
        sdev := sdev + sqr((x[i] - mean));
        num := nodata;
    end;
    sdev := sqrt(sdev/(num-1.0))
end;

PROCEDURE getcor(var cor : real; var x1,x2 : variable;
m1,m2,sdev1,sdev2 : real) ;
{ calculate the correlation between two data set }
var
    i, num                     : integer;
begin
    cor := 0.0;
    if (not nomiss) then  {there are missing data in data set}
    begin
        for i := 1 to nodata do
        begin
             if mask[i] then
             cor := cor + (x1[i]-m1)*(x2[i]-m2);
         end;
   num := nonm;
   end;
    if nomiss then  { no missing data in data set }
    begin
        for i := 1 to nodata do
        cor := cor + (x1[i]-m1)*(x2[i]-m2);
        num := nodata;
    end;
    cor := cor/sdev1/sdev2/(num - 1.0);
end;



PROCEDURE nyqist(var ap,ph : real);
begin
    if (ap >= 0.0) then ph := 0.0
    else  ph := 0.5;
    ap := abs(ap);
end;

FUNCTION atan2(ph,aa : real): real;
begin
    if (aa < 0) then
    begin
        if (ph > 0) then atan2 := pi + arctan(ph/aa);
        if (ph < 0) then atan2 := arctan(ph/aa) - pi;
    end
    else   atan2 := arctan(ph/aa);
end;

PROCEDURE period ;
{--------- calculate periodogram ------------------}
var
    a,bandno                           : real;

begin

{-------- initialize -------------}
    twopi := 8.0 * arctan(1.0);
    pi := 4.0 * arctan(1.0);

{-------- set odd = true if the # of points is odd -----------}
    if ((noyear mod 2) = 1) then odd := true
    else odd := false;

    if wr_print then  {----- to the printer --------}
    begin
  write(f_prt,'    Band  Frequency     Period        R^2  Amplitude');
    writeln(f_prt,'      Phase    F-ratio');
  write(f_prt,'    ----  ---------     ------        ---  ---------');
 writeln(f_prt,'      -----    -------');
    end;
    if refe then
    setph;
    {----- loop over integer bands -------}
    for iband := ibdbeg to ibdend do
    begin

       {---- loop over band interval -------}
        for ibd := 1 to ibdint do
        begin

        {--- do not calculate for zero freq and above final band }
        if ((iband > 0) or (ibd >1)) and
           ((iband < ibdend) or (ibd = 1)) then
             begin

             {-- calculate the freq and period -----}
             bandno := iband+(ibd-1.0)/ibdint;
             freq := bandno/noyear;

             {--- calculate the sin and cosin data ---}
             for idata := 1 to nodata do
             begin
                 arg := (time[idata] - time[1] + factor)*freq*twopi;
                 coswork[idata] := cos(arg);
                 sinwork[idata] := sin(arg);
             end;

             {-- calculate means and sd of sin and cosin ------}
             getmn(avg[2],coswork);  {-- mean of cosin --}
             getmn(avg[3],sinwork);  {-- mean of sin ----}
             getsd(sd[2],coswork,avg[2]); {-- sd of cosin--}
             getsd(sd[3],sinwork,avg[3]); {-- sd of sin ---}

             {-- calculate correlations ------------------------}
             getcor(c[1],data,coswork,avg[1],avg[2],sd[1],sd[2]);
        {-- special handling of Nyquist freq ----------}
             if ((iband = nbands) and (not odd)) then
             begin
                 amp := c[1]/sd[2]*sd[1];
                 rsq := c[1]*c[1];
                 fratio := (rsq*(nonm-3)) / ((1-rsq)*2);
                 nyqist(amp,pha);
             end
             else
             begin
             getcor(c[2],data,sinwork,avg[1],avg[3],sd[1],sd[3]);
             getcor(c[3],coswork,sinwork,avg[2],avg[3],sd[2],sd[3]);
             {-- calculate regression coefficients--------}
             a := 1.0-c[3]*c[3];
             amp := (c[1]-c[2]*c[3])/a/sd[2]*sd[1];
             pha := -(c[2]-c[1]*c[3])/a/sd[3]*sd[1];

             {-- calculate amp, pha, and rsq --------------}
             rsq := (c[1]*c[1]+c[2]*c[2]-2.0*c[1]*c[2]*c[3])
                     /a;
             fratio := (rsq*(nonm-3)) / ((1-rsq)*2);

             a := amp;
             amp := sqrt(a*a+pha*pha);
             pha := atan2(pha,a)/twopi;
            end;

             write(bandno:10:5,'(',freq:7:5,') ');

             {------- to a printer ---------}
             if wr_print then
             begin
             write(f_prt,bandno:8:4,freq:11:5,(1.0/freq):11:5);
             writeln(f_prt,rsq:11:5,amp:11:5,pha:11:5,fratio:11:5);
             end;


             {------- to a output file -------}
             if wr_file then
             begin
             write(f_out,bandno:8:4,' ',freq:11,' ',(1.0/freq):11);
writeln(f_out,' ',rsq:11,' ',amp:11,' ',pha:11,' ',fratio:11);
             end;

             end; {-------- end calculation ----------------}
        end;   {--- end loop of ibd (band interval) ---------}

    end;       {--- end loop of iband (integer band) --------}

    if wr_print then close(f_prt);
    if wr_file  then close(f_out);

end;

{----------------- main program ---------------------}
begin
    clrscr;
    gotoxy(1,1);
    writeln('*************************************************');
    writeln('   Digital Fourier Transform Spectral Analysis');
    writeln('*************************************************');
    writeln;

{--------  open input data file for reading ------------------}
    write('Specify the data file containing the time series: ');
repeat
    readln(f_in_name);
    if length(f_in_name) = 0 then
       exit;
    if pos('.',f_in_name) = 0 then   { without typing .PRN }
    f_in_name := concat(f_in_name,'.','PRN');
         ok := checkfname(f_in_name);
         if not ok then
         begin
              writeln('Cannot find file: ',f_in_name);
              write('Retype file name (or <Enter> to quit) : ');
         end;
until ok;
    assign(f_in,f_in_name);
    reset(f_in);

{-------- read in TITLE, MISSING, EVEN or UNEVEN, and COLUMN of data ---}
    write('Title for printer output : ');
    readln(f_id);
    write('Value of missing data : ');
    readln(xmiss);
    write('Even or Uneven data (e/u) : ');
    eorue := readchar(['E','U']);
    column := 0;
    repeat
         write('Which column of data to analyze : ');
         readln(column);
    until column >= 1;


{----- checking input FILE,TITLE,MISSING,EVEN or UNEVEN,and COLUMN of data }
    writeln;
    writeln('****** CHECKING **********************************');
    writeln('INPUT DATA FILE: ',f_in_name);
    writeln('FILE TITLE: ',f_id);
    writeln('MISSING DATA: ',xmiss);
    writeln('EVEN or UNEVEN DATA: ',eorue);
    writeln('COLUMN of DATA to ANALYZE: ',column);
    writeln('**************************************************');
    write('continue... (y/n): ');
    if not yorno then
    begin
        writeln('**** Program is aborted by the user ****');
        halt;
    end;

{----- read in first column of data : TIME (time record)  and
       column of data to analyze : TEMDATA (data ) ---------}
    number := 0 ;
    writeln('Reading data, please wait ...');
    while not eof(f_in) do
    begin
         number := number + 1;
         read(f_in,time[number]);
         for month := 1 to column do
         read(f_in,temdata[number]);
         readln(f_in);
    end;
    close(f_in);

{------- Moving window only applies to EVEN spacing data ------}
  if eorue = 'E' then
  begin
{---------------- Moving Window or Not ---------------}
    writeln;
    writeln('************* Moving window ***********');
    write('Moving window (only applies to EVEN spacing data) ...(y/n): ');
    if yorno then
          begin
             window;
             nodata := leng;
          end
    else
          begin
             nodata := number;
             lag    := number;
          end;
    noyear := nodata;
  end

  else
     begin
         writeln('Total number of years can be different from total number ');
         write('of data points. Enter number of years : ');
         readln(noyear);
         nodata := number;
         lag    := number;
     end;

    nbands := noyear div 2;
    harmonic; { freq spcing }
    range;    { spectra range  }
    refeph;   { set phase reference time at 1900 or first data point
                yes = 1900, no = first data point   }
    factor := 0.0;

    {---- open output file(s) --------------}
    writeln;
    writeln('******** Send result to the output device(s) *********');
    writeln('There are two possible output files:');
    writeln('1. Printer file for printer');
    writeln('2. Output file for Time Series Reconstruction (TSR)');
    writeln;

    outfile := 0;   {-- number of outfile -----}
    count := 0;     {--- lag -------------}
    while (count + nodata) <= number do
    begin
    for i := 1 to nodata  do
        data[i] := temdata[i+count];

    outfile := outfile + 1;

{------- set extension of output file as *.001 to *.999 -------}
    str(outfile,ext);
    if outfile < 10 then ext := concat('00',ext)
    else if
           outfile < 100 then ext := concat('0',ext);
    if outfile = 1 then {--- asking output file name only at first time ----}
    outprt1;
    outprt2;
    if outfile = 1 then
    outtrc1;
    outtrc2;

    xxi := round(time[1]) + count;  {------- initial time step ------}
    xxf := xxi + noyear -1;  {------- final time step --------}
    getmsk; { put data in a logical array }

{-- calculate mean and sd of data -------}
    getmn(avg[1],data);
    getsd(sd[1],data,avg[1]);

{-- write out some output ----------------}
    if wr_print then   {--- to the printer -------}
    begin
    writeln(f_prt);
    writeln(f_prt,'Digital Fourier Transform Spectral Analysis');
    writeln(f_prt,'*******************************************');
    writeln(f_prt);
    writeln(f_prt,'This is the output printer file : ',prt_name);
    writeln(f_prt,'Input data file  : ',f_in_name);
    writeln(f_prt,'Input file   ID  : ',f_id);
    writeln(f_prt,'Initial time step: ',xxi:6);
    writeln(f_prt,'Final   time step: ',xxf:6);
    if eorue = 'E' then
    writeln(f_prt,'Data spacing     : ','Equal')
    else
    writeln(f_prt,'Data spacing     : ','Unequal');
    writeln(f_prt,'Value of missing data: ',xmiss:3);
    writeln(f_prt);
    writeln(f_prt,'Column to be analyzed    : ',column:3);
    writeln(f_prt,'Number of data read      : ',nodata:3);
    writeln(f_prt,'Number of missing data   : ',(nodata-nonm):3);
    writeln(f_prt,'Number of data to be used: ',nonm:3);
    writeln(f_prt);
    writeln(f_prt,'Data   average : ',avg[1]:12:5);
    writeln(f_prt,'Data std. dev. : ',sd[1]:12:5);
    writeln(f_prt); writeln(f_prt);


{-------- Fisher significance limits for Harmonic R^2 value ----}
writeln(f_prt,'Fisher (1929) Significance Limits for Harmonic R^2 Values');
   writeln(f_prt);

   fs1 := nbands;
   fs2 := 1.0/(fs1-1);
   f99a := 1.0 - exp(fs2*ln(0.01/fs1));
   f99b := 1.0 - exp(fs2*ln(0.01));
   f95a := 1.0 - exp(fs2*ln(0.05/fs1));
   f95b := 1.0 - exp(fs2*ln(0.05));
   f90a := 1.0 - exp(fs2*ln(0.1/fs1));
   f90b := 1.0 - exp(fs2*ln(0.1));
writeln(f_prt,'                   a posteriori      a priori');
   writeln(f_prt);
writeln(f_prt,'99% Sig Lvl:',f99a:18:8,f99b:18:8);
writeln(f_prt,'95% Sig Lvl:',f95a:18:8,f95b:18:8);
writeln(f_prt,'90% Sig Lvl:',f90a:18:8,f90b:18:8);
writeln(f_prt,'White Noise:',(1.0/fs1):36:8);
     writeln(f_prt); writeln(f_prt);
    if rangeall then
    writeln(f_prt,'DFT spectral analysis for the entire spectra')
    else
    begin
    writeln(f_prt,'DFT spectral analysis for some cut-off freg:');
    writeln(f_prt,'Low band # : ',ibdbeg:3,'     High band # : ',ibdend:3);
    end;
    writeln(f_prt);
    end;
    if refe then xxi := 1900; {reference point is 1900}
    if wr_file then    {------- to a output file ---------}
    begin
    writeln(f_out,f_id);
    if rangeall then
    write(f_out,xxi:6,' ',(ibdint*(ibdend-ibdbeg)):6,' ',nodata:6)
    else
    write(f_out,xxi:6,' ',(ibdint*(ibdend-ibdbeg+1)):6,' ',nodata:6);
    writeln(f_out,' ',column:3,' ',avg[1]:10,' ',sd[1]:10);
    end;

{-- calculate periodogram, and write result to output device(s) ---}
    writeln;
    writeln(xxi,xxf:10);
    period;
    count := count + lag;
    end;   {end of while }

    clrscr;
    gotoxy(1,10);
    writeln('**************************************************');
    writeln('PROGRAM IS DONE');
    if wr_print then
       if (outfile > 1) then
           writeln('PRINTER FILES ARE ',prt_name1,'.001 - ',prt_name)
       else
           writeln('PRINTER FILE IS ',prt_name);
    if wr_file then
       if (outfile > 1) then
           writeln('OUTPUT FILES  ARE ',out_name1,'.001 - ',out_name)
       else
           writeln('OUTPUT FILE IS ',out_name);
    writeln('**************************************************');
end.