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

#if @HELP_RCEPS

    RCEPS

    Purpose: Calculates the real cepstrum

    Syntax:  RCEPS(s, n)

             s   - input series or array

             n   - an optional integer, the number of samples to use.
                   If n > length(s), the series is zero padded. Defaults
                   to length(s).


    Returns: A real series or array


    Example:
             W1: gtri(100, 1, 1/100)^3
             W2: w1-delay(w1, 60)/2
             W3: rceps(w1)
             W4: rceps(w1, 512)

             W2 adds a synthesized echo at 60 seconds to the data of
             W1.

             W3 displays a small peak at t == 60 indicating the
             presence of the echo. W4 performs the same calculation
             with the data padded to 512 samples.

    Remarks:

            The complex cepstrum of a series is essentially
            IDFT(log(DFT(s))).  However, the complex log calculation
            requires the evaluation of of the continuous phase
            component.  RCEPS ignores the phase component and calculates:

            Real(ifft(log(mag(fft(s)))))

    See Also:
             Cceps
             Iceps


    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


/* real cepstrum - no phase handling */
rceps(s, n)
{
        local f, r;

        if (argc < 2)
        {
                if (argc < 1) error("rceps - input series required");
                
                n = length(s);
        }
        
        if (n < 1) error("rceps - length must be >= 1");

        if (n != length(s))
        {
                f = fft(extract(s, 1, n));
        }
        else
        {
                f = fft(s);
        }

        /* compute IDFT(log(DFT(s)) - no phase handling */
        r = real(ifft(log(mag(f))));

        return(r);
}