View Raw SPL
/*****************************************************************************
*                                                                            *
*   ICCEPS.SPL     Copyright (C) 1999 DSP Development Corporation            *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:       Randy Race                                                 *
*                                                                            *
*   Synopsis:     Inverse complex cepstrum calculation                       *
*                                                                            *
*   Revisions:     9 Jul 1999  RRR  Creation - from CEPSTRUM.MAC             *
*                                                                            *
*****************************************************************************/

#if @HELP_ICCEPS

    ICCEPS

    Purpose: Calculates the inverse complex cepstrum

    Syntax:  ICCEPS(s, d)

             s   - input series or array

             d   - an optional integer, lag value for phase correction,
                   defaults to 0.


    Returns: A real series or array


    Example:
             W1: gtri(100, 1, 1/100)^3
             W2: w1-delay(w1, 60)/2
             W3: w1+w2
             W4: cceps(w3)
             W5: icceps(w4)

             A a synthesized echo at 60 seconds is added to the data of
             W1. The cepstrum is calculated in W4 and the inverse
             cepstrum in W5. The inverse cepstrum is identical to the
             original data except for a period shift of 50 samples;

    Example:

             W6: (c, d) = cceps(w3);c
             W7: icceps(w6, d)

             Same as W5 except the 50 sample period shift is now
             corrected.

    Remarks:

            The complex cepstrum of a series is essentially
            IDFT(log(DFT(s))).  Because the log is used, some
            information is lost and ICCEPS cannot always reconstruct
            the original data. For more information, see CCEPS.


    See Also:
             Cceps
             Rceps


    References:

            [1] Oppenheim & Shafer
                Discrete-Time Signal Processing
                Prentice Hall, 1989
                pp 788-792

            [2] IEEE Press
                Programs for Digital Signal Processing
                IEEE Press, 1979
                Section 7
#endif


/* inverse complex cepstrum */
icceps(s, d)
{
        local f, c;

        if (argc < 2)
        {
                if (argc < 1) error("icceps - input series required");
                
                d = zeros(1, numcols(s));
        }
        
        /* make d a series */
        d = {d};

        f = fft(s);

        setdeltax(f, 1.0);

        /* generate phase correction */
        r = gline(length(s), 1.0, 1 / (length(s)), 0);

        /* use *^ in case we have multi column data */
        p = 2 * pi * r *^ d;

        /* inverse cepstrum */
        c = real(ifft(exp(real(f) + i * (imag(f) + p))));

        /* reset dx */
        setdeltax(c, deltax(s));

        return(c);
}