/*****************************************************************************
*                                                                            *
*   DST.SPL      Copyright (C) 2012 DSP Development Corporation              *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Calculates the Discrete Sine Transform                      *
*                                                                            *
*   Revisions:   26 Apr 2012  RRR  Creation                                  *
*                                                                            *
*****************************************************************************/


#if @HELP_DST

    DST

    Purpose: Calculates the Discrete Sine Transform

    Syntax:  DST(s, n)

              s - input series or array

              n - optional integer, transform length (defaults to length
                  of input)

    Returns: A series or array.

    Example:
             idst(dst(gsin(1000, 1/1000, 20)))

             returns 1000 samples of a 20 Hz sinewave.

    Example:
             dst(gsin(1000, 1/1000, 20), 8192)

             Transforms a sinewave, where the input is zero padded to
             length 8192 before the DST is calculated.

    Remarks:
             The transform is applied to each column if the input is an
             array.

    See Also:
             Dct
             Fft
             Idct
             Idst
#endif


/* compute the discrete sine transform */
ITERATE dst(a, n)
{
        local aa, y, yy, b;

        if (argc < 2)
        {
                n  = length(a);
                aa = refseries(a);
        }
        else
        {
                aa = extract(a, 1, n);
        }

        /* form transformation series */
        y = {0.0, aa, 0.0, reverse(-aa)};

        /* modified FFT */
        yy = fft(y);

        /* compute DST coefficients */
        b = i * yy[2..(n + 1)] / 2;

        if (isreal(a))
        {
                b = real(b);
        }

        /* set correct deltax and offset */
        setdeltax(b, 1 / (2 * n * deltax(a)));
        setxoffset(b, 0.0);

        return(b);
}