Language Reference

RCEPSTRUM (x) ;

This function is supported by the IML procedure and the iml action.

The RCEPSTRUM function computes the real cepstrum of a time series. The real cepstrum of the time series x is defined as

where is the magnitude of the Fourier transform of x and .

The input argument to the RCEPSTRUM function is as follows:

The RCEPSTRUM function returns a column vector that contains the real cepstrum of the input data. The length of the output is the same as the length of the input data.

The following example shows how to use the RCEPSTRUM function to estimate the period of a periodic signal. First an input signal is generated that consists of four segments. Each segment is 1 second long; the first 0.75 seconds contain the sum of two sinusoidal signals, and the next 0.25 seconds contain all zeros. This segment is repeated four times, and Gaussian random noise is added to form a pseudoperiodic signal. The output complex spectrum in Figure 361 shows a spike around second, which is a good estimation of the period.

sf  = 100;    /* sampling frequency */
T = 1.0/sf;   /* sampling period */
N = 400;
t = do(0, (N-1)*T, T);
t1 = t[, 1:N*3/16];

f1 = 15;
f2 = 25;
pi = constant("pi");
x1 = 5*sin(2*pi*f1*t1);
x2 = sin(2*pi*f2*t1);
x3 = x1 + x2;        /* sum of two sinusoidal signals */
x4 = j(1, N/16, 0);
x5 = repeat(x3 || x4, 1, 4);

call randseed(1);
noise_mean = 0;
noise_std  = 0.15;
xn = randfun(N, "Normal", noise_mean, noise_std);
x = x5` + xn;  /* input signal */

y = rcepstrum(x);

title "Input Signal";
call series(t,x) grid= {X Y}
                 label={"Time (seconds)" "Amplitude"};

title "Output Real Cepstrum";
call series(t,y) grid= {X Y}
                 label={"Time (seconds)" "Amplitude"};

Figure 361: Real Cepstrum Example