/*****************************************************************************
*                                                                            *
*   SPECGRAM.SPL   Copyright (C) 1997-2021 DSP Development Corporation       *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Calculates and plots the magnitude of the STFT              *
*                                                                            *
*   Revisions:    9 Jul 1997  RRR  Creation                                  *
*                 9 Jul 1999  RRR  added fftlen                              *
*                18 Oct 2021  RRR  faster welchseg                           *
*                                                                            *
*****************************************************************************/

/* spectrogram windows */
#define HAMMWND 0


#if @HELP_SPECGRAM

    SPECGRAM

    Purpose: Calculates the 2D Spectrogram as an image.

    Syntax:  SPECGRAM(series, len, lap, fftlen, swin, ampflag)

                 series - Series, the input series.

                    len - Optional. An integer, the segment length, 
                          defaults to 512.

                    lap - Optional, An integer, the overlap length, 
                          defaults to len / 2

                 fftlen - Optional, An integer, the FFT length, 
                          defaults to len.

                   swin - Optional. An integer, the windowing function, 
                          defaults to 0, HAMMING window.

                ampflag - Optional. An integer, the window correction method:

                                0: do not correct amplitude (default)
                                1: correct amplitude
                                2: correct RMS amplitude
                                3: correct mean squared

    Returns: An array, the Amplitude values in Frequency vs Time format.

    Example:

             W1: {gsin(1000,.001,200),gsin(1000,.001,400),gsin(1000,.001,300)}
             W2: Spectrum(W1)
             W2: Specgram(W1, 128);rainbow

             W1 consists of three concatenated sinusoids of 200, 400 and 300
             Hertz.

             W2 shows the frequency spectrum with peaks at 200, 400, and 300
             Hertz.

             W3 divides W1 into columns of 128 points that overlap by 64
             points. The Spectrum (i.e. magnitude of the FFT) of each
             column is calculated and the result is displayed as an
             image with the rainbow colormap.

             The image in W3 shows how the frequency of the series in
             W1 changes over time, clearly showing the distinct 200,
             400 and 300 Hertz components and the times when the components
             occurred.


             W3: Specgram(W1, 128, 64, 1024)

             Same as above, but a 1024 point zero padded FFT is used.

     Remarks:
             See WINFUNC for a list of available windowing functions.

             The SPECTRUM function displays the frequency content of the
             data where the SPECGRAM functions displays the frequencies and
             the times at which those frequency occurs.

             For example, when applied to music, the spectrum (or FFT) only
             indicates the notes and the amplitudes of the notes of a
             given song. The specgram is more like a musical score, displaying
             the notes, amplitudes of the notes and the times at which those
             notes were played.

     See Also:
             FFT
             Ravel
             Sonogram
             Spectrum
#endif



/* Magnitude of Short Term Fourier Transform as Density plot */
SERIES specgram2(s, len, lap, fftlen, swin, fixamp)
{
        local f;

        /* display help if no input */
        if (argc < 1)
        {
                eval("help('specgram')");
                return;
        }

        if (argc < 6)
        {
                if (argc < 5)
                {
                        if (argc < 4)
                        {
                                if (argc < 3)
                                {
                                        if (argc < 2)
                                        {
                                                if (argc < 1) error("specgram - input series required");

                                                len = 512;
                                        }

                                        /* default overlap to half */
                                        lap = int(len / 2);
                                }

                                fftlen = len;
                        }

                        /* default to Hamming window */
                        swin = HAMMWND;
                }

                fixamp = 0;
        }

        /* calculate Spectrum on each column of length len */
        f = welchseg(s, winfunc(swin, len, fixamp), lap, fftlen, rate(s), 11)';

        /* units */
        setvunits(f, getvunits(s), -1);

        /* restore x offset */
        setxoffset(f, xoffset(s));

        /* set to density plot */
        setplotstyle(f, 0);
        setplottype(f, 3);

        if (outargc == 0)
        {
                /* stick in current window */
                f;

                /* set current window to density plot */
                setplotstyle(0);
                setplottype(3);
        }
        else
        {
                return(f);
        }
}