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


#if @HELP_IDCT

    IDCT

    Purpose: Calculates the Inverse Discrete Cosine Transform

    Syntax:  IDCT(s, n)

              s - input series or array

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

    Returns: A series or array.

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

             returns a 20 Hz cosine wave.

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

    See Also:
             Dct
             Dct2
             Fft
             Idct2
             Ifft

    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 inverse discrete cosine transform */
ITERATE idct(b, n)
{
        local bb, y, yy, ww, a;

        if (argc < 2)
        {
                n  = length(b);
                bb = refseries(b);
        }
        else
        {
                bb = extract(b, 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 */
                yy = {ww * bb, 0, -i * ww[2..n] * reverse(bb[2..n])};

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

                y = ifft(yy);

                /* extract inverse DCT */
                a = extract(y, 1, n);
        }

        else
        {
                /*
                 *  even case
                 */
                 
                /* compute precorrection factor */
                ww[1] /= 2;

                /* weights */                
                yy = ww * bb;

                /* units */
                sethunits(yy, gethunits(bb));

                /* compute x tilde using equation (5.93) in Jain */
                y = ifft(yy);

                /*
                 *  reorder elements of each column according to equations (5.93)
                 *  and (5.94) in Jain
                 */
                 
                a = merge(extract(y, 1, int(n / 2)), rev(extract(y, int(n / 2 + 1), n / 2)));
        }

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

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

        return(a);
}