oxygenise + doc + mfft fixes

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

NAMESPACE

-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)

NAMESPACE.old

+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)

R/approx_ts.R

-  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]] 

R/arspec.R

-## 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

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")

R/mem.R

+
+#' Maximum entropy method 
+
+#' @export mem
 mem <- function (A,method='burg',order.max=90,...)
 {
   deltat<-deltat(A)

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)
-#}
+#}:

R/phase_randomize.R

-## 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)

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)),"*")

inst/REFERENCES.bib

+@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}, 
+} 

man/approx_ts.Rd

+% 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
+}

man/arspec.Rd

+% 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)
+}

man/cwt_morlet.Rd

+% 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
+}

man/mem.Rd

+% 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
+}

man/mfft.Rd

+% 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}
+}

man/phase_randomize.Rd

+% 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
+}

man/ssa.Rd

+% 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}
 }

old_man/ssa.Rd

+\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))
+}

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);
 
 }