rceps: Real cepstrum

In gjmvanboxtel/gsignal: Signal Processing

Real cepstrum

Description

Return the real cepstrum and minimum-phase reconstruction of a signal

Usage

rceps(x, minphase = FALSE)

Arguments

x	input data, specified as a real vector.
minphase	logical (default: FALSE) indication whether to compute minimum-phase reconstructed signal

Details

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.

Value

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

Author(s)

Paul Kienzle, [e-mail],
Mike Miller.
Conversion to R by Geert van Boxtel, [e-mail].

References

https://en.wikipedia.org/wiki/Minimum_phase

See Also

cceps

Examples

## 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.

gjmvanboxtel/gsignal documentation built on Nov. 22, 2023, 8:19 p.m.