From 05f96535ba5a66b422f2c2dad7b94f1c98ba4313 Mon Sep 17 00:00:00 2001
From: mcrucifix <michel.crucifix@uclouvain.be>
Date: Tue, 8 Oct 2024 13:13:14 +0200
Subject: [PATCH] wrap up + doc
---
R/data_description.R | 34 +++++++++++++++++
R/mfft_support.R | 72 ++++++++++++++++++++++++++++++++++++
man/harmonic_sample.Rd | 15 ++++++++
man/harmonic_sample_noisy.Rd | 15 ++++++++
man/mfft_complex.Rd | 42 +++++++++++++++++++++
man/pseudo_log.Rd | 15 ++++++++
6 files changed, 193 insertions(+)
create mode 100644 R/data_description.R
create mode 100644 R/mfft_support.R
create mode 100644 man/harmonic_sample.Rd
create mode 100644 man/harmonic_sample_noisy.Rd
create mode 100644 man/mfft_complex.Rd
create mode 100644 man/pseudo_log.Rd
diff --git a/R/data_description.R b/R/data_description.R
new file mode 100644
index 0000000..787a357
--- /dev/null
+++ b/R/data_description.R
@@ -0,0 +1,34 @@
+#' harmonic_sample
+#'
+#' A simple time series for testing purposes
+#' cos(t*0.13423167+0.00) + 1.3 * cos(t*0.119432+2.314) + 0.134994 + 0.4*cos(t*0.653167) + 0.11 * cos(t*0.78913498)
+#' for t from 0 to 1023. It is coded as a `mfft_deco` objcet and self-describes its spectral decomposition.
+#'
+#' @docType data
+#' @keywords datasets
+#' @name harmonic_sample
+#' @usage data(harmonic_sample)
+NULL
+
+#' harmonic_sample_noisy
+#'
+#' Same has harmonic sample but with an added random time series
+#' cos(t*0.13423167+0.00) + 1.3 * cos(t*0.119432+2.314) + 0.134994 + 0.4*cos(t*0.653167) + 0.11 * cos(t*0.78913498) + rnorm(1024)*0.12
+#' for t from 0 to 1023. It is coded as a `mfft_deco` objcet and self-describes its spectral decomposition.
+#'
+#' @docType data
+#' @keywords datasets
+#' @name harmonic_sample_noisy
+#' @usage data(harmonic_sample_noisy)
+NULL
+
+#' pseudo_log
+#'
+#' A truncated version (squared-shape) version of a harmonic signal for mimicking a sedimentary log
+#' sign( cos(t*0.13423167+0.00) + 1.3 * cos(t*0.119432+2.314) + 0.134994 + 0.4*cos(t*0.653167) + 0.11 * cos(t*0.78913498) + 0.8))
+#' for t from 0 to 8191. It is coded as a `mfft_deco` objcet and self-describes its spectral decomposition.
+#' @docType data
+#' @keywords datasets
+#' @name pseudo_log
+#' @usage data(pseudo_log)
+NULL
diff --git a/R/mfft_support.R b/R/mfft_support.R
new file mode 100644
index 0000000..22df4ae
--- /dev/null
+++ b/R/mfft_support.R
@@ -0,0 +1,72 @@
+#' MFFT reconstruction
+#' @rdname mfft_deco
+#' @export reconstruct_mfft
+reconstruct_mfft <- function(M){
+ if (!(attr(M,"class")=="mfft_deco")) stop ("object is not a MFFT decomposition")
+ xdata <- attr(M,"xdata")
+ nfreq <- attr(M,"nfreq")
+ times <- seq(length(xdata))*dt(xdata) + startx(xdata)
+ reconstructed <- lapply(seq(nfreq), function(i) M$Amp[i]*cos(M$Freq[i]*times + M$Phase[i]) )
+}
+
+#' MFFT ANOVA
+#' not ready. do not use.
+#' @rdname mfft_deco
+#' @export mfft_anova
+mfft_anova <- function(M){
+ if (!(attr(M,"class")=="mfft_deco")) stop ("object is not a MFFT decomposition")
+ xdata <- attr(M,"xdata")
+ nfreq <- attr(M,"nfreq")
+ N <- length(xdata)
+ times <- seq(length(xdata))*deltat(xdata) + start(xdata)
+ reconstructed <- sapply(seq(nfreq), function(i) M$Amp[i]*cos(M$Freq[i]*times + M$Phase[i]) )
+ cum_reconstruct <- apply(reconstructed, 1, cumsum)
+ residual_vars <- apply(apply(cum_reconstruct, 1, function(i) xdata-i) , 2, function(j) sum(j^2))
+ var0 <- sum(xdata^2)
+ master_vars <- c(var0, residual_vars[-length(residual_vars)])
+ p2s <- 2*seq(nfreq)
+ p1s <- c(0, p2s[-length(p2s)])
+ F <- (master_vars - residual_vars)/residual_vars * (N - p2s)/(p2s-p1s)
+}
+
+
+#' @rdname mfft_deco
+#' @export
+as.data.frame.mfft_deco <- function(x) {data.frame(Freq=x$Freq, Amp=x$Amp, Phases=x$Phases)}
+
+
+#' @rdname mfft_deco
+#' @export
+plot.mfft_deco <- function (M,periods=FALSE,...){
+# O <- order(M$Freq)
+ plot(abs(M$Freq), abs(M$Amp),'h',ylab="Amplitudes", xlab="", ...)
+ if (periods) {
+ frequencies <- pretty(range(M$Freq/(2*pi)))
+ labels <- as.character(1/frequencies)
+ if (0 %in% frequencies) labels[which(frequencies == 0)] = "∞"
+ axis(1, line=3, at=2*pi*frequencies, labels=labels)
+ mtext("Rate", 1, 2)
+ mtext("Period", 1, 4)
+ } else {
+ mtext("Rate", 1, 3)
+ }
+ points(abs(M$Freq), abs(M$Amp),'p',...)
+}
+
+#' @rdname mfft_deco
+#' @export
+lines.mfft_deco <- function (M,...){
+# O <- order(M$Freq)
+ lines(abs(M$Freq), abs(M$Amp),'h',...)
+ points(abs(M$Freq), abs(M$Amp),'p',...)
+}
+
+
+#' @rdname mfft_deco
+#' @export
+print.mfft_deco <- function (M,...){
+ print.data.frame(cbind(as.data.frame(M), Period=2*pi/M$Freq))
+}
+
+
+
diff --git a/man/harmonic_sample.Rd b/man/harmonic_sample.Rd
new file mode 100644
index 0000000..beaa588
--- /dev/null
+++ b/man/harmonic_sample.Rd
@@ -0,0 +1,15 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/data_description.R
+\docType{data}
+\name{harmonic_sample}
+\alias{harmonic_sample}
+\title{harmonic_sample}
+\usage{
+data(harmonic_sample)
+}
+\description{
+A simple time series for testing purposes
+cos(t*0.13423167+0.00) + 1.3 * cos(t*0.119432+2.314) + 0.134994 + 0.4*cos(t*0.653167) + 0.11 * cos(t*0.78913498)
+for t from 0 to 1023. It is coded as a `mfft_deco` objcet and self-describes its spectral decomposition.
+}
+\keyword{datasets}
diff --git a/man/harmonic_sample_noisy.Rd b/man/harmonic_sample_noisy.Rd
new file mode 100644
index 0000000..cc57b10
--- /dev/null
+++ b/man/harmonic_sample_noisy.Rd
@@ -0,0 +1,15 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/data_description.R
+\docType{data}
+\name{harmonic_sample_noisy}
+\alias{harmonic_sample_noisy}
+\title{harmonic_sample_noisy}
+\usage{
+data(harmonic_sample_noisy)
+}
+\description{
+Same has harmonic sample but with an added random time series
+cos(t*0.13423167+0.00) + 1.3 * cos(t*0.119432+2.314) + 0.134994 + 0.4*cos(t*0.653167) + 0.11 * cos(t*0.78913498) + rnorm(1024)*0.12
+for t from 0 to 1023. It is coded as a `mfft_deco` objcet and self-describes its spectral decomposition.
+}
+\keyword{datasets}
diff --git a/man/mfft_complex.Rd b/man/mfft_complex.Rd
new file mode 100644
index 0000000..0afdc3b
--- /dev/null
+++ b/man/mfft_complex.Rd
@@ -0,0 +1,42 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/mfft_complex.R
+\name{mfft_complex}
+\alias{mfft_complex}
+\title{Frequency Modified Fourier transform for Complex Numbers}
+\usage{
+mfft_complex(xdata, nfreq = 30, minfreq = NULL, maxfreq = NULL, correction = 1)
+}
+\arguments{
+\item{xdata}{The data provided either as a time series (advised), or as a vector.
+may be complex}
+
+\item{nfreq}{is the number of frequencies returned, must be smaller that the length of xdata.}
+
+\item{minfreq, maxfreq}{If provided, bracket the frequencies to be probed. Note this are
+more precisely angular velocities (2\pi / period), expressed in time-inverse units
+with the time resolution encoded in `xdata` if the latter is a time series.}
+
+\item{correction}{:
+ Modified Fourier Transform if correction = 0;
+ Frequency Modified Fourier Transform if correction = 1;
+ FMFT with additional non-linear correction if correction = 2
+(while the first algorithm is app. 3 times faster than the third one,
+the third algorithm should be in general much more precise).
+The computed frequencies are in the range given by minfreq and maxfreq.}
+}
+\value{
+a `mfft_deco` object, based on a data.frame with columns "Freq", "Ampl" and "Phases".
+ note that because of a language glitch (to be fixed), "Freq" actually means "Rate"
+}
+\description{
+Implementation of the the Frequency Modified Fourier Transform
+(Sidlichovsky and Nesvorny 1997, Cel. Mech. 65, 137).
+Given a quasi--periodic complex signal X + iY, the algorithm
+estimates the frequencies (f_j), amplitudes (A_j) and phases
+(psi_j) in its decomposition.
+}
+\author{
+Michel Crucifix for the R code, and David Nesvorny for most of the supporting C code doing the
+actual computations
+\insertRef{sidlichovsky97aa}{gtseries}
+}
diff --git a/man/pseudo_log.Rd b/man/pseudo_log.Rd
new file mode 100644
index 0000000..590aef4
--- /dev/null
+++ b/man/pseudo_log.Rd
@@ -0,0 +1,15 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/data_description.R
+\docType{data}
+\name{pseudo_log}
+\alias{pseudo_log}
+\title{pseudo_log}
+\usage{
+data(pseudo_log)
+}
+\description{
+A truncated version (squared-shape) version of a harmonic signal for mimicking a sedimentary log
+sign( cos(t*0.13423167+0.00) + 1.3 * cos(t*0.119432+2.314) + 0.134994 + 0.4*cos(t*0.653167) + 0.11 * cos(t*0.78913498) + 0.8))
+for t from 0 to 8191. It is coded as a `mfft_deco` objcet and self-describes its spectral decomposition.
+}
+\keyword{datasets}
--
GitLab