add mfft

.Rbuildignore

+.gitignore
+test
+man_old
+^data-raw$
+.Rprofile
+TODO
+^LICENSE\.md$
+.Renviron

.Renviron

+# needed for WSL
+TZ="Europe/Brussels"

.Rprofile

+library(devtools)
+library(usethis)
+library(roxygen2)
+# load_all()

.gitignore

+.Rproj.user
+.Rhistory
+.Rdata
+.httr-oauth
+.DS_Store
+*.o
+*.so
+*.pdf
+config.log
+config.status

DESCRIPTION

@@ -3,7 +3,7 @@ Type: Package
 Title: Time Series analysis for galcial cycles
 Version: 1.04
 Date: 2022-05-28
-Author: Michel Crucifix
+Author: Michel Crucifix [aut, cre] (<https://orcid.org/0000-0002-3437-4911>)
 Maintainer: Michel Crucifix <michel.crucifix@uclouvain.be>
 Description: More about what it does (maybe more than one line)
 License: What license is it under?

NAMESPACE

-exportPattern("^[[:alpha:]]+")
+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/mfft.R

+# I need to add lots of comments about licencing
+# the fact that the hanning window is now made here
+# 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)
+{
+
+  # please at the moent do not play with 'flag'
+  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 (is.null(minfreq)){
+      my_minfreq = pi/ndata
+      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)
+
+
+     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)
+
+     Freq <- OUTVEC[seq(nfreq) ]/dt
+     Ampl <- OUTVEC[nfreq+seq(nfreq) ] * N/ndata
+     Phase <- OUTVEC[2*nfreq+seq(nfreq) ]
+
+     # 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(OUTVEC)
+# dyn.load('fmft.so')
+# ndata=8192
+
+#for (exp in c('d')){
+#  namein <- file.path('La2010',sprintf("La2010%s_alkhqp3L.dat",exp));
+#  nameout <- sprintf("La2010%s_out.RData",exp);
+#  A <- read.table(namein)
+#  deltat=1000
+#  colnames(A) <- c('times','a','l','k','h','q','p')
+#  # ndata = as.integer64(ndata)
+#
+#  times <- seq(ndata)
+#
+#  nfreq <- as.integer(30)
+#  flag =  2 
+#
+#  minfreq <- as.double ( -1e6 / 180. / 3600. * pi /5. ) ;
+#  maxfreq <- as.double ( 1.e6 / 180. / 3600. * pi /5. ) ;
+#
+#  signal1 = as.double (rep(0, nfreq*3))
+#  signal2 = as.double (rep(0, nfreq*3))
+#  signal3 = as.double (rep(0, nfreq*3))
+#
+#  ntrunks <- floor(length(A$times)/ndata)
+#
+#  # EMFFT <- sapply (seq(0,floor(length(A$times)/ndata)-1), function(i)
+#
+#
+#  timemin <- sapply (seq(0,ntrunks-1), function(i) { A$times [ndata* i + 1] });
+#  timemax <- sapply (seq(0,ntrunks-1), function(i) { A$times [ndata* i + ndata ] });
+#
+#  EMFFT <- sapply (seq(0,ntrunks-1), function(i)
+#                       {
+#                         xdata <- A$h[ndata*i + (1:ndata)]
+#                         ydata <- A$k[ndata*i + (1:ndata)]
+#
+#                         OUT = .C("fmft", nfreq, minfreq, maxfreq, as.integer(flag), as.integer(ndata), xdata, ydata, signal1,signal2,signal3, DUP=FALSE)
+#
+#                         # OUTVEC <- matrix(OUT[[8]], 10,3)
+#                         OUTVEC <- t(matrix(OUT[[8]], 3, nfreq))
+#
+#                         return(OUTVEC)
+#
+#                       })
+#
+#
+#  # repack
+#
+#  # repack
+#
+#
+#  Freq <- EMFFT[seq(nfreq), ]*360*60*60/2/pi/deltat
+#  Ampl <- EMFFT[nfreq+seq(nfreq), ]
+#  Phas <- EMFFT[2*nfreq+seq(nfreq), ]
+#
+#
+#  snake <- function(x1, x2, tol=2.e-1, mynfreq=nfreq){
+#    xnew <- c();
+#    xamp <- c();
+#    xpha <- c();
+#    xf1  <- x1[seq(mynfreq)]
+#    xf2  <- x2[seq(mynfreq)]
+#    xa1  <- x1[mynfreq+seq(mynfreq)]
+#    xa2  <- x2[mynfreq+seq(mynfreq)]
+#    xp1  <- x1[2*mynfreq+seq(mynfreq)]
+#    xp2  <- x2[2*mynfreq+seq(mynfreq)]
+#
+#    while(length(xf1) > 0) {
+#      i <- which.min(abs(xf2 - xf1[1])+0*abs(xa2-xa1[1]));
+#      xnew <- c(xnew, xf2[i])
+#      xamp <- c(xamp, xa2[i])
+#      xpha <- c(xpha, xp2[i])
+#      xf2 <- xf2[-i]
+#      xf1 <- xf1[-1]
+#      xa2 <- xa2[-i]
+#      xa1 <- xa1[-1]
+#      xp2 <- xp2[-i]
+#      xp1 <- xp1[-1]
+#    }
+#    return(c(xnew, xamp, xpha))
+#  }
+#
+#  mylist <- apply(EMFFT,2,as.numeric, simplify=FALSE)
+#
+#  terms <- simplify2array(Reduce(snake, mylist, mylist[[1]], accumulate=TRUE))
+#
+#  myfreq <- terms[seq(nfreq),-1] * 360*60*60/2/pi/deltat
+#  myamp <- terms[nfreq+seq(nfreq),-1]
+#  myphase <- terms[2*nfreq+seq(nfreq),-1]
+#
+#
+#  rownames (myfreq) <- timemin
+#  rownames (myamp)  <- timemin
+#  rownames (myphase)<- timemin
+#
+#
+#  OUT = list(Freq=myfreq, Amp=myamp, Phase=myphase)
+#  save(OUT, file=nameout)
+#}

src/Makevars

+# Kindly supplied by Dirk Eddelbuettel
+# set by configure
+GSL_CFLAGS = -I/usr/local/include
+GSL_LIBS   = -L/usr/local/lib -lgsl -lgslcblas -lm
+
+# combine to standard arguments for R
+PKG_CPPFLAGS =  $(GSL_CFLAGS) -I.
+PKG_LIBS = $(GSL_LIBS)
+
+

src/Makevars.in

+# Kindly supplied by Dirk Eddelbuettel
+# set by configure
+GSL_CFLAGS = @GSL_CFLAGS@
+GSL_LIBS   = @GSL_LIBS@
+
+# combine to standard arguments for R
+PKG_CPPFLAGS =  $(GSL_CFLAGS) -I.
+PKG_LIBS = $(GSL_LIBS)
+
+

src/Makevars.win

+# PKG_LIBS=-LF:/MinGW/usr/local/lib -lgsl -lgslcblas
+# CPPFLAGS=-I$(R_HOME)/include -IF:/MinGW/usr/local/include
+# PKG_CPPFLAGS=-IF:/MinGW/usr/local/include
+
+# lines below supplied by Brian Ripley and Uwe Ligges
+
+PKG_CPPFLAGS=-I$(LIB_GSL)/include
+PKG_LIBS=-L$(LIB_GSL)/lib -lgsl -lgslcblas
+

src/nrutil.c

+/* NRUTIL: from "Numerical Recipes in C", pp. 705-709 */
+#include <stdio.h>
+#include <stdlib.h>
+#include <malloc.h>
+
+/* nrerror() := Numerical Recipes standard error handler */
+void nrerror(char *error_text)
+{
+	fprintf(stderr, "Numerical Recipes run-time error...\n");
+	fprintf(stderr, "%s\n", error_text);
+	fprintf(stderr, "...now exiting to system...\n");
+	exit(1);
+}
+
+/* *vector() := Allocates a float vector with range [nl..nh] */
+float *vector(int nl, int nh)
+{
+	float *v;
+
+	v = (float *) malloc( (unsigned)(nh-nl+1) * sizeof(float) );
+	if (!v)
+		nrerror("\nallocation failure in vector()");
+	return v-nl;
+}
+
+float Square(float x)
+{
+	return x * x;
+}
+
+int *ivector(int nl, int nh)
+{
+	int *v;
+	v = (int *)malloc((unsigned) (nh-nl+1)*sizeof(int));
+	if (!v) nrerror("\nallocation failure in ivector()");
+	return v-nl;
+}
+
+long *lvector(int nl, int nh)
+{
+	long *v;
+	v = (long *)malloc((unsigned) (nh-nl+1)*sizeof(long));
+	if (!v) nrerror("\nallocation failure in ivector()");
+	return v-nl;
+}
+
+/* *dvector() := Allocates a double vector with range [nl..nh] */
+double *dvector(int nl, int nh)
+{
+	double *v;
+
+	v = (double *) malloc( (unsigned)(nh-nl+1) * sizeof(double) );
+	if (!v)
+		nrerror("allocation failure in dvector()");
+	return v-nl;
+}
+
+double **dmatrix(int nrl, int nrh, int ncl, int nch)
+{
+	int i;
+	double **m;
+
+	m = (double **) malloc( (unsigned)(nrh - nrl + 1) * sizeof(double *) );
+	if (!m)
+		nrerror("allocation failure 1 in dmatrix()");
+	m -= nrl;
+
+	for (i = nrl; i <= nrh; i++)
+	{
+		m[i] = (double *) malloc( (unsigned)(nch - ncl + 1)*sizeof(double) );
+		if (!m[i])
+			nrerror("allocation failure 2 in dmatrix()");
+		m[i] -= ncl;
+	}
+	return m;
+}
+
+float **matrix(int nrl, int nrh, int ncl, int nch)
+{
+	int i;
+	float **m;
+
+	m = (float **) malloc( (unsigned)(nrh - nrl + 1) * sizeof(float *) );
+	if (!m)
+		nrerror("allocation failure 1 in matrix()");
+	m -= nrl;
+
+	for (i = nrl; i <= nrh; i++)
+	{
+		m[i] = (float *) malloc( (unsigned)(nch - ncl + 1)*sizeof(float) );
+		if (!m[i])
+			nrerror("allocation failure 2 in matrix()");
+		m[i] -= ncl;
+	}
+	return m;
+}
+
+/* free_vector() := Frees a float vector allocated by vector() */
+void free_vector(float *v, int nl, int nh)
+{
+	free( (char *)(v + nl) );
+}
+
+void free_ivector(int *v, int nl, int nh)
+{
+	free( (char *)(v + nl) );
+}
+
+void free_lvector(long *v, int nl, int nh)
+{
+	free( (char *)(v + nl) );
+}
+
+void free_dmatrix(double **m, int nrl, int nrh, int ncl, int nch)
+{
+	int i;
+
+	for (i = nrh; i >= nrl; i--)
+		free( (char *)(m[i] + ncl) );
+	free( (char *)(m + nrl) );
+}
+
+/* free_dvector() := Frees a double vector allocated by dvector() */
+void free_dvector(double *v, int nl, int nh)
+{
+	free( (char *)(v + nl) );
+}
+
+void free_matrix(float **m, int nrl, int nrh, int ncl, int nch)
+{
+	int i;
+
+	for (i = nrh; i >= nrl; i--)
+		free( (char *)(m[i] + ncl) );
+	free( (char *)(m + nrl) );
+}
+

src/nrutil.h

+extern void nrerror(char *error_text);
+extern float Square(float x);
+extern double *dvector(int nl, int nh);
+extern double **dmatrix(int nrl, int nrh, int ncl, int nch);
+extern float **matrix(int nrl, int nrh, int ncl, int nch);
+extern float *vector(int nl, int nh);
+extern int *ivector(int nl, int nh);
+extern long *lvector(int nl, int nh);
+extern void free_dvector(double *v, int nl, int nh);
+extern void free_dmatrix(double **m, int nrl, int nrh, int ncl, int nch);
+extern void free_vector(float *v, int nl, int nh);
+extern void free_ivector(int *v, int nl, int nh);
+extern void free_lvector(long *v, int nl, int nh);
+extern void free_matrix(float **m, int nrl, int nrh, int ncl, int nch);