Overview

Beyond its causal-inference and survey-statistics surfaces, rmorie includes signal-processing functions: periodogram, multitaper and Welch spectral estimates, the short-time Fourier transform, the Hilbert envelope, cepstra and cepstral pitch estimation. They apply to audio and biomedical signals and to any regularly sampled series, such as daily incident counts.

A synthetic two-tone signal

library(rmorie)

set.seed(42)
fs <- 1000                           # sampling rate (Hz)
t  <- seq(0, 1, by = 1 / fs)
sig <- sin(2 * pi * 50 * t) +
       0.5 * sin(2 * pi * 120 * t) +
       0.3 * rnorm(length(t))

Basic FFT

spec <- stats::fft(sig)
n    <- length(sig)
freq <- (0:(n / 2 - 1)) * fs / n
mag  <- Mod(spec)[1:(n / 2)]
peak_freqs <- freq[order(mag, decreasing = TRUE)[1:5]]
peak_freqs
#> [1]  49.95005 119.88012 120.87912  50.94905  48.95105

The largest peaks sit next to 50 Hz and 120 Hz, the two components of the synthetic signal.

Spectral density estimates

morie_dsp_psd_periodogram() returns the raw periodogram, morie_psdmt() a multitaper estimate (nw = 4 by default), and welch() Welch’s averaged-segment estimate. Each returns the frequencies and the power spectral density:

per <- morie_dsp_psd_periodogram(sig, fs = fs)
mt  <- morie_psdmt(sig, fs = fs)
wl  <- welch(sig, fs = fs, nperseg = 256L)

# frequency of the largest PSD value in each band (Hz)
peaks <- function(f, p) c(
  low  = f[f > 25 & f < 85][which.max(p[f > 25 & f < 85])],
  high = f[f > 95 & f < 145][which.max(p[f > 95 & f < 145])]
)
rbind(periodogram = peaks(per$freqs, per$psd),
      multitaper  = peaks(mt$frequencies, mt$psd),
      welch       = peaks(wl$extra$freqs, wl$filtered))
#>                  low     high
#> periodogram 49.95005 119.8801
#> multitaper  49.95005 120.8791
#> welch       50.78125 121.0938
mt$n_tapers
#> [1] 7

All three put a peak next to 50 Hz and next to 120 Hz. Welch’s estimate trades frequency resolution (here fs / nperseg, about 3.9 Hz) for lower variance, so its peaks fall on the nearest bins of that coarser grid. The Welch PSD values are in the filtered element and its frequencies in extra$freqs.

Short-time Fourier transform

morie_rangayyan_stft() returns the spectrogram matrix Sxx (frequencies by time segments):

st <- morie_rangayyan_stft(sig, fs = fs, nperseg = 128L)
dim(st$Sxx)
#> [1] 65 14
range(st$times)
#> [1] 0.064 0.896

Envelope and cepstral pitch

morie_dsp_hilbert_envelope() returns the amplitude envelope of a signal. morie_cepstrum_pitch() estimates a fundamental frequency from the real cepstrum; on a harmonic signal at 150 Hz:

env <- morie_dsp_hilbert_envelope(sig)
length(env)
#> [1] 1001

fs2 <- 8000
tt  <- seq(0, 0.5, by = 1 / fs2)
voiced <- rowSums(sapply(1:10, function(k) sin(2 * pi * k * 150 * tt) / k))
morie_cepstrum_pitch(voiced, fs = fs2)$f0
#> [1] 150.9434

Where to go next

• morie_dsp_complex_cepstrum() and morie_dsp_homomorphic() cover cepstral and homomorphic processing.
• Spectrogram(), IStft() and reassigned_spectrogram() are further time-frequency tools; see their help pages.