From 0e961c4888c312e9f2ed4608cc36d3996e799b2f Mon Sep 17 00:00:00 2001
From: mcrucifix <michel.crucifix@uclouvain.be>
Date: Tue, 8 Oct 2024 20:37:18 +0200
Subject: [PATCH] golden search in C
---
NAMESPACE | 1 +
R/mfft_real.R | 74 +++++++++++++++++++++++-----------
man/mfft.Rd | 6 +--
src/fmft.c | 107 +++++++++++++++++++++++++++++++++++++++++++++++++-
4 files changed, 160 insertions(+), 28 deletions(-)
diff --git a/NAMESPACE b/NAMESPACE
index 64c9449..d9c5d53 100644
--- a/NAMESPACE
+++ b/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)
diff --git a/R/mfft_real.R b/R/mfft_real.R
index aaeac6a..6353354 100644
--- a/R/mfft_real.R
+++ b/R/mfft_real.R
@@ -1,7 +1,9 @@
-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)])
diff --git a/man/mfft.Rd b/man/mfft.Rd
index b50a031..7fc6a0c 100644
--- a/man/mfft.Rd
+++ b/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".
diff --git a/src/fmft.c b/src/fmft.c
index 84bef58..0b773e2 100644
--- a/src/fmft.c
+++ b/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;
+}
+
+
+
--
GitLab