#!/bin/sh

if [ $# != 4 ]; then
  echo "usage: fft <input file> <column> <1st line> <section length>"
  echo "The <column>, ignoring the <1st line> first lines, will be read " \
    "from the <input file>, which must contain tabular numeric data in ASCII." \
    "It will be broken in sections with length <section length>.  Each " \
    "section will be Fourier-transformed and written to one column of the " \
    "output file <input file>.fft. The fist column of" \
    "the input file is assumed to contain equidistant times" \
    "in seconds.  From the times in the first and last line used, the" \
    "frequency interval will be derived, and equidistant frequencies written" \
    "to the first column of the output file." | fmt
  echo "<column> is one-based, but <1st line> is zero-based (it is the number" \
    "of lines to skip).  All of the input file must be made up of numerical" \
    "data, even the lines at the top which are ignored.  The <section length>" \
    "should be smaller than the number of lines in the input file (after" \
    "subtracting <1st line>).  Otherwise, or if it is 0, the whole column" \
    "in the input file is used (again, after <1st line>).  If <section length>"\
    "is within range but does not divide the number of available lines," \
    "the surplus lines at the end are left over." | fmt
  exit -1;
fi

outfile="$1.fft"

octave <<EOF
filedata= load( "$1" );
seclen= $4;
sz= size(filedata)(1) - $3;
if (seclen > 0 && sz > seclen)
  sz= seclen*floor(sz/seclen);
else
  seclen= sz;
endif
freqstep= 1/((filedata($3+sz, 1) - filedata($3+1, 1))/(sz-1)*sz);
timedata= filedata($3+1:$3+sz, $2);
timedata= reshape(timedata, [seclen sz/seclen]);
freqdata= fft(timedata);
freqsz= floor(seclen/2)+1;
freqdata= freqdata(1:freqsz, :);     # ignore redundant complex-conjugate half
freqdata= 2*abs(freqdata);	# account for deleted half and take modulus
freqdata(1, :) /= 2;		# undo duplication for constant coefficient
if ( mod(seclen, 2) == 0 )
  freqdata(freqsz, :) /= 2;	# and for Nyquist coefficient
endif
freqdata= [ freqstep*((1:freqsz)-1)' freqdata ];
save -ascii $outfile freqdata
EOF