Real cepstrum
Return the real cepstrum and minimum-phase reconstruction of a signal
rceps(x, minphase = FALSE)
x |
input data, specified as a real vector. |
minphase |
logical (default: |
Cepstral analysis is a nonlinear signal processing technique that is applied most commonly in speech and image processing, or as a tool to investigate periodic structures within frequency spectra, for instance resulting from echos/reflections in the signal or to the occurrence of harmonic frequencies (partials, overtones).
The cepstrum is used in many variants. Most important are the power cepstrum,
the complex cepstrum, and real cepstrum. The function rceps implements
the real cepstrum by computing the inverse of the log-transformed FFT while
discarding phase, i.e.,
rceps(x) <- ifft(log(Mag(fft(x))))
The real cepstrum is related to the power spectrum by the relation pceps = 4 * rceps^2.
The function rceps() can also return a minimum-phase reconstruction of
the original signal. The concept of minimum phase originates from filtering
theory, and denotes a filter transfer function with all of its poles and
zeroes in the Z-transform domain lie inside the unit circle on the complex
plane. Such a transfer function represents a stable filter.
A minimum-phase signal is a signal that has its energy concentrated near the front of the signal (near time 0). Such signals have many applications, e.g. in seismology and speech analysis.
If minphase equals FALSE, the real cepstrum is returned
as a vector. If minphase equals TRUE, a list is returned
containing two vectors; y containing the real cepstrum, and
ym containing the minimum-phase reconstructed signal
Paul Kienzle, pkienzle@users.sf.net,
Mike Miller.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.
## Simulate a speech signal with a 70 Hz glottal wave
f0 <- 70; fs = 10000 # 100 Hz fundamental, 10 kHz sampling rate
a <- Re(poly(0.985 * exp(1i * pi * c(0.1, -0.1, 0.3, -0.3))))
s <- 0.05 * runif(1024)
s[floor(seq(1, length(s), fs / f0))] <- 1
x <- filter(1, a, s)
## compute real cepstrum and min-phase of x
cep <- rceps(x, TRUE)
hx <- freqz(x, fs = fs)
hxm <- freqz (cep$ym, fs = fs)
len <- 1000 * trunc(min(length(x), length(cep$ym)) / 1000)
time <- 0:(len-1) * 1000 / fs
op <- par(mfcol = c(2, 2))
plot(time, x[1:len], type = "l", ylim = c(-10, 10),
xlab = "Time (ms)", ylab = "Amplitude",
main = "Original and reconstructed signals")
lines(time, cep$ym[1:len], col = "red")
legend("topright", legend = c("original", "reconstructed"),
lty = 1, col = c(1, 2))
plot(time, cep$y[1:len], type = "l",
xlab = "Quefrency (ms)", ylab = "Amplitude",
main = "Real cepstrum")
plot (hx$w, log(abs(hx$h)), type = "l",
xlab = "Frequency (Hz)", ylab = "Magnitude",
main = "Magnitudes are identical")
lines(hxm$w, log(abs(hxm$h)), col = "red")
legend("topright", legend = c("original", "reconstructed"),
lty = 1, col = c(1, 2))
phx <- unwrap(Arg(hx$h))
phym <- unwrap(Arg(hxm$h))
range <- c(round(min(phx, phym)), round(max(phx, phym)))
plot (hx$w, phx, type = "l", ylim = range,
xlab = "Frequency (Hz)", ylab = "Phase",
main = "Unwrapped phase")
lines(hxm$w, phym, col = "red")
legend("bottomright", legend = c("original", "reconstructed"),
lty = 1, col = c(1, 2))
par(op)
## confirm the magnitude spectrum is identical in the signal
## and the reconstruction and that there are peaks in the
## cepstrum at 14 ms intervals corresponding to an F0 of 70 Hz.Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.