diff --git a/DESCRIPTION b/DESCRIPTION
index e08bf880023288260cf8d51f04c865d7a9c05dfe..882f36c7149f7e16c9d506285738de9d0ba7bdd7 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -8,3 +8,6 @@ Maintainer: Michel Crucifix <michel.crucifix@uclouvain.be>
Description: More about what it does (maybe more than one line)
License: What license is it under?
LazyLoad: yes
+RoxygenNote: 7.2.3
+Imports: Rdpack
+RdMacros: Rdpack
diff --git a/NAMESPACE b/NAMESPACE
index e300e1602c4fb526f95cbb1e2c99a9885486f4be..3d31d16518e06486b61eb08afa0a6bdc0594f725 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,10 +1,13 @@
-export(cwt_morlet)
-export(mfft)
+# Generated by roxygen2: do not edit by hand
+
+S3method(plot,SSAObject)
+S3method(plot,wavelet)
export(approx_ts)
+export(arspec)
+export(cwt_morlet)
export(mem)
-export(phase_randomize)
+export(mfft)
+export(powerspectrum.wavelet)
export(ssa)
+importFrom(Rdpack,reprompt)
useDynLib(gtseries)
-S3method(plot,wavelet)
-S3method(plot,memObject)
-importFrom(grDevices, colorRampPalette)
diff --git a/NAMESPACE.old b/NAMESPACE.old
new file mode 100644
index 0000000000000000000000000000000000000000..e300e1602c4fb526f95cbb1e2c99a9885486f4be
--- /dev/null
+++ b/NAMESPACE.old
@@ -0,0 +1,10 @@
+export(cwt_morlet)
+export(mfft)
+export(approx_ts)
+export(mem)
+export(phase_randomize)
+export(ssa)
+useDynLib(gtseries)
+S3method(plot,wavelet)
+S3method(plot,memObject)
+importFrom(grDevices, colorRampPalette)
diff --git a/R/approx_ts.R b/R/approx_ts.R
index b28816a1258ae1aae169bac77289731c6f5f1b83..942359b7c0add232425c579d0fcf6fb3cfc797bb 100644
--- a/R/approx_ts.R
+++ b/R/approx_ts.R
@@ -1,4 +1,24 @@
- approx_ts <- function (data,tcoord="x",dcoord="y",scale=1,n=2048,thin=FALSE,spline=FALSE,xlim=NULL)
+#' Produce an evenly-spaced time series by simple interpolation
+#'
+#' Uses standard R routine 'approx' or 'spline' to generate an evenly-spaced time series
+#' based on data supplied on a time coordinate grid.
+#' @param data The data input, typically a dataframe,
+#' with two columns
+#' @param tcoord The name of the column name of the time coordinate. Defaults to "x"
+#' @param dcoord The name of the column name of the data coordinate. Defaults to "y"
+#' @param scale Time scaling (the coordinate "tcoord" will be divided by `scale`.
+#' @param thin Discards data that are sampled with a time space less than 0.5 times the time resolution in the output
+#' @param n (defaults to 2048) is the number of data points in the output
+
+#' @param spline : if True, the standard `spline` routine will be used, otherwise, `approx'
+#' @param xlim : if provided, brackets the time output
+#' @return a time series R object
+#' @author Michel Crucifix
+
+
+#' @importFrom Rdpack reprompt
+#' @export approx_ts
+approx_ts <- function (data,tcoord="x",dcoord="y",scale=1,n=2048,thin=FALSE,spline=FALSE,xlim=NULL)
{
local({
x <- data[tcoord][[1]]
diff --git a/R/arspec.R b/R/arspec.R
index e6807c2c988a6bd80d400f63578acaca4a459037..48b7c57f4ff3484ba201a4d59614cd68d21fee7b 100644
--- a/R/arspec.R
+++ b/R/arspec.R
@@ -1,7 +1,10 @@
-## Calculates the theoretical density frequency of an
-## autoregressive process, according to equation
-## (3.3) of Akaike (1969)
+#' Theoretical Spectrum of autoregressive spectrum
+#'
+#' Calculates the theoretical density frequency of an
+#' autoregressive process, according to equation
+#' (3.3) of Akaike (1969)
+#' @export arspec
arspec <- function (f,ar){
Exp <- exp(-1i*2*pi*outer(f,seq_along(ar),"*"))
1./abs(1-Exp%*%ar)^2
diff --git a/R/cwt_morlet.R b/R/cwt_morlet.R
index d000e0e72b0032c2442463a3e005c23f6219dded..e6c871b0adf76feb090edb3485646d3f9a638b5d 100644
--- a/R/cwt_morlet.R
+++ b/R/cwt_morlet.R
@@ -17,6 +17,7 @@
## Tchawitchia P., Wavelets, functions and operators, ch. 3 in
## Wavelets : theory and applications, Erlenbacher et al. eds, Oxford University Press 1996
+
col_wavelet <- colorRampPalette(c('darkblue','lightblue','grey','white','red'))(100)
search_ridge <- function (A,first_guess,k0=5.6,tol=1e-6, itmax=20,B=seq(along=A),window=c(0,Inf))
@@ -117,6 +118,9 @@ cross_morlet <- function(A, B, ...)
(CA * Conj(CB)) / (Mod(CA) * Mod(CB) )
}
+#' Continous Morlet Wavelet Transform
+
+#' @export cwt_morlet
cwt_morlet <- function (A,inter=20,k0=5.6,amin=1,amax=Inf,calcmask=TRUE,scale=NA,deriv=FALSE)
{
y <- A
@@ -192,6 +196,8 @@ cwt_morlet <- function (A,inter=20,k0=5.6,amin=1,amax=Inf,calcmask=TRUE,scale=NA
wave
}
+#' @rdname cwt_morlet
+#' @export
plot.wavelet <- function (wave,resx=400,resy=300,xlab="Time",ylab="Period",scaling_correction=0,col=col_wavelet,legend=FALSE,Mode=Mod,plotMask=TRUE,...)
{
@@ -216,7 +222,8 @@ plot.wavelet <- function (wave,resx=400,resy=300,xlab="Time",ylab="Period",scali
par(oma=c(2,2,2,5))}
}
-
+#' @rdname cwt_morlet
+#' @export
powerspectrum.wavelet <- function (wave,...)
{
aa <- attr(wave,"scale")
diff --git a/R/mem.R b/R/mem.R
index b921ff82940c46a8b57bff7dcb094dd42bb0b908..36fb6b45d944ee7212fa2c57aba3a1be16475446 100644
--- a/R/mem.R
+++ b/R/mem.R
@@ -1,3 +1,7 @@
+
+#' Maximum entropy method
+
+#' @export mem
mem <- function (A,method='burg',order.max=90,...)
{
deltat<-deltat(A)
diff --git a/R/mfft.R b/R/mfft.R
index 9827cce5e8db8a763c93c2499599dcfeae7aa649..e0a1579c17224ccad94d3fe48c24eebf3742c6c9 100644
--- a/R/mfft.R
+++ b/R/mfft.R
@@ -3,92 +3,94 @@
# check that it provides the same output as the original code
# espacially with the new window
-
-
-mfft <- function(xdata, ydata=NULL, minfreq=NULL, maxfreq=NULL, flag=2, window="hanning", nfreq=30)
+#' Modified Fourier transform
+#'
+#' 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.
+#' @useDynLib gtseries
+#' @param xdata The data provided either as a time series (advised), or as a vector.
+#' may be complex
+#' @param 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.
+#' @param flag :
+#' Modified Fourier Transform if flag = 1;
+#' Frequency Modified Fourier Transform if flag = 2;
+#' FMFT with additional non-linear correction if flag = 3
+#' (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.
+#' @param nfreq is the number of frequencies returned, must be smaller that the length of xdata.
+#' @return STILL NEED TO BE DESCRIBED
+#' @references
+#' @author Michel Crucifix for the R code, and David Nesvorny for most of the supporting C code doing the
+#' actual computations
+#' \insertRef{sidlichovsky97aa}{gtseries}
+#
+#' @export mfft
+mfft <- function(xdata, minfreq=NULL, maxfreq=NULL, flag=1, nfreq=30)
{
- # please at the moent do not play with 'flag'
+ xdata = as.ts(xdata)
+ dt = deltat(xdata)
+ print('deltat')
+ print(dt)
+ startx = start(xdata)
+ ydata = Im(xdata)
+ xdata = Re(xdata)
ndata = length(xdata);
- if (is.null(ydata)) {
- ydata = rep(0, ndata)} else {
- if (length(ydata) != ndata) stop("I am confused: you provide a ydata that does not match the length of xdata")
- }
- if (window=="hanning")
- {
- hanning = (1-cos(2*pi*(seq(ndata)-1)/(ndata-1)))/2.
- # windowing
- xdata = xdata*hanning
- ydata = ydata*hanning
- } else if (window=="none") {} else stop ("this window is not supported, sorry")
- # padding
- N = 2^ceiling(log(ndata)/log(2))
- if (N > ndata){
- xdata[(ndata+1):N] = 0
- ydata[(ndata+1):N] = 0
- }
-
- xdata = ts(xdata)
- ydata = ts(ydata, start=start(xdata), deltat=deltat(xdata))
- dt = deltat(xdata)
+
+ if (!(flag == 1 || flag == 2 || flag == 3)) stop("flag must be either 1, 2, or 3")
+
+# if (window=="hanning")
+# {
+# hanning = (1-cos(2*pi*(seq(ndata)-1)/(ndata-1)))/2.
+# windowing
+# xdata = xdata*hanning
+# ydata = ydata*hanning
+# } else if (window=="none") {} else stop ("this window is not supported, sorry")
+# padding
+# N = 2^ceiling(log(ndata)/log(2))
+# if (N > ndata){
+# xdata[(ndata+1):N] = 0
+# ydata[(ndata+1):N] = 0
+# }
+#
if (is.null(minfreq)){
- my_minfreq = pi/ndata
+ my_minfreq = -pi
my_maxfreq = pi} else {
my_minfreq = minfreq * dt
my_maxfreq = maxfreq * dt}
- # to do: separate minfreq from maxfreq
-
signal1 = as.double(rep(0,3*nfreq))
signal2 = as.double(rep(0,3*nfreq))
signal3 = as.double(rep(0,3*nfreq))
-
- print('length signal1')
- print(length(as.double(signal1)))
- print(length(as.double(signal2)))
- print(length(as.double(signal3)))
- print(nfreq)
-
- print(c("xdata:",xdata[1],xdata[2], xdata[3]))
- print(c("ydata:",ydata[1],ydata[2], ydata[3]))
-
Infreq=as.integer(nfreq)
Iminfreq = as.double(my_minfreq)
Imaxfreq = as.double(my_maxfreq)
Ixdata = as.double(xdata)
Iydata = as.double(ydata)
+ OUT = .C("fmft", Infreq, Iminfreq, Imaxfreq, as.integer(flag),
+ as.integer(ndata), Ixdata, Iydata, signal1,signal2,signal3, DUP=TRUE)
- print('length xdata')
- print(length(Ixdata))
- print(length(Iydata))
- print(length(signal1))
- print(N)
-
- print("N")
- print(N)
- print(Ixdata)
- print("Iydata")
- print(Iydata)
-# print (c("fmft", Infreq, Iminfreq, Imaxfreq, as.integer(flag), as.integer(N), Ixdata, Iydata))
- OUT = .C("fmft", Infreq, Iminfreq, Imaxfreq, as.integer(flag), as.integer(N), Ixdata, Iydata, signal1,signal2,signal3, DUP=TRUE)
-# OUT = .C("fmft", nfreq, minfreq, maxfreq, as.integer(flag), as.integer(ndata), xdata, ydata, signal1,signal2,signal3, DUP=FALSE)
-
- OUTVEC <- t(matrix(OUT[[8]], 3, nfreq))
- print("OUTVEC")
- print(OUTVEC)
+ OUTVEC <- t(matrix(OUT[[7+flag]], 3, nfreq))
+
Freq <- OUTVEC[seq(nfreq) ]/dt
- Ampl <- OUTVEC[nfreq+seq(nfreq) ] * N/ndata
- Phase <- OUTVEC[2*nfreq+seq(nfreq) ]
+ Ampl <- OUTVEC[nfreq+seq(nfreq) ]
+ Phase <- OUTVEC[2*nfreq+seq(nfreq) ] - startx*Freq
# if this is a real vector, there will be positive and negative frequencies corresponding to the same amplitude
# take care and verify that this actually works this way
- return(list(Freq=Freq, Ampl=Ampl, Phase=Phase))
+ return(data.frame(Freq=Freq, Ampl=Ampl, Phase=Phase))
}
#
@@ -193,4 +195,4 @@ mfft <- function(xdata, ydata=NULL, minfreq=NULL, maxfreq=NULL, flag=2, window="
#
# OUT = list(Freq=myfreq, Amp=myamp, Phase=myphase)
# save(OUT, file=nameout)
-#}
+#}:
diff --git a/R/phase_randomize.R b/R/phase_randomize.R
index 6243cd5bcb1ac58f74fbdf6f3997b56c78eca170..1164e741e84051f6380b309b02712add0b412c46 100644
--- a/R/phase_randomize.R
+++ b/R/phase_randomize.R
@@ -1,5 +1,12 @@
-## randomize the phases of a real signal, for non-linear time series analysis
-
+#' Phase Randomizer
+#'
+#' randomizes the phases of a real signal, for non-linear time series analysis
+#' It works by computing the FFT of the signal, then randomize the phases of all
+#' complex numbers, but preserving the complex conjugate relationship guarantee
+#' that the inverse transform will be real
+#' @param x the input time series
+#' @return x the output time series.
+#'
phase_randomize <- function (x)
{
N <- length(x)
diff --git a/R/ssa.R b/R/ssa.R
index 7830e7585c70eb985294debbefc454439094ebc2..d08d2d54b74fdd29a76c1fcca620a2c3b04eebfd 100644
--- a/R/ssa.R
+++ b/R/ssa.R
@@ -8,8 +8,70 @@
#----------------------------------------------------------------------
# using astronomical forcing as sample data
-ssa <- function(X=X,N=length(X),M=M,Nrec=10)
-{
+#' Singular spectrum analysis
+#'
+#' @description Singular spectrum analysis as referred and described in Ghil et al.
+#' and in the \url{http://research.atmos.ucla.edu/tcd/ssa/}
+#' results have been validated and shown to produce exactly the same results
+#' as the Spectra toolking on 25.9.08
+
+#' @param X the input series (typically a numeric)
+#' @param N if provided, the length of the input series to beconsiderd
+#' @param M the length of the sliding window
+#' @param Nrec: number of components to be output
+#' @references
+#' \insertRef{Ghil02aa}{gtseries}
+#' \insertRef{VAUTARD89aa}{gtseries}
+#' @return a SSA object, compatible with the supplied plot function
+#' @examples
+#'
+#' bbridge <- function (n=256, amp) {
+#' rw1 <- cumsum(rnorm(n,1))*amp / (n)
+#' rw2 <- cumsum(rnorm(n,1))*amp / (n)
+#' rw2 <- rw2 * rw1[n] / rw2[n]
+#' bbridge <- rw1 - rw2
+#' }
+#'
+#' N <- 256
+#'
+#' Phase1 <- bbridge(N,2*pi)
+#' Phase2 <- bbridge(N,2*pi)
+#' Amp1 <- 1 + bbridge(N,0.5)
+#' Amp2 <- 1 + bbridge(N,0.5)
+#'
+#' t = seq(N)
+#'
+#' P1 = 40
+#' P2 = 90
+#'
+#' sig1 <- sin(2*pi*t/P1 + Phase1)*Amp1
+#' sig2 <- sin(2*pi*t/P2 + Phase2)*Amp2
+#'
+#' noise <- bbridge(N, 0.3) + rnorm(N,0.2)
+#'
+#' signal <- noise + sig1 + sig2
+#'
+#' plot (signal, type='l')
+#'
+#' M <- 50
+#'
+#' SSAsignal <- ssa (signal, M=M, Nrec=M)
+#'
+#' NRec <- c(10,30,M)
+#'
+#' for (i in seq(along=NRec)) {
+#' rec <- apply(SSAsignal$PCA[seq(NRec[i]),], 2, sum)
+#' lines(rec,col=i+1)
+#' }
+#'
+#' # full reconstruction (modulo the mean of the signal)
+#' fullrec <- apply(SSAsignal$PCA[seq(M),], 2, sum)
+#'
+#' # show that the difference is small
+#' plot(fullrec - signal + mean(signal))
+#'
+#' @export ssa
+ssa <- function(X=X,N=length(X),M=M,Nrec=10) {
# construct Toeplitz Matrix for reference signal (eq. 6)
# according to Vautard and Ghil, 1989
@@ -90,6 +152,8 @@ ssa <- function(X=X,N=length(X),M=M,Nrec=10)
+#' @rdname ssa
+#' @export
plot.SSAObject <- function (SSAObject,...)
{
yat <- outer(seq(1:Nrec),10^(seq(-3,5)),"*")
diff --git a/inst/REFERENCES.bib b/inst/REFERENCES.bib
new file mode 100644
index 0000000000000000000000000000000000000000..1fff25de3b343ede86e5856e9188fd477c1774a1
--- /dev/null
+++ b/inst/REFERENCES.bib
@@ -0,0 +1,35 @@
+@article{Ghil02aa,
+ author = {Ghil, M. and Allen, M. R. and Dettinger, M. D. and Ide, K. and Kondrashov, D. and Mann, M. E. and Robertson, A. W. and Saunders, A. and Tian, Y. and Varadi, F. and Yiou, P. },
+ editor = {},
+ year = {2002},
+ pages = {769-794},
+ doi = {10.1029/2000RG000092},
+ abstract = {The analysis of univariate or multivariate time series provides crucial information to describe, understand, and predict climatic variability. The discovery and implementation of a number of novel methods for extracting useful information from time series has recently revitalized this classical field of study. Considerable progress has also been made in interpreting the information so obtained in terms of dynamical systems theory. In this review we describe the connections between time series analysis and nonlinear dynamics, discuss signal-to-noise enhancement, and present some of the novel methods for spectral analysis. The various steps, as well as the advantages and disadvantages of these methods, are illustrated by their application to an important climatic time series, the Southern Oscillation Index. This index captures major features of interannual climate variability and is used extensively in its prediction. Regional and global sea surface temperature data sets are used to illustrate multivariate spectral methods. Open questions and further prospects conclude the review.},
+ keywords = {time-series analysis, SSA, dynamical systems},
+ title = {Advanced spectral methods for climatic time series},
+ volume = {40},
+ journal = {Reviews of Geophysics},
+}
+@article{VAUTARD89aa,
+ author = {Vautard, R. and Ghil, M. },
+ editor = {},
+ year = {1989},
+ pages = {395-424},
+ date = {MAY 1989},
+ title = {Singular spectrum analysis in nonlinear dynamics, with applications to paleoclimatic time-series},
+ volume = {35},
+ journal = {Physica D},
+ number = {3},
+}
+@article{sidlichovsky97aa,
+ author = {Šidlichovský, M. and Nesvorný, D. },
+ editor = {},
+ year = {1997},
+ title = {Frequency Modified Fourier Transform and its Application to Asteroids},
+ pages = {137–148},
+ journal = {The Dynamical Behaviour of our Planetary System},
+ doi = {10.1007/978-94-011-5510-6_9},
+ url = {http://dx.doi.org/10.1007/978-94-011-5510-6_9},
+ publisher = {Springer Netherlands},
+ isbn = {9789401155106},
+}
diff --git a/man/approx_ts.Rd b/man/approx_ts.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..98c135c248a6dd245aa8f2728d45b44963cd47ce
--- /dev/null
+++ b/man/approx_ts.Rd
@@ -0,0 +1,45 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/approx_ts.R
+\name{approx_ts}
+\alias{approx_ts}
+\title{Produce an evenly-spaced time series by simple interpolation}
+\usage{
+approx_ts(
+ data,
+ tcoord = "x",
+ dcoord = "y",
+ scale = 1,
+ n = 2048,
+ thin = FALSE,
+ spline = FALSE,
+ xlim = NULL
+)
+}
+\arguments{
+\item{data}{The data input, typically a dataframe,
+with two columns}
+
+\item{tcoord}{The name of the column name of the time coordinate. Defaults to "x"}
+
+\item{dcoord}{The name of the column name of the data coordinate. Defaults to "y"}
+
+\item{scale}{Time scaling (the coordinate "tcoord" will be divided by `scale`.}
+
+\item{n}{(defaults to 2048) is the number of data points in the output}
+
+\item{thin}{Discards data that are sampled with a time space less than 0.5 times the time resolution in the output}
+
+\item{spline}{: if True, the standard `spline` routine will be used, otherwise, `approx'}
+
+\item{xlim}{: if provided, brackets the time output}
+}
+\value{
+a time series R object
+}
+\description{
+Uses standard R routine 'approx' or 'spline' to generate an evenly-spaced time series
+based on data supplied on a time coordinate grid.
+}
+\author{
+Michel Crucifix
+}
diff --git a/man/arspec.Rd b/man/arspec.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..8ff82b8a1e6f78ecda443e3624387d65fc61c1dc
--- /dev/null
+++ b/man/arspec.Rd
@@ -0,0 +1,13 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/arspec.R
+\name{arspec}
+\alias{arspec}
+\title{Theoretical Spectrum of autoregressive spectrum}
+\usage{
+arspec(f, ar)
+}
+\description{
+Calculates the theoretical density frequency of an
+autoregressive process, according to equation
+(3.3) of Akaike (1969)
+}
diff --git a/man/cwt_morlet.Rd b/man/cwt_morlet.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..c0833e02dea96b2a4ddb552099745874019965e6
--- /dev/null
+++ b/man/cwt_morlet.Rd
@@ -0,0 +1,38 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/cwt_morlet.R
+\name{cwt_morlet}
+\alias{cwt_morlet}
+\alias{plot.wavelet}
+\alias{powerspectrum.wavelet}
+\title{Continous Morlet Wavelet Transform}
+\usage{
+cwt_morlet(
+ A,
+ inter = 20,
+ k0 = 5.6,
+ amin = 1,
+ amax = Inf,
+ calcmask = TRUE,
+ scale = NA,
+ deriv = FALSE
+)
+
+\method{plot}{wavelet}(
+ wave,
+ resx = 400,
+ resy = 300,
+ xlab = "Time",
+ ylab = "Period",
+ scaling_correction = 0,
+ col = col_wavelet,
+ legend = FALSE,
+ Mode = Mod,
+ plotMask = TRUE,
+ ...
+)
+
+powerspectrum.wavelet(wave, ...)
+}
+\description{
+Continous Morlet Wavelet Transform
+}
diff --git a/man/mem.Rd b/man/mem.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..2dde547cd4935fa0d936d99c3be442815239fde1
--- /dev/null
+++ b/man/mem.Rd
@@ -0,0 +1,11 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/mem.R
+\name{mem}
+\alias{mem}
+\title{Maximum entropy method}
+\usage{
+mem(A, method = "burg", order.max = 90, ...)
+}
+\description{
+Maximum entropy method
+}
diff --git a/man/mfft.Rd b/man/mfft.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..f547f6e31dbc060a3851ee2daf3e6beb72d77e0d
--- /dev/null
+++ b/man/mfft.Rd
@@ -0,0 +1,41 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/mfft.R
+\name{mfft}
+\alias{mfft}
+\title{Modified Fourier transform}
+\usage{
+mfft(xdata, minfreq = NULL, maxfreq = NULL, flag = 1, nfreq = 30)
+}
+\arguments{
+\item{xdata}{The data provided either as a time series (advised), or as a vector.
+may be complex}
+
+\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{flag}{:
+ Modified Fourier Transform if flag = 1;
+ Frequency Modified Fourier Transform if flag = 2;
+ FMFT with additional non-linear correction if flag = 3
+(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.}
+
+\item{nfreq}{is the number of frequencies returned, must be smaller that the length of xdata.}
+}
+\value{
+STILL NEED TO BE DESCRIBED
+}
+\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/phase_randomize.Rd b/man/phase_randomize.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..e27a947be5097d93774c82b37c82cb2cfd3e19ce
--- /dev/null
+++ b/man/phase_randomize.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/phase_randomize.R
+\name{phase_randomize}
+\alias{phase_randomize}
+\title{Phase Randomizer}
+\usage{
+phase_randomize(x)
+}
+\arguments{
+\item{x}{the input time series}
+}
+\value{
+x the output time series.
+}
+\description{
+randomizes the phases of a real signal, for non-linear time series analysis
+It works by computing the FFT of the signal, then randomize the phases of all
+complex numbers, but preserving the complex conjugate relationship guarantee
+that the inverse transform will be real
+}
diff --git a/man/ssa.Rd b/man/ssa.Rd
index 8a389fafa4d125f838ae337bb5dfe5d6a91b322e..858a583a34f1333031fb3c2456ceb43266447b25 100644
--- a/man/ssa.Rd
+++ b/man/ssa.Rd
@@ -1,31 +1,33 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/ssa.R
\name{ssa}
\alias{ssa}
-\title{Singar Spectrum Analysis}
-\description{
- Singular spectrum analysis based on Ghil et al. 2002}
-
+\alias{plot.SSAObject}
+\title{Singular spectrum analysis}
\usage{
- (X = X, N = length(X), M = M, Nrec = 10)
+ssa(X = X, N = length(X), M = M, Nrec = 10)
+
+\method{plot}{SSAObject}(SSAObject, ...)
}
\arguments{
- \item{X}{numeric vector}
- \item{N}{number of points to be considered in the vector}
- \item{M}{window length}
- \item{Nrec}{number of components to keep for the reconstruction}
+\item{X}{the input series (typically a numeric)}
+
+\item{N}{if provided, the length of the input series to beconsiderd}
+
+\item{M}{the length of the sliding window}
+
+\item{Nrec:}{number of components to be output}
}
-\details{}
\value{
- an \code{SSAObject}
+a SSA object, compatible with the supplied plot function
}
-\author{Michel Crucifix}
-\references{
-Yiou P., E. Baert and M. F. Loutre (1996), Spectral analysis of climate data, Surveys in Geophysics, (17) 619-663 doi:10.1007/BF01931784
-
-Ghil M., M. R. Allen, M. D. Dettinger, K. Ide, D. Kondrashov, M. E. Mann, A. W. Robertson, A. Saunders, Y. Tian, F. Varadi and P. Yiou (2002), Advanced spectral methods for climatic time series, Reviews of Geophysics, (40) 769-794 doi:10.1029/2000RG000092
+\description{
+Singular spectrum analysis as referred and described in Ghil et al.
+and in the \url{http://research.atmos.ucla.edu/tcd/ssa/}
+results have been validated and shown to produce exactly the same results
+as the Spectra toolking on 25.9.08
}
-
\examples{
-require (gtseries)
bbridge <- function (n=256, amp) {
rw1 <- cumsum(rnorm(n,1))*amp / (n)
@@ -69,6 +71,11 @@ lines(rec,col=i+1)
# full reconstruction (modulo the mean of the signal)
fullrec <- apply(SSAsignal$PCA[seq(M),], 2, sum)
-# show that the difference is smal
+# show that the difference is small
plot(fullrec - signal + mean(signal))
+
+}
+\references{
+\insertRef{Ghil02aa}{gtseries}
+\insertRef{VAUTARD89aa}{gtseries}
}
diff --git a/old_man/ssa.Rd b/old_man/ssa.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..8a389fafa4d125f838ae337bb5dfe5d6a91b322e
--- /dev/null
+++ b/old_man/ssa.Rd
@@ -0,0 +1,74 @@
+\name{ssa}
+\alias{ssa}
+\title{Singar Spectrum Analysis}
+\description{
+ Singular spectrum analysis based on Ghil et al. 2002}
+
+\usage{
+ (X = X, N = length(X), M = M, Nrec = 10)
+}
+\arguments{
+ \item{X}{numeric vector}
+ \item{N}{number of points to be considered in the vector}
+ \item{M}{window length}
+ \item{Nrec}{number of components to keep for the reconstruction}
+}
+\details{}
+\value{
+ an \code{SSAObject}
+}
+\author{Michel Crucifix}
+\references{
+Yiou P., E. Baert and M. F. Loutre (1996), Spectral analysis of climate data, Surveys in Geophysics, (17) 619-663 doi:10.1007/BF01931784
+
+Ghil M., M. R. Allen, M. D. Dettinger, K. Ide, D. Kondrashov, M. E. Mann, A. W. Robertson, A. Saunders, Y. Tian, F. Varadi and P. Yiou (2002), Advanced spectral methods for climatic time series, Reviews of Geophysics, (40) 769-794 doi:10.1029/2000RG000092
+}
+
+\examples{
+require (gtseries)
+
+bbridge <- function (n=256, amp) {
+ rw1 <- cumsum(rnorm(n,1))*amp / (n)
+ rw2 <- cumsum(rnorm(n,1))*amp / (n)
+ rw2 <- rw2 * rw1[n] / rw2[n]
+ bbridge <- rw1 - rw2
+}
+
+N <- 256
+
+Phase1 <- bbridge(N,2*pi)
+Phase2 <- bbridge(N,2*pi)
+Amp1 <- 1 + bbridge(N,0.5)
+Amp2 <- 1 + bbridge(N,0.5)
+
+t = seq(N)
+
+P1 = 40
+P2 = 90
+
+sig1 <- sin(2*pi*t/P1 + Phase1)*Amp1
+sig2 <- sin(2*pi*t/P2 + Phase2)*Amp2
+
+noise <- bbridge(N, 0.3) + rnorm(N,0.2)
+
+signal <- noise + sig1 + sig2
+
+plot (signal, type='l')
+
+M <- 50
+
+SSAsignal <- ssa (signal, M=M, Nrec=M)
+
+NRec <- c(10,30,M)
+
+for (i in seq(along=NRec)) {
+rec <- apply(SSAsignal$PCA[seq(NRec[i]),], 2, sum)
+lines(rec,col=i+1)
+}
+
+# full reconstruction (modulo the mean of the signal)
+fullrec <- apply(SSAsignal$PCA[seq(M),], 2, sum)
+
+# show that the difference is smal
+plot(fullrec - signal + mean(signal))
+}
diff --git a/src/fmft.c b/src/fmft.c
index ccafefa5474c546434527af0a11b0364da459ccf..d7e9cb30e6c5c7d4e931b27f2e6d399831d30b7d 100644
--- a/src/fmft.c
+++ b/src/fmft.c
@@ -6,9 +6,21 @@ estimates the frequencies (f_j), amplitudes (A_j) and phases
X(t) + iY(t) = Sum_j=1^N [ A_j * exp i (f_j * t + psi_j) ] */
+/* adapted from https://www.boulder.swri.edu/~davidn/fmft/fmft.html */
+/* most of the code has been authored by David Nesvornky and released on the above website */
+/* adaptations by Michel Crucifix (2023) for
+ * - avoid references to Numerical Recipies: replace Fourier transform and sorting by corresponding gsl routines
+ * - generalise to data length not a power of 2
+ * - technical adaptations for compatibility with R */
+
+/* left overs of development phase */
+/*#define N2FLAG */
+/*#define N2FLAG2*/
+
#define FMFT_TOL 1.0e-10 /* MFT NOMINAL PRECISION */
#define FMFT_NEAR 0. /* MFT OVERLAP EXCLUSION PARAMETER */
+#include <stdbool.h>
#include <stdio.h>
#include <math.h>
@@ -22,20 +34,28 @@ X(t) + iY(t) = Sum_j=1^N [ A_j * exp i (f_j * t + psi_j) ] */
#define TWOPI (2.*PI)
-static int itemp;
-static unsigned long ultemp;
static float ftemp;
static double dtemp;
-#define FSQR(a) ((ftemp=(a)) == 0.0 ? 0.0 : ftemp*ftemp)
#define DSQR(a) ((dtemp=(a)) == 0.0 ? 0.0 : dtemp*dtemp)
-
#define SHFT3(a,b,c) (a)=(b);(b)=(c)
#define SHFT4(a,b,c,d) (a)=(b);(b)=(c);(c)=(d)
-#define ISWAP(a,b) itemp=(a);(a)=(b);(b)=itemp
-#define ULSWAP(a,b) ultemp=(a);(a)=(b);(b)=ultemp
-#define FSWAP(a,b) ftemp=(a);(a)=(b);(b)=ftemp
+bool fastflag;
+
+bool isPowerofTwo (size_t n){
+ if (n == 0)
+ return 0;
+ while (n != 1){
+ if (n % 2 != 0)
+ return (0);
+ n = n/2;
+ }
+ return 1;
+}
+
+
+
void window(double *x, double *y, double *xdata, double *ydata, size_t ndata);
@@ -69,6 +89,7 @@ void dindex(unsigned long n, double arr[], unsigned long indx[]);
};
+
int fmft(int *localnfreq, double *localminfreq, double *localmaxfreq, int *localflag,
int *localndata, double *localxdata, double *localydata,
struct component *signal1, struct component *signal2, struct component *signal3);
@@ -124,10 +145,19 @@ be a power of 2), are the input data X(j-1) and Y(j-1).
int flag = *localflag;
size_t ndata = *localndata;
+
+ fastflag = isPowerofTwo(ndata) ;
+
+ if (fastflag)
+ (printf("ndata is power of two: we will be faster ! \n"));
+ if (ndata <= 2){
+ printf("at least 2 data needed - output non-reliable"); return(0);
+ }
+ if (ndata <= nfreq){
+ printf("nfreq must be smaller than nata"); return(0);
+ }
+
/* printf("prelimarg %d, %zu %d", nfreq, ndata, flag);*/
- printf("prelimarg %d, %d %d", nfreq, *localndata, *localflag);
- printf("xdata: %.2f, %.2f, %.2f", localxdata[0],localxdata[1],localxdata[2]);
- printf("ydata: %.2f, %.2f, %.2f", localydata[0],localydata[1],localydata[2]);
/* ALLOCATION OF VARIABLES */
@@ -502,11 +532,33 @@ CALLS FFT AND RETURNS POWER SPECTRAL DENSITY */
z[2*j-1] = y[j];
}
+// #ifdef N2FLAG2
+
+if (fastflag){
+
gsl_fft_complex_radix2_backward (z,1,ndata);
-
+
+} else {
+
+ gsl_fft_complex_wavetable * wavetable;
+ gsl_fft_complex_workspace * workspace;
+
+ wavetable = gsl_fft_complex_wavetable_alloc (ndata);
+ workspace = gsl_fft_complex_workspace_alloc (ndata);
+
+
+ gsl_fft_complex_backward (z, 1, ndata,
+ wavetable, workspace);
+ gsl_fft_complex_wavetable_free (wavetable);
+ gsl_fft_complex_workspace_free (workspace);
+
+}
+
for(j=1;j<=ndata;j++)
- powsd[j] = FSQR(z[2*j-2]) + FSQR(z[2*j-1]);
+ powsd[j] = DSQR(z[2*j-2]) + DSQR(z[2*j-1]);
+
+/* return(0);*/
}
double bracket(float *powsd, size_t ndata)
@@ -541,14 +593,14 @@ double bracket(float *powsd, size_t ndata)
maxpow = powsd[1];
}
- if(maxpow == 0) nrerror("DFT has no maximum ...");
+ if(maxpow == 0) printf("DFT has no maximum ...");
- if(maxj < ndata/2) freq = -(maxj-1);
- if(maxj > ndata/2) freq = -(maxj-ndata-1);
+ if(maxj < ndata/2 - 1) freq = -(maxj-1);
+ if(maxj > ndata/2 - 1) freq = -(maxj-ndata-1);
return (TWOPI*freq / ndata);
- /* negative signs and TWOPI compensate for the Numerical Recipes
+ /* negative signs and TWOPI compensate for the
definition of the DFT */
}
@@ -622,13 +674,16 @@ double phisqr(double freq, double xdata[], double ydata[], size_t ndata)
}
void phifun(double *xphi, double *yphi, double freq,
- double xdata[], double ydata[], long n)
+ double xdata[], double ydata[], long n){
/* COMPUTES THE FUNCTION PHI */
+ /* THIS REQUIRES DATA LENGTH TO BE 2^N */
+ /* BECAUSE IT PERFORMS A DYADIC DECOMPOSITION */
-{
- long i, j, nn;
- double c, s, *xdata2, *ydata2;
+ long i;
+ double c, s;
+
+ double *xdata2, *ydata2;
xdata2 = dvector(1, n);
ydata2 = dvector(1, n);
@@ -641,6 +696,12 @@ void phifun(double *xphi, double *yphi, double freq,
ydata2[i] = ydata[i];
}
+
+/*#ifdef N2FLAG*/
+
+if (fastflag) {
+
+ long j, nn;
nn = n;
while(nn != 1){
@@ -657,12 +718,33 @@ void phifun(double *xphi, double *yphi, double freq,
}
}
-
+
*xphi = 2*xdata2[1] / (n-1);
*yphi = 2*ydata2[1] / (n-1);
- free_dvector(xdata2,1,n);
- free_dvector(ydata2,1,n);
+
+
+} else {
+
+/*xdata2[1] = xdata[1] ; ydata2[1] = ydata[1] ;*/
+/*xdata2[n] = xdata[n] ; ydata2[n] = ydata[n] ;*/
+
+ *xphi = 0;
+ *yphi = 0;
+
+ for(i=1;i<=n;i++){
+
+ c = 2*cos(-(i-1)*freq)/(n-1);
+ s = 2*sin(-(i-1)*freq)/(n-1);
+
+ *xphi += c*xdata2[i] - s*ydata2[i];
+ *yphi += c*ydata2[i] + s*xdata2[i];
+ }
+
+}
+
+free_dvector(xdata2,1,n);
+free_dvector(ydata2,1,n);
}