R/spect.R

In limads/eegr: Evoked potential and spectral analysis tools for EEG data

#' Short-time fast Fourier transform. Returns the normalized spectrum
short_fft <- function(epoch) { 
	fft_epoch <- fft(epoch) / length(epoch)
	fft_epoch[1:as.integer(floor(length(fft_epoch)/2))]
}

#' For a given sigal, return the possible frequency values in Hz
get_freqs <- function(signal, srate) {
	time <- length(signal) / srate
	fres <- 1 / time
	seq(0, fres*(length(signal)/2) - fres, fres)
}

#' Returns the spectrum, vector-averaged across epochs that match cond_rgx
#' Calculate spectral power at specified frequencies or frequency range
#' @param eeg EEG session
#' @param cond_rgx Regex that picks epochs over which function will be applied
#' @return EEG frequency spectrum (complex values) (tibble)
spect <- function(eeg, cond_rgx) { 

	use_epochs <- extract_epochs(eeg, cond_rgx)
	sz <- epoch_size(use_epochs)
	use_fft <- eeg@spec_epochs[str_detect(names(eeg@spec_epochs), cond_rgx)]
	
	spec_mat <- matrix(0+0i, nrow=length(use_fft), ncol=sz)

	for(i in 1:length(use_epochs)) {
		spec_mat[i,] <- use_fft[[i]]
	}

	spec_vals <- apply(spec_mat, 2, mean)[1:as.integer(floor(sz/2))]
	freqs <- get_freqs(use_epochs[[1]], srate=eeg@srate)
	
	return(tibble(
		freq  = freqs,
		spect = spec_vals,
	))
}

#' Given a spectrum returned by spect, returns the phase
#' @param spect Spectrum returned by spect
#' @param deg Returns phase in degrees (TRUE) or radians (FALSE)
#' @return EEG phase
spect_phase <- function(spect, deg=FALSE) { 
	phase <- spect %>% pull(spect) %>% Arg
	if(deg) {
		return(180 / pi) * phase
	} else {
		phase
	}
}

#' Given a spectrum returned by spec, returns the power.
#' @param spect Spectrum returned by spect
#' @param power Returns spectrum amplitude (FALSE) or amplitude squared (TRUE)
#' @param norm Returns normalized power (TRUE) or not (FALSE)
#' @return EEG power
spect_power <- function(spect, power=TRUE, norm=FALSE) {
	ampl <- spect %>% pull(spect) %>% Mod
	if(norm) {
		return(ampl / max(ampl))
	} else {
		if(power) {
			return(ampl^2)
		} else{
			return(ampl)
		}
	}
}

#' Given a vector of complex Fourier coefficients, calculate the T2-Circ statistic
#' (from Victor & Mast, 1991)
calculate_tcirc <- function(z){
  n  <-  length(z)
  zest <-  mean(z)
  vind  <-  sum(abs(z - zest)^2)
  tcirc  <-  (n-1) * abs(zest)^2 / vind
  return(tcirc)
}

#' Given a vector of complex Fourier coefficients, calculate the p-value of
#' the T2-Circ statistic
calculate_pval <- function(z){
  tcirc <- calculateTCirc(z)
  n     <- length(z)
  pval  <- 1 - pf(tcirc, 2, 2*n-2)
  return(pval)
}

#' Given a vector of Fourier coefficients, calculate the confidence region
#' of the T2-circ statistic.
calculate_conf <- function(z){
  tcirc <- calculateTCirc(z)
}