/*****************************************************************************
*                                                                            *
*   DCT.SPL      Copyright (C) 1998-2012 DSP Development Corporation         *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Calculates the Discrete Cosine Transform                    *
*                                                                            *
*   Revisions:   16 Mar 1998  RRR  Creation                                  *
*                25 Apr 2012  RRR  scaling                                   *
*                                                                            *
*****************************************************************************/


#if @HELP_DCT

    DCT

    Purpose: Calculates the Discrete Cosine Transform

    Syntax:  DCT(s, n)

              s - input series or array

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

    Returns: A series or array

    Example:
             dct(gcos(100, 1/100, 20))

             returns a series with a peak at 20 Hz.

    Example:
             dct(gcos(100, 1/100, 20), 1024)

             Same as above, but the input is zero padded to length 1024
             before the DCT is calculated.

    Remarks:
             The transform is applied to each column if the input is an
             array.  The DCT is often used in image processing
             applications to perform image compression.

    See Also:
             Dct2
             Dst
             Fft
             Idct
             Idct2

    References:

    [1] Jae S. Lim, "Two-dimensional Signal and Image Processing",
        pp. 148-162.  Implements an even-symmetrical DCT.

    [2] Jain, "Fundamentals of Digital Image Processing", pp. 150-153.

    [3] Wallace, "The JPEG Still Picture Compression Standard",
        Communications of the ACM, April 1991.

#endif


/* compute the discrete cosine transform */
ITERATE dct(a, n)
{
        local aa, y, yy, ww, b;

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

        /* weights */
        ww     = exp(-i * (0..n - 1) * pi / (2 * n)) / sqrt(2 * n);
        ww[1] /= sqrt(2);

        if (n % 2 == 1)
        {
                /*
                 *  odd case
                 */
                 
                /* even symmetric series */
                y = {aa, reverse(aa)};

                yy = fft(y);

                /* compute DCT coefficients */
                b = ww * yy[1..n];
        }
        else
        {
                /*
                 *  even case
                 */

                /* re-order the odd and even elements */
                y = {aa[1..2..n], aa[n..-2..2]} ;

                yy = fft(y);

                /* Compute DCT using equation (5.92) in Jain */
                b = 2 * ww * yy;
        }

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

        /* set correct deltax */
        setdeltax(b, 1 / (2 * n * deltax(a)));

        /* X units */
        sethunits(b, gethunits(yy));

        return(b);
}