//-------------------------------------------------------------------//

// Synopsis:   Compute the short-time periodic windowed FFT of a signal.

// Syntax:     spectrogram ( SEQ, FRAMESIZE, NOVERLAP )
//             spectrogram ( SEQ, FRAMESIZE, NOVERLAP, ZEROPAD )

// Description:

//      For each frame, window it with a Hanning window, compute the
//      FFT, remove the negative frequencies, and take the log of
//      the resulting elements.  Then hop ahead by (FRAMESIZE -
//      NOVERLAP) elements to the next frame. 

//      If the ZEROPAD argument is supplied, that many zeros are
//      appended to the frame before computing the FFT.
//      ZEROPAD + FRAMESIZE should be a power of two.

//      The resulting array has (FRAMESIZE+ZEROPAD)/2 rows and a
//      number of columns determined by the framesize and overlap
//      values. Low frequencies come first, i.e. they have low row
//      indices. 

// Sample seq:     seq = readb("somedatafile.rb");
//                 spect = spectrogram(seq, 512, 0);
//                 plmesh(<< x=1:spect.nr; y=1:spect.nc; z=spect >>);

// This file is a translation of spectrogram.m from the Osprey toolbox.
// The results of this function may differ from Matlab's slightly 
// since the Hanning window is computed differently.

require window

//-------------------------------------------------------------------//

spectrogram = function (seq, frameSize, nOverlap, zeroPad)
{
  local (seq, frameSize, nOverlap, zeroPad)

  if (!exist (zeroPad)) { zeroPad = 0; }

  inrows = seq.nr;
  incols = seq.nc;

  if (inrows != 1 && incols != 1)
  {
    error ("Spectrogram requires a 1-dimensional sequence");
  }

  nfft    = frameSize + zeroPad;
  seqsize = max (inrows, incols);
  outcols = floor(1 + (seqsize - frameSize) / (frameSize - nOverlap));
  outrows = nfft / 2;
  result  = zeros (outrows, outcols);

  // column vector
  win  = window (frameSize, "hann");

  for (i in 1:outcols)
  {
    start = 1 + (frameSize - nOverlap) * (i - 1);
    frame = reshape (seq[start:(start + frameSize - 1)], frameSize, 1);

    // column vector
    spectrum = fft (frame .* win, nfft); 
    result[;i] = spectrum [1:outrows];
  }

  // move inside loop if memory problems
  result = log (abs (result));

  return result;
};