NAMESPACE

@@ -16,6 +16,7 @@ export(mem)
 export(mfft_anova)
 export(mfft_complex)
 export(mfft_real)
+export(mfft_real_C)
 export(periodogram)
 export(powerspectrum.wavelet)
 export(reconstruct_mfft)

R/mfft_real.R

 
-Q <- function(wT) {sin(wT)/(wT)*(pi^2)/(pi^2-(wT)^2)}
-Qprime <- function(y) {ifelse(y==0, 0, pi^2/(pi^2-y^2)/y*(cos(y)+(sin(y)/y)*(3*y^2-pi^2)/(pi^2-y^2)))}
-Qsecond0 <- 2/pi^2 - 1./3. 
+pisquare <- pi^2
+
+Q <- function(wT) {sin(wT)/(wT)*(pisquare)/(pisquare-(wT*wT))}
+Qprime <- function(y) {ifelse(y==0, 0, pisquare/(pisquare-y*y)/y*(cos(y)+(sin(y)/y)*(3*y*y-pisquare)/(pisquare-y*y)))}
+Qsecond0 <- 2/pisquare - 1./3. 
 
 mfft_analyse <- function(xdata, nfreq, fast = TRUE, nu = NULL, minfreq=NULL, maxfreq=NULL){
   
@@ -63,27 +65,48 @@ mfft_analyse <- function(xdata, nfreq, fast = TRUE, nu = NULL, minfreq=NULL, max
     # or either they were not provided, but we are in this case where it was set in "m-1" because
     # we identified a non-null frequency
     # in both configurations, we look at a frequency with the fourier transform
-    fbase <- freqs[which.max(power(hx))]
-    brackets <- c(fbase-pi/N, fbase+pi/N);
+     fbase <- freqs[which.max(power(hx))]
+     brackets <- c(fbase-pi/N, fbase+pi/N);
     # and then further corrects the brackets if minfreq and maxfreq where provided
-    brackets[1] <- max(minfreq, brackets[1])
-    brackets[2] <- max(maxfreq, brackets[1])
-    thx <- t(hx)
+     brackets[1] <- max(minfreq, brackets[1])
+     brackets[2] <- min(maxfreq, brackets[2])
+     thx <- t(hx)
     
     # after profiling, the fastest seems the first option below
     # this is the weak link
     # coding the function in c might be an option
     
-    fmax = cmna::goldsectmax(function(f) 
-                             {
-                              ft <- f*t
-                             (thx %*% cos(f*t))^2 + (thx %*% sin(f*t))^2}, 
-                             brackets[1], brackets[2], tol=1.e-10, m=9999)
-    
-   #  fmax = cmna::goldsectmax(function(f) {
-   #                           ft <- f*t
-   #                           (sum(hx * cos(ft)))^2 + (sum(hx %*% sin(ft)))^2}, 
-   #                           brackets[1], brackets[2], tol=1.e-10, m=9999)
+    #
+     tomax <-  function(t) {
+       function(f) {
+       ft <- f*t
+       a <- hx %*% cbind(cos(ft), sin(ft))
+       a[1]*a[1] + a[2]*a[2]
+     } 
+     }
+
+
+#     print ("x[[m]][1,2]")
+#     print (as.double(x[[m]][1:2]))
+
+    localxdata <- as.double(x[[m]])
+#     print ('call it')
+#    OUT <- .C("fastgsec", as.double(minfreq), as.double(maxfreq), 
+#              as.integer(N), localxdata , as.double(rep(0,N)), 
+#              outfreq = 0., DUP=TRUE)
+#     print ('called it')
+
+#     plot (Mod(fft(x[[m]])[0:(N/2)])^2 , type='l')
+
+   fmax = cmna::goldsectmax(tomax(t), 
+                            brackets[1], brackets[2],
+                            tol=1.e-10, m=9999)
+   
+
+    # fmax = cmna::goldsectmax(function(f) {
+    #                          ft <- f*t
+    #                          (sum(hx * cos(ft)))^2 + (sum(hx %*% sin(ft)))^2}, 
+    #                          brackets[1], brackets[2], tol=1.e-10, m=9999)
 
    #  fmax = cmna::goldsectmax(function(f) {
    #                           ft <- f*t
@@ -91,11 +114,15 @@ mfft_analyse <- function(xdata, nfreq, fast = TRUE, nu = NULL, minfreq=NULL, max
    #                           brackets[1], brackets[2], tol=1.e-10, m=9999)
 
 
+    #print (sprintf("fmaxlocal: %.2f ; fmax with C: %.2f", 
+    #               fmax, OUT$outfreq))
+
+#    fmax <- OUT$outfreq
 
     if (fmax > freqs[2]/2){
       # if we really identified a frequency
       # we are in this case where the frequency was not a priori provided
-      # we se set this phase to zero, and the next one to pi/2, and also set the frequencies accordinly
+      # we se set this phase to zero, and the next one to pi/2, and also set the frequencies accordingly
     phase[m] <- 0.
     nu[m]  <- fmax
     phase[m+1] <- pi/2. 
@@ -190,13 +217,14 @@ mfft_analyse <- function(xdata, nfreq, fast = TRUE, nu = NULL, minfreq=NULL, max
 
     S[m] = 0. 
     for (j in seq(m))  S[m] = S[m] + Prod[j]*A[m,j]
+
     # for (j in seq(m))  S[m] = S[m] + A[m,j] * hprod(x[[1]], f[[j]])
 
     # les Sm sont les projections dans la nouvelle base
     # not necessary, for verification only
     # computes the B and verify they are orthonormal
-     B[[m]]=0
-     for (j in seq(m)) B[[m]] = B[[m]] + A[m,j] * f[[j]]
+    # B[[m]]=0
+    # for (j in seq(m)) B[[m]] = B[[m]] + A[m,j] * f[[j]]
 
 
     #print ('S[m] computed two different ways')
@@ -217,8 +245,8 @@ mfft_analyse <- function(xdata, nfreq, fast = TRUE, nu = NULL, minfreq=NULL, max
 
 
     x[[m+1]] = x[[m]]
-    # for (j in seq(m))  x[[m+1]] = x[[m+1]] - S[m] * A[m,j]*f[[j]]
-    x[[m+1]] = x[[m+1]] - S[m] * B[[m]]
+    for (j in seq(m))  x[[m+1]] = x[[m+1]] - S[m] * A[m,j]*f[[j]]
+    #x[[m+1]] = x[[m+1]] - S[m] * B[[m]]
 #     print('residu final')
 #     print(m)
 #     print (x[[m+1]][seq(20)])

man/mfft.Rd

@@ -4,12 +4,14 @@
 \alias{mfft}
 \title{Modified Fourier transform}
 \usage{
-mfft(xdata, minfreq = NULL, maxfreq = NULL, correction = 1, nfreq = 30)
+mfft(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. 
@@ -21,8 +23,6 @@ The parameter is passed on to `mfft_real` or `mfft_complex` depending on
 quite the same for both. This needs to be fixed. 
 one approach would be to define scalars like 'fft', 'mfft', 'mfft_linear', 
 'mfft_second_order_correction'.}
-
-\item{nfreq}{is the number of frequencies returned, must be smaller that the length of  xdata.}
 }
 \value{
 a `mfft_deco` object, based on a data.frame with columns "Freq", "Ampl" and "Phases". 

src/fmft.c

@@ -681,7 +681,7 @@ double phisqr(double freq, double xdata[], double ydata[], size_t ndata)
   xphi = yphi = 0;
   
   phifun(&xphi, &yphi, freq, xdata, ydata, ndata);
-  
+  // printf ("xphi etc %.f %.f %.f %.f \n',  xphi, yphi, xphi*xphi + yphi*yphi"); //
   return xphi*xphi + yphi*yphi;
 }
 
@@ -1337,9 +1337,11 @@ double bracket_real(float *powsd, size_t ndata)
   maxj = 0;
   maxpow = 0;
  
+/*  printf("I am calling bracket real  %i \n", ndata);*/
 /* probe positive frequencies in priority */
 
   for(j=2;j<=ndata/2-2;j++)
+/*      printf ("powsd[%d]  %.g  \n ", j, powsd[j]);*/
     if(powsd[j] > powsd[j-1] && powsd[j] > powsd[j+1])
       if(powsd[j] > maxpow){ 
 	maxj = j;
@@ -1356,7 +1358,7 @@ double bracket_real(float *powsd, size_t ndata)
   if(maxpow == 0) printf("DFT has no maximum ...");
  
    freq = (maxj-1);
-   printf("freq = %i \n",freq);
+/*   printf("freq = %d \n",freq);*/
 /*  if(maxj < ndata/2 - 1) freq = -(maxj-1);  */
 /*  if(maxj > ndata/2 - 1) freq = -(maxj-ndata-1);*/
   return (TWOPI*freq / ndata);
@@ -1505,3 +1507,104 @@ free_dvector(xdata2,1,n);
 
 }
 
+int fastgsec(double *localminfreq, double *localmaxfreq,  
+	 int *localndata, double *localxdata, double *localydata,  double *outfreq);
+
+int fastgsec(double *localminfreq, double *localmaxfreq,  
+	 int *localndata, double *localxdata, double *localydata, 
+   double *outfreq)  
+{
+  int nearfreqflag;
+  long j;
+  float *powsd;
+  double *xdata, *ydata, *x, *y;
+  double centerf, leftf, rightf;
+  double f, A, psi;
+
+  FILE *fp;
+
+  double minfreq = *localminfreq;
+  double maxfreq = *localmaxfreq;
+  size_t ndata = *localndata;
+ 
+  bool fastflag;
+
+  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);
+  }
+
+/*  printf("prelimarg %d, %zu %d", nfreq, ndata, flag);*/
+
+  /* ALLOCATION OF VARIABLES */
+
+/*  xdata = dvector(1,ndata);*/
+/*  ydata = dvector(1,ndata);*/
+  x = dvector(1,ndata);
+  y = dvector(1,ndata);
+
+  powsd = vector(1, ndata);
+  
+
+   xdata = localxdata -1;  // -1 because dvector vs *double
+   ydata = localydata -1;
+  
+    /* MULTIPLY THE SIGNAL BY A WINDOW FUNCTION, STORE RESULT IN x AND y */
+    window(x, y, xdata, ydata, ndata);
+    
+
+    /* COMPUTE POWER SPECTRAL DENSITY USING FAST FOURIER TRANSFORM */
+    power(powsd, x, y, ndata);
+
+
+      /* CHECK IF THE FREQUENCY IS IN THE REQUIRED RANGE */
+      while((centerf = bracket_real(powsd, ndata)) < minfreq || centerf > maxfreq) {
+
+	/* IF NO, SUBSTRACT IT FROM THE SIGNAL */
+	leftf = centerf - TWOPI / ndata;
+	rightf = centerf + TWOPI / ndata;
+	
+	f =  golden(phisqr, leftf, centerf, rightf, x, y, ndata);
+	
+	amph(&A, &psi, f, x, y, ndata);
+	
+	for(j=1;j<=ndata;j++){
+	  xdata[j] -= A*cos( f*(j-1) + psi );
+	  ydata[j] -= A*sin( f*(j-1) + psi );
+	}
+
+	window(x, y, xdata, ydata, ndata);
+
+	power(powsd, x, y, ndata); 
+      }   
+
+
+    leftf = centerf - TWOPI / ndata;
+    rightf = centerf + TWOPI / ndata;
+
+    /* DETERMINE THE FIRST FREQUENCY */
+    f = golden(phisqr, leftf, centerf, rightf, x, y, ndata);
+    
+    *outfreq = f ; 
+
+    /* COMPUTE AMPLITUDE AND PHASE */
+    amph(&A, &psi, f, x, y, ndata);
+    
+    /* SUBSTRACT THE FIRST HARMONIC FROM THE SIGNAL */
+    for(j=1;j<=ndata;j++){
+      xdata[j] -= A*cos( f*(j-1) + psi );
+    }  
+
+  free_dvector(x, 1, ndata);
+  free_dvector(y, 1, ndata);
+  free_vector(powsd, 1, ndata);
+  
+
+  return 1;
+}
+
+
+