## ----setup, include = FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE ) old_par <- par("mfrow", "mai", "pty") options(width = 999) ## ----R_version, include = TRUE, echo = FALSE----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- out <- version cat( "R version: ", out$major, ".", out$minor, "\n", "Generated on: ", format(Sys.time(), "%d-%B-%Y"), sep = "" ) ## ----library------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- library(MatchingPursuit) ## ----empi_locate, include = TRUE, echo = TRUE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- empi_locate() ## ----empi_install, include = TRUE, echo = TRUE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- # empi_install() ## ----empi_check, include = TRUE, echo = TRUE----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- empi_check() ## ----7_non_stationary, include = TRUE, echo = TRUE, warning = FALSE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ fs <- 1024 T <- 1 t <- seq(0, T - 1 / fs, 1 / fs) N <- length(t) # 7 non-stationary signals. x1 <- sin(2 * pi * (10 + 40 * t) * t) # linear chirp x2 <- sin(2 * pi * (20 * t^2) * t) # nonlinear chirp x3 <- (1 + 0.5 * sin(2 * pi * 2 * t)) * sin(2 * pi * 30 * t) # AM x4 <- sin(2 * pi * 50 * t + 5 * sin(2 * pi * 3 * t)) # FM x5 <- exp(-2 * t) * sin(2 * pi * 60 * t) # decreasing amplitude x6 <- sin(2 * pi * (5 + 20 * sin(2 * pi * t)) * t) # frequency modulated sine wave x7 <- t * sin(2 * pi * 40 * t) # increasing amplitude signal <- data.frame(x = x1 + x2 + x3 + x4 + x5 + x6 + x7) ## ----7_non_stationary_plot, include = TRUE, echo = FALSE, warning = FALSE, fig.width = 7, fig.height = 7----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- range <- range(signal) par(mfrow = c(8, 1), pty = "m", mai = c(0.2, 0.4, 0.2, 0.1)) plot(t, signal$x, type = "l", col = "blue", main = "", xlab = "", ylab = "") plot(t, x1, type = "l", ylab = "", xlab = "", xaxt = "n", yaxt = "s", main = "Linear chirp") plot(t, x2, type = "l", ylab = "", xlab = "", xaxt = "n", yaxt = "s", main = "Nonlinear chirp") plot(t, x3, type = "l", ylab = "", xlab = "", xaxt = "n", yaxt = "s", main = "AM") plot(t, x4, type = "l", ylab = "", xlab = "", xaxt = "n", yaxt = "s", main = "FM") plot(t, x5, type = "l", ylab = "", xlab = "", xaxt = "n", yaxt = "s", main = "Decreasing amplitude") plot(t, x6, type = "l", ylab = "", xlab = "", xaxt = "n", yaxt = "s", main = "Frequency modulated sine wave") plot(t, x7, type = "l", ylab = "", xlab = "", xaxt = "n", yaxt = "s", main = "Increasing amplitude") par(old_par) ## ----sample1_csv, include = TRUE, echo = TRUE, warning = FALSE----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- # The sample1.csv file contains exactly the same data as shown in Step 1. file <- system.file("extdata", "sample1.csv", package = "MatchingPursuit") # The first line of the file contains two values: # the sampling rate in Hz (1024 Hz here) and the signal duration # in seconds (1 s here). out <- read.csv(file, header = FALSE) head(out) signal <- read_csv_signals(file) signal ## ----mp_omp_execute_1, include = TRUE, echo = TRUE, warning = FALSE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ fit_mp <- mp_omp_execute( mode = "mp", signal = signal, n_nonzero_coefs = 25 ) summary(fit_mp) ## ----sample1_plot_1, include = TRUE, echo = TRUE, warning = FALSE, fig.width = 7, fig.height = 7------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- plot(fit_mp) ## ----sample1_plot_2, include = TRUE, echo = TRUE, warning = FALSE, fig.width = 7, fig.height = 7------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- out <- tf_map( x = fit_mp, channel = 1, freq_divide = 4, atom_centers = "numbers" ) ## ----empi_execute, include = TRUE, echo = TRUE, warning = FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- # sample1_empi_out <- empi_execute( # signal = signal, # empi_options = "-o local --gabor -i 25", # ) # saveRDS(sample1_empi_out, file = "sample1.rds") ## ----sample1_plot_3, include = TRUE, echo = TRUE, warning = FALSE, fig.width = 7, fig.height = 7------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- file <- system.file("extdata", "sample1.rds", package = "MatchingPursuit") sample1_empi_out <- readRDS(file = file) plot(sample1_empi_out) ## ----read_edf_params, include = TRUE, echo = TRUE, warning = FALSE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- file <- system.file("extdata", "EEG.edf", package = "MatchingPursuit") # Read signal parameters and display them in a tabular form. eeg <- read_edf_signals(file) signal_eeg <- eeg$signal sampling_frequency <- eeg$sampling_frequency time <- seq(0, nrow(signal_eeg) - 1) / sampling_frequency ## ----eeg_filtering, include = TRUE, echo = TRUE, warning = FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- # Filter parameters that will be used (quite typical in filtering EEG signals). fc <- design_filters( sampling_frequency = sampling_frequency, notch = c(49, 51), lowpass = 40, highpass = 1, ) # Filtering input signals. signal_eeg_f <- signal_eeg for (m in 1:ncol(signal_eeg_f)) { signal_eeg_f[, m] = signal::filtfilt(fc$notch, signal_eeg[, m]) # 50Hz notch filter signal_eeg_f[, m] = signal::filtfilt(fc$lowpass, signal_eeg_f[, m]) # Low pass IIR Butterworth signal_eeg_f[, m] = signal::filtfilt(fc$highpass, signal_eeg_f[, m]) # High pass IIR Butterwoth } ## ----read_edf_signal_resampling, include = TRUE, echo = TRUE, warning = FALSE-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- signal_eeg_f_r <- resample_signal(signal = signal_eeg_f, p = 1, q = 2) time_128 <- seq(0, nrow(signal_eeg_f_r) - 1) / (sampling_frequency / 2) sampling_frequency_r <- 128 ## ----filtering_and__resampling_plot, include = TRUE, echo = FALSE, warning = FALSE, fig.width = 7, fig.height = 5-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- old_par <- par("mfrow", "mai", "mgp") par(mfrow = c(3, 1), pty = "m", mai = c(0.35, 0.5, 0.4, 0.1), mgp = c(1.5, 0.5, 0)) rg <- range(c(signal_eeg[, 1], signal_eeg_f[, 1], signal_eeg_f_r[, 1])) plot( time, signal_eeg[, 1], type = "l", xlab = "Time (s)", ylab = "Amplitude", main = "Original EEG signal, channel Fp1", col = "blue", ylim = rg, panel.first = grid() ) abline(h = 0, col = "gray") plot( time, signal_eeg_f[, 1], type = "l", xlab = "Time (s)", ylab = "Amplitude", main = "Filtered EEG signal, channel Fp1", col = "blue", ylim = rg, panel.first = grid() ) abline(h = 0, col = "gray") plot( time_128, signal_eeg_f_r[, 1], type = "l", xlab = "Time (s)", ylab = "Amplitude", main = "Filtered and downsampled EEG signal (128 Hz), channel Fp1", col = "blue", ylim = rg, panel.first = grid() ) abline(h = 0, col = "gray") par(old_par) ## ----eeg_montage, include = TRUE, echo = TRUE, warning = FALSE----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- # Pairs of signals for bipolar montage (so called "double banana"). pairs <- list( c("Fp2", "F4"), c("F4", "C4"), c("C4", "P4"), c("P4", "O2"), c("Fp1", "F3"), c("F3", "C3"), c("C3", "P3"), c("P3", "O1"), c("Fp2", "F8"), c("F8", "T4"), c("T4", "T6"), c("T6", "O2"), c("Fp1", "F7"), c("F7", "T3"), c("T3", "T5"), c("T5", "O1"), c("Fz", "Cz"), c("Cz", "Pz") ) # Make the bipolar montage. signal_eeg_f_r_m <- eeg_montage( signal_eeg_f_r, montage_type = c("bipolar"), bipolar_pairs = pairs ) # Original signal (first 6 rows, first 6 channels). signal_eeg_f_r[1:6, 1:6] # Signal after banana montage (first 6 rows, first 6 channels). signal_eeg_f_r_m[1:6, 1:6] ## ----eeg_empi_execute, include = TRUE, echo = TRUE, warning = FALSE, fig.width = 7, fig.height = 7----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- # The empi_options parameter is NULL, so the EMPI program is # run with the parameters "-o local --gabor -i 50" # # To make the RDS file smaller (CRAN requirement), we select only one channel. # sig <- as_sig(signal_eeg_f_r_m[, 2], sampling_frequency_r) # eeg_empi_out <- empi_execute ( # signal = sig, # empi_options = NULL # ) # # saveRDS(eeg_empi_out, file = "eeg_empi_out.rds") ## ----eeg_tf_map, include = TRUE, echo = TRUE, warning = FALSE, fig.width = 7, fig.height = 7----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- file <- system.file("extdata", "eeg_empi_out.rds", package = "MatchingPursuit") eeg_empi_out <- readRDS(file = file) plot(eeg_empi_out) ## ----plot_eeg_filt_resamp_mont, include = TRUE, echo = TRUE, warning = FALSE, fig.width = 7, fig.height = 7-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- edf_processed <- structure( list( signal = as.data.frame(signal_eeg_f_r_m), sampling_frequency = sampling_frequency_r, time = time_128, signal_names = colnames(signal_eeg_f_r_m), record_name = basename(file) ), class = "edf" ) plot( x = edf_processed, begin = 0, end = 10, panel_height = NULL, rainbow = FALSE, bg_colour = "white", txt_col = "blue", zero_line = TRUE, main = "EEG after filtering, resampling, and double banana montage" ) ## ----plot_orig_eeg, include = TRUE, echo = TRUE, warning = FALSE, fig.width = 7, fig.height = 7-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- plot( x = eeg, begin = 0, end = 10, panel_height = NULL, rainbow = FALSE, bg_colour = "white", txt_col = "blue", zero_line = TRUE, main = "Original EEG before preprocessing" ) ## ----read_wfdb_signals, include = TRUE, echo = TRUE, warning = FALSE----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- file <- system.file("extdata", "00001_lr.hea", package = "MatchingPursuit") out_ecg <- read_wfdb_signals(file) head(out_ecg$signal) out_ecg$sampling_frequency out_ecg$lead_names out_ecg$record_name ## ----ecg_empi_execute, include = TRUE, echo = TRUE, warning = FALSE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ # ecg_empi_out <- empi_execute( # signal = out_ecg # ) # # saveRDS(ecg_empi_out, file = "00001_lr.rds") ## ----00001_lr, include = TRUE, echo = TRUE, warning = FALSE, fig.width = 7, fig.height = 7------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- # Read the previously generated decomposition result. file <- system.file("extdata", "00001_lr.rds", package = "MatchingPursuit") ecg_empi_out <- readRDS(file = file) # Create time-frequency map based on atoms. out <- tf_map( x = ecg_empi_out, channel = 1, verbose = TRUE ) ## ----plot_ecg, include = TRUE, echo = TRUE, warning = FALSE, fig.width = 7, fig.height = 7------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- plot( x = out_ecg, begin = 0, end = 10, panel_height = 1, zero_line = FALSE, small_squares = TRUE ) ## ----chirp_def, include = TRUE, echo = FALSE, warning = FALSE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ n <- 1280 f <- 128 t1 = n / f t <- seq(from = 0, to = n - 1, by = 1) / f g1 <- gabor_atom(n, f, mean = 2, phase = 0, sigma = 1, frequency = 50, normalization = F) g2 <- gabor_atom(n, f, mean = 8, phase = 0, sigma = 1.5, frequency = 25, normalization = F) g3 <- gabor_atom(n, f, mean = 2, phase = 0, sigma = 2, frequency = 30, normalization = F) g4 <- gabor_atom(n, f, mean = 7, phase = 0, sigma = 0.5, frequency = 15, normalization = F) imp <- rep(0, n) imp[n / 2] <- 4 sine <- 0.2 * sin(2 * pi * 3 * t) chirp <- signal::chirp(t = t, f0 = 0, t1 = t1, f1 = 50, form = c("linear"), phase = 0) signal <- g1$gabor + g2$gabor + g3$gabor + g4$gabor + imp + chirp + sine ## ----chirp_plot, include = TRUE, echo = FALSE, warning = FALSE, fig.width = 7, fig.height = 7---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- par(mfcol = c(8, 1), pty = "m", mai = c(0.2, 0.1, 0.2, 0.1)) ylim <- c(-2, 2) plot(t, signal, ylim = ylim, xaxt = "s", yaxt = "n", bty = "o", col = "blue", type = "l", xlab = "", ylab = "", main = "The sum of the below functions") plot(t, g1$gabor, ylim = ylim, xaxt = "n", yaxt = "n", bty = "o", col = "brown", type = "l", xlab = "", ylab = "", main = "Gabor: f = 50Hz, mean = 2, sigma = 1, phase = 0") plot(t, g2$gabor, ylim = ylim, xaxt = "n", yaxt = "n", bty = "o", col = "brown", type = "l", xlab = "", ylab = "", main = "Gabor: f = 25Hz, mean = 8, sigma = 1.5, phase = 0") plot(t, g3$gabor, ylim = ylim, xaxt = "n", yaxt = "n", bty = "o", col = "brown", type = "l", xlab = "", ylab = "", main = "Gabor: f = 30Hz, mean = 2, sigma = 2, phase = 0") plot(t, g4$gabor, ylim = ylim, xaxt = "n", yaxt = "n", bty = "o", col = "brown", type = "l", xlab = "", ylab = "", main = "Gabor: f = 15Hz, mean = 7, sigma = 0.5, phase = 0") plot(t, imp, ylim = ylim, xaxt = "n", yaxt = "n", bty = "o", col = "brown", type = "l", xlab = "", ylab = "", main = "Unit impulse") plot(t, sine, ylim = ylim, xaxt = "n", yaxt = "n", bty = "o", col = "brown", type = "l", xlab = "", ylab = "", main = "Sine wave: 3 Hz") plot(t, chirp, ylim = ylim, xaxt = "n", yaxt = "n", bty = "o", col = "brown", type = "l", xlab = "", ylab = "", main = "Linear-frequency chirp (0Hz - 50Hz)") par(old_par) ## ----chirp_plot_tf_1, include = TRUE, echo = TRUE, warning = FALSE, fig.width = 7, fig.height = 6------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ # sig_file <- system.file("extdata", "sample2.csv", package = "MatchingPursuit") # signal <- read_csv_signals(sig_file, col_names_in_csv = FALSE) # sample2_empi_out <- empi_execute ( # signal = signal # ) # saveRDS(sample2_empi_out, file = "sample2.rds") ## ----chirp_plot_tf_2, include = TRUE, echo = TRUE, warning = FALSE, fig.width = 7, fig.height = 6------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ file <- system.file("extdata", "sample2.rds", package = "MatchingPursuit") sample2_empi_out <- readRDS(file = file) out <- tf_map( x = sample2_empi_out, channel = 1, freq_divide = 1 ) ## ----gabor_atom, include = TRUE, echo = FALSE, fig.width = 7, fig.height = 5, fig.align = 'center'----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- N <- 512 fs <- 256 # normalization = T --> signal's norm = 1 par(mfrow = c(2,2), pty = "m", mai = c(0.4, 0.4, 0.3, 0.2)) n <- 4 sigmas <- c(0.5, 0.2, 0.8, 0.5) frequencies <- c(14, 8, 4, 1) phases <- c(0, 1, 1.5, -2) means = c(0.5, 0.8, 1, 1.5) s <- rep(0, N) for (i in c(1:n)) { sigma <- sigmas[i] frequency <- frequencies[i] phase <- phases[i] mean <- means[i] main <- latex2exp::TeX(paste( "$\\mu=$", means[i], ", ", "$\\sigma=$", sigmas[i], ", ", "$\\f=$", frequencies[i], ", ", "$\\phi=$", phases[i], sep = "" ) ) gb <- gabor_atom(N, fs, mean, phase, sigma, frequency, normalization = FALSE) plot( gb$time, gb$gauss, type="l", ylim = c(-1, 1), col = "red", xlab= "", ylab= "", xaxt = "t", yaxt = "t", bty = "or", cex.axis = 1, lwd = 2, main = main) lines(gb$time, gb$cosine, type = "l", col = "grey") lines(gb$time, gb$gabor, col="blue", lwd = 2) s <- s + gb$gabor } par(old_par) ## ----restore_par, include = FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- par(old_par) ## ----cores_example_1, include = TRUE, echo = TRUE, warning = FALSE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- dictionary <- matrix( c( 0.098, -0.308, -0.342, -0.894, 0.928, 0.674, -0.270, 0.283, 0.145, -0.326, 0.810, -0.114, -0.170, -0.466, 0.377, -0.036, 0.281, -0.357, 0.105, -0.327 ), nrow = 5, byrow = TRUE ) # Although mp_core() and omp_core() normalize dictionary atoms internally, # the atoms are normalized explicitly here so that the coefficients used to # construct the signal have a direct interpretation. dictionary <- sweep(dictionary, 2, sqrt(colSums(dictionary^2)), "/") # Signal constructed from two correlated dictionary atoms signal <- 1.5 * dictionary[, 1] - 0.8 * dictionary[, 2] ## ----cores_example_2, include = TRUE, echo = TRUE, warning = FALSE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- fit_mp <- mp_core( dictionary = dictionary, signal = signal, n_nonzero_coefs = 2 ) fit_omp <- omp_core( dictionary = dictionary, signal = signal, n_nonzero_coefs = 2 ) ## ----cores_example_3, include = TRUE, echo = TRUE, warning = FALSE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- fit_mp$support fit_omp$support fit_mp$coefs fit_omp$coefs sum(fit_mp$residual^2) sum(fit_omp$residual^2) ## ----omp_read_signal, include = TRUE, echo = TRUE, warning = FALSE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- sig_file <- system.file("extdata", "sample1.csv", package = "MatchingPursuit") signal <- read_csv_signals(sig_file, col_names_in_csv = FALSE) ## ----mp_omp_execute_2, include = TRUE, echo = TRUE, warning = FALSE, fig.width = 7, fig.height = 7----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- fit <- mp_omp_execute( mode = 'omp', # or "mp" for classical Matching Pursuit signal = signal, n_nonzero_coefs = 50, topk = 10000, verbose = FALSE ) plot(fit) ## ----generate_xml_dict, include = TRUE, echo = TRUE, warning = FALSE----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- sampling_frequency <- signal$sampling_frequency signal_length <- nrow(signal$signal) duration <- nrow(signal$signal) / sampling_frequency xml_file <- tempfile(fileext = ".xml") dict <- generate_xml_dict( N = signal_length, file = xml_file ) ## ----read_gabor_dict, include = TRUE, echo = TRUE, warning = FALSE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- atoms_dict <- read_gabor_dict( xml_file = xml_file, sampling_frequency = sampling_frequency, duration = duration, full_atoms_in_signal = FALSE, verbose = FALSE ) dim(atoms_dict) head(atoms_dict) ## ----topk_gabor_atoms, include = TRUE, echo = TRUE, warning = FALSE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ dict_topk <- topk_gabor_atoms( atoms_dict = atoms_dict, signal = signal, topk = 10000, verbose = TRUE ) ## ----validation_1, include = TRUE, echo = TRUE, warning = FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- set.seed(1) N <- 256 t <- (0:(N - 1)) / N dictionary <- cbind( sin(2 * pi * 3 * t), cos(2 * pi * 3 * t), sin(2 * pi * 7 * t), cos(2 * pi * 7 * t), sin(2 * pi * 12 * t), cos(2 * pi * 12 * t), sin(2 * pi * 20 * t), cos(2 * pi * 20 * t) ) # Unit-L2 normalization of dictionary atoms. dictionary <- sweep(dictionary, 2, sqrt(colSums(dictionary^2)), "/") colnames(dictionary) <- c( "sin_3", "cos_3", "sin_7", "cos_7", "sin_12", "cos_12", "sin_20", "cos_20" ) ## ----validation_2, include = TRUE, echo = TRUE, warning = FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- true_coef <- c(0, 1.5, 0, 0, -0.8, 0, 0, 0.5) true_atoms <- c("cos_3", "sin_12", "cos_20") signal <- as.vector(dictionary %*% true_coef) ## ----validation_3, include = TRUE, echo = TRUE, warning = FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- out <- omp_core( dictionary = dictionary, signal = signal, n_nonzero_coefs = 3 ) colnames(out$selected_atoms) out$coefs ## ----validation_4, include = TRUE, echo = TRUE, warning = FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- setequal(colnames(out$selected_atoms), true_atoms) relative_error <- sqrt(sum((signal - out$reconstruction)^2)) / sqrt(sum(signal^2)) relative_error ## ----validation_5, include = TRUE, echo = TRUE, warning = FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- signal_noisy <- signal + rnorm(N, sd = 0.1) out_noisy <- omp_core( dictionary = dictionary, signal = signal_noisy, n_nonzero_coefs = 3 ) colnames(out_noisy$selected_atoms) out_noisy$coefs setequal(colnames(out_noisy$selected_atoms), true_atoms) ## ----validation_6, include = TRUE, echo = TRUE, warning = FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- relative_residual_error <- sqrt( sum((signal_noisy - out_noisy$reconstruction)^2) ) / sqrt(sum(signal_noisy^2)) relative_residual_error relative_clean_error <- sqrt( sum((signal - out_noisy$reconstruction)^2) ) / sqrt(sum(signal^2)) relative_clean_error ## ----fig_1, echo=FALSE, out.width="100%", fig.align="center"------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- knitr::include_graphics("../man/figures/main_workflows.png") ## ----generate_dict, include = TRUE, echo = TRUE, warning = FALSE--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- # Generate a dictionary for a 512-sample signal xml_file <- tempfile(fileext = ".xml") dict <- generate_xml_dict ( N = 512, file = xml_file, max_window_length = "3N" ) dict ## ----read_gabor_dict_2, include = TRUE, echo = TRUE, warning = FALSE----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- # Read the dictionary specification and retain only atoms # fully contained within the analyzed signal dict_full <- read_gabor_dict( xml_file = xml_file, sampling_frequency = 256, duration = 2, full_atoms_in_signal = TRUE, verbose = FALSE ) dict_padded <- read_gabor_dict( xml_file = xml_file, sampling_frequency = 256, duration = 2, full_atoms_in_signal = FALSE, verbose = FALSE ) dim(dict_full) dim(dict_padded)