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);
}