diff --git a/NAMESPACE b/NAMESPACE index d3f8f06306512f59be233b348ae1a4e40cdc8ca9..64c9449fab41f390b41ccbd8e09e9ec5a539bde9 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -13,12 +13,9 @@ export(arspec) export(cwt_morlet) export(hilbert_extension) export(mem) -export(mfft) export(mfft_anova) +export(mfft_complex) export(mfft_real) -export(mfft_real_C) -export(mfft_real_quatro) -export(mfft_real_ter) export(periodogram) export(powerspectrum.wavelet) export(reconstruct_mfft) @@ -32,8 +29,6 @@ importFrom(graphics,lines) importFrom(graphics,mtext) importFrom(graphics,par) importFrom(graphics,text) -importFrom(gsignal,hann) -importFrom(pracma,gramSchmidt) importFrom(sfsmisc,axTexpr) importFrom(sfsmisc,eaxis) importFrom(stats,ar) diff --git a/R/mfft.R b/R/mfft_complex.R similarity index 52% rename from R/mfft.R rename to R/mfft_complex.R index a17726fd11f72c9682d2c016db126e1b3ef7a83c..9d9d5cc8e5b57bc3f96afb67204027844f51eb03 100644 --- a/R/mfft.R +++ b/R/mfft_complex.R @@ -1,9 +1,43 @@ -# 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 +#' 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. The algorithm was generalised +#' by M. Crucifix when the input signal is real. +#' `mfft` is a wrapper thar redirects the data +#' to `mfft_real` or `mfft_complex` depending on the nature of +#' the input data. +#' TODO: this needs more work because the two routines +#' do not yet take exactly the same inputs. In particular, the meaning +#' the `correction` parameter is not the same +#' TODO: add Laskar reference and say that the correction=0 version +#' corresponds to Laskar's routine. +#' @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. +#' The computed frequencies are in the range given by minfreq and maxfreq. +#' @param correction : +#' The parameter is passed on to `mfft_real` or `mfft_complex` depending on +#' `xdata` is real or complex, but the meaning of the parameter is not +#' 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'. +#' @param nfreq is the number of frequencies returned, must be smaller that the length of xdata. +#' @return a `mfft_deco` object, based on a data.frame with columns "Freq", "Ampl" and "Phases". +#' note that because of a language glitch (to be fixed), "Freq" actually means "Rate" +#' @author Michel Crucifix for the R code, and David Nesvorny for most of the supporting C code doing the +#' Frequency Modified Fourier transform for Complex Numbers +#' \insertRef{sidlichovsky97aa}{gtseries} +mfft <- function(xdata, nfreq=30, minfreq=NULL, maxfreq=NULL, correction=1) { + if (is.complex(xdata)) return(mfft_complex(xdata, nfreq, minfreq, maxfreq, correction)) else + return(mfft_real(xdata, nfreq, minfreq, maxfreq, correction)) +} -#' Modified Fourier transform +#' Frequency Modified Fourier transform for Complex Numbers #' #' Implementation of the the Frequency Modified Fourier Transform #' (Sidlichovsky and Nesvorny 1997, Cel. Mech. 65, 137). @@ -18,64 +52,52 @@ #' @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 +#' @param correction : +#' Modified Fourier Transform if correction = 0; +#' Frequency Modified Fourier Transform if correction = 1; +#' FMFT with additional non-linear correction if correction = 2 #' (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 +#' @return a `mfft_deco` object, based on a data.frame with columns "Freq", "Ampl" and "Phases". +#' note that because of a language glitch (to be fixed), "Freq" actually means "Rate" #' @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) { +#' @export mfft_complex +mfft_complex <- function(xdata, nfreq=30, minfreq=NULL, maxfreq=NULL, correction=1) { - xdata = stats::as.ts(xdata) - dt = deltat(xdata) - startx = stats::start(xdata)[1] - ydata = Im(xdata) - xdata = Re(xdata) - ndata = length(xdata); + xdata <- stats::as.ts(xdata) + dt <- deltat(xdata) + startx <- stats::start(xdata)[1] + ydata <- Im(xdata) + xdata <- Re(xdata) + ndata <- length(xdata); + flag <- correction + 1 # this is an ugly translation towards Sidlichovsky's convention 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 - my_maxfreq = pi} else { - my_minfreq = minfreq * dt - my_maxfreq = maxfreq * dt} + my_minfreq <- -pi + my_maxfreq <- pi} else { + my_minfreq <- minfreq * dt + my_maxfreq <- maxfreq * dt} - signal1 = as.double(rep(0,3*nfreq)) - signal2 = as.double(rep(0,3*nfreq)) - signal3 = as.double(rep(0,3*nfreq)) + signal1 <- as.double(rep(0,3*nfreq)) + signal2 <- as.double(rep(0,3*nfreq)) + signal3 <- as.double(rep(0,3*nfreq)) - Infreq=as.integer(nfreq) - Iminfreq = as.double(my_minfreq) - Imaxfreq = as.double(my_maxfreq) - Ixdata = as.double(xdata) - Iydata = as.double(ydata) + 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), + OUT <- .C("fmft", Infreq, Iminfreq, Imaxfreq, as.integer(flag), as.integer(ndata), Ixdata, Iydata, signal1,signal2,signal3, DUP=TRUE) OUTVEC <- t(matrix(OUT[[7+flag]], 3, nfreq)) @@ -84,17 +106,21 @@ mfft <- function(xdata, minfreq=NULL, maxfreq=NULL, flag=1, nfreq=30) { 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 - OUT = data.frame(Freq=Freq, Amp=Ampl, Phases=Phase) - attr(OUT,"class") = "mfft_deco" - attr(OUT,"nfreq") = nfreq - attr(OUT,"xdata") = xdata + # if this is input is a real vector, there will be positive and negative frequencies + # corresponding to the same amplitude + # however, better use mfft_real in this case + OUT <- data.frame(Freq=Freq, Amp=Ampl, Phases=Phase) + + attr(OUT,"class") <- "mfft_deco" + attr(OUT,"nfreq") <- nfreq + attr(OUT,"xdata") <- xdata + return(OUT) } + # # diff --git a/R/mfft_real.R b/R/mfft_real.R index 9ca1735d491c0e9dbb4d543daa3501fbeeb2c501..aaeac6a3d2982e61e4e7116b028bafefceee82f5 100644 --- a/R/mfft_real.R +++ b/R/mfft_real.R @@ -1,69 +1,299 @@ -calc_amp <- function(x, nfreq, Fs){ - N <- length(x) - t <- seq(N) - xx = matrix(rep(1,N),N,1) - freqs = 2.*pi*seq(0,(N-1))/N + +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. + +mfft_analyse <- function(xdata, nfreq, fast = TRUE, nu = NULL, minfreq=NULL, maxfreq=NULL){ + + # nu can be provided, in which case the frequencies + # are considered to be given + if (!is.null(nu)) + { + # we will need two versions of nfreq, one for the cosine and one for the cosine + # if if there is a zero frequency we only need it once + nutemp <- unlist ( lapply( nu, function(a) if (a==0) return (a) else return ( c(a,-a))) ) + phase <- unlist ( lapply( nu, function(a) if (a==0) return (0) else return ( c(0, pi/2))) ) + nu <- nutemp + if (length(nu) < 2*nfreq) { + # this will typically occur if one zero frequency was provided + nu[2*nfreq] <- NA + phase[2*nfreq] <- NA + } + } else + { + nu <- rep(NA,2*nfreq) + phase <- rep(NA,2*nfreq) + } + + N <- length(xdata) + T <- N / 2 + hann <- function(N) {return (1-cos(2*pi*seq(0,(N-1))/(N)))}; + hN <- hann(N) + t <- seq(N)-1 power <- function(x) (Mod(fft(x))^2)[1:floor((N-1)/2)] - residu <- x - mean(x) - for (j in seq(nfreq)) { + + hprod <- function(x,y) sum(x * hN * y)/N + + # B_i = sum[m<= i](A_im) f_m the orthogonal base signals + # x_m what remains of the signal at time step m + + S <- rep(0,2*nfreq) + Prod <- rep(0,2*nfreq) + amp <- rep(NA,2*nfreq) + A <- matrix(0,2*nfreq,2*nfreq) + Q <- matrix(0,2*nfreq,2*nfreq) + f <- list() + x <- list() + B <- list() + # S[m] = hprod (f[m], B[m]) + freqs = 2.*pi*seq(0,(N-1))/N + x[[1]] = xdata + + + for (m in seq(2*nfreq)){ + hx <- hN*x[[m]] + # there is a little tweak here: + # if we have reached 2*nfreq and not nu[m] is na, then + # it means that we have encountered a constant somewhere, and that last step should be skipped. + if ( ! ( m == 2*nfreq && is.na(nu[m]))) { + # ok we are on business + if (is.na(nu[m])){ + # is the "m" frequency already provided ? + # either there were provided from the start and developped in the first lines of this routine + # 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); + # 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) + + # 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) + + # fmax = cmna::goldsectmax(function(f) { + # ft <- f*t + # Mod(sum(hx * exp(1i*ft)))}, + # brackets[1], brackets[2], tol=1.e-10, m=9999) + + + + 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 + phase[m] <- 0. + nu[m] <- fmax + phase[m+1] <- pi/2. + nu[m+1] <- -fmax + } else { + # again we are still in thes case where a freqency was not a priori provided + # but the more particular case where the frequency is (with tolerance) considered as zero + # f[[m]] will effectively be constant. + phase[m] <- 0. + nu[m] <- 0. + }} + + f[[m]] <- cos(nu[m]*t + phase[m]) + + if (fast) # we use analytical forms of the products + { + + for (i in seq(m)) + { + + num <- (nu[m] - nu[i])*T + nup <- (nu[m] + nu[i])*T + phim <- (phase[m] - phase[i]) + phip <- (phase[m] + phase[i]) + + Qm <- ifelse(num == 0, 1 , Q(num)) + Qp <- ifelse(nup == 0, 1 , Q(nup)) + cosm <- cos(phim) * cos(num) * Qm + sinm <- sin(phim) * sin(num) * Qm + cosp <- cos(phip) * cos(nup) * Qp + sinp <- sin(phip) * sin(nup) * Qp + + Q[m,i] = 0.5 * (cosm + cosp - sinm - sinp) + } + } else { + for (i in seq(m)) + {Q[m,i] = hprod(f[[m]],f[[i]]) # so remember the convetion + # these are symmetric matrices + # and we only populate the j <= m + } + } + + # + # before normalisation + # B[[m]] = sum_j=1,m A[m,j]*f[j] et + # B[[m]] = ( f[m] - sum_(1,m-1) (f[m] %*% B[i]) * B[i] + # B[[m]] = ( f[m] - sum_(1,m-1) + # (f[m] * sum(j=1,i A_j,i f[i]) ) * + # ( sum (A_j=1, i) A_j,i f[i] ) + + #/ ||B[[m]]|| + + # one can then verify that: + + # B[[m]] %*% B[j]] = + # ( f[m] - sum_(1,m-1) (f[m] %*% B[i]) * B[i] ) %*% B[j] = + # ( f[m] - (f[m] %*% B[j]) * B[j] ) %*% B[j] = + # ( f[m] %*% B[j] - f[m] %*% B[j]) = 0 + + # before normalisation + A[m,] = 0 + A[m, m] = 1. + if (m>1){ + # fmbi = f[m] %*% B[i] + fmbi <- rep(0,(m-1)) + for (j in seq(m-1)) for (i in seq(j)) fmbi[j] = fmbi[j] + A[j,i]*Q[m,i] + for (j in seq(m-1)) for (i in seq(j,(m-1))) A[m,j] = A[m,j] - fmbi[i]*A[i,j] + } + # so, at this stage, we have B[[m]] = sum [A[m,j]] * f[j] + + # B[[2]] = (A[2,1]*f[1] + A[2,2]*f[2]) + # B[[3]] = (A[3,1]*f[1] + A[3,2]*f[2] + A[3,3] * f[3] ) + + + norm = 0 + if (m > 1){ + for (i in seq(m)) { + norm = norm + (A[m,i]*A[m,i]*Q[i,i]) + if (i>1) for (j in seq(i-1)) { norm = norm + 2*A[m,i]*A[m,j]*Q[i,j] + }} + } else { + norm <- A[m,m]*A[m,m]*Q[m,m] + } -# temps de calcul peut etre gagne en C - - - if (Fs[j] == (-1) ){ - hresidu = gsignal::hann(N)*residu - fbase <- freqs[which.max(power(hresidu))] - brackets <- c(fbase-pi/N, fbase+pi/N); - thresidu <- t(hresidu) - fmax = cmna::goldsectmax(function(f) - (thresidu %*% cos(f*t))^2 + (thresidu %*% sin(f*t))^2, - brackets[1], brackets[2], tol=1.e-10, m=9999) - Fs[j] = fmax - } - xx <- cbind(xx, cos(Fs[j]*t), sin(Fs[j]*t)) + + A[m,] = A[m,] / sqrt(norm) + -# cette partie doit pouvoir etre acceleree en tirant partit du fait -# que l'expression analytitque de xx%*xx est connue -# voir non-working code fmft_real.C -# Il faut utiliser l'integrale connue int_0^2pi cos(x)hann(x) -# on utilise: cos(a)cos(b) = 0.5*(cos(a+b) - cos(a-b)) etc pour toutes les comb. cossin -# et int_0^2T cos(wx) = 1/W(sin[2T]-sin[0]) , ou si on applique le produit scalaire avec fenetre Hann -# et int_0^2T cos(wx)Hann[x] = (pi^2)(pi^2-w^2T^2)(sin(wT)/wT) * (coswT) -# et int_0^2T sin(wx)Hann[x] = (pi^2)(pi^2-w^2T^2)(sin(wT)/wT) * (sin(wT)) - # sin(wT)/(wT) = 1 - # et on se rappelle que la decomposition gramschmidt - # c'est simplement, avec x_i le ie vecteur - # xn_1 = x_1 / ||x_1|| - # xn_(i+1) = x_(i+1) - sum_(j=1,i) x_(i+1)%*%x(i)/||x_(i+1)|| -# en C, cela doit reduire le temps de calcul considerablement - - - B <- pracma::gramSchmidt(xx)$Q - Coefs = t(xx) %*% B - aa <- as.numeric(residu) %*% B - - signal <- B %*% t(aa) - residu <- residu - signal + Prod[m] = hprod(x[[1]],f[[m]]) + # les prods sont corrects + + 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]] + + + #print ('S[m] computed two different ways') + #print (m) + #print (S[m]) + #print (hprod(x[[1]], B[[m]])) + #print ('end test') + + # for (j in seq(m)) print (sprintf("B%i%i : %3g", i, j, hprod(B[[j]], B[[m]]))) + # if you are curious, the amplitude of the signal will be + # amp[j] = contribution of the sum( S[m]*B[m]) to the f[[m]] + # amp[j] = sum(m in seq(j)) S(m) * A[m,j] + + # for (i in seq(m)) amp[j] = amp[j]+S[m]*A[m,j] + # and the residual signal + + # x[[m+1]] = x[[m]] - S[m] * B[[m]] + + + 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]] +# print('residu final') +# print(m) +# print (x[[m+1]][seq(20)]) + + } + } - aa <- as.numeric(x) %*% B - sol = solve(t(Coefs), t(aa)) - sol1 <- matrix(sol[-1],2, nfreq) + #xtest <- 0 + #xtest2 <- 0 + #coeftests <- rep(0,m) + #phi1 <- 0.4 ; phi2 <- 0.1234 + # coefreals <- c(cos(phi1), -sin(phi1), cos(phi2), -sin(phi2)) + # for (j in seq(m)) xtest <- xtest + coefreals[j] * f[[j]] + # for (j in seq(m)) for (i in seq(j)) xtest2 <- xtest2 + S[j] * A[j,i] * f[[i]] + # for (j in seq(m)) for (i in seq(j)) coeftests[i] <- coeftests[i] + S[j] * A[j,i] + # print ((xtest - x[[1]])[seq(20)]) + # print ((xtest2 - x[[1]])[seq(20)]) + # print ("cofreals") + # print (coefreals) + # print ("coftests") + # print (coeftests) + + # donc: j'ai demontre que les deux xtests sont les memes + # coefreals est la bonne solution + # et je sais ussi que B = sum A_mj f_j + # donc je devroais simplement avour les Sj A_mj - Freqs <- c(0,Fs) - amps <- c(sol[1], sqrt ( apply(sol1^2, 2, sum) )) - Phases <- -c(0, ( apply(sol1, 2, function(i) Arg(i[1] + 1i*i[2]) ))) - - - OUT = data.frame(Freq=Freqs, Amp=amps, Phases=Phases) - attr(OUT,"class") = "mfft_deco" - return(OUT) -} + # at htis point I would like to make a littel check + # make A full + # is A Q t(A) = BtB orthonormal ? + # for (m in seq(2,2*nfreq)) for (j in seq(m-1)) A[j,m]=A[m,j] + # for (m in seq(2,2*nfreq)) for (j in seq(m-1)) Q[j,m]=Q[m,j] + # print (A*t(Q)*t(A)) + # print ('END TEST') + + #print (A) + #print (A %*% S) + # amplitudes + + mmax <- 2*nfreq + if (is.na(nu[mmax])) mmax <- mmax - 1 + + amp[1:mmax] = 0; + for (m in seq(mmax)) for (j in seq(m)) amp[j] <- amp[j] + S[m] * A[m,j] + + amp2m <- amp; + + for (m in seq(mmax)){ + # if the perivous "m" was already of the same frequency, we can merge amplitudes and set phase + if ((m > 1) && (nu[m-1] == -nu[m])){ + + phase[m] <- Arg (-1i * amp[m] + amp[m-1]) + amp[m] <- sqrt(amp[m-1]^2 + amp[m]^2) + amp[m-1] <- NA + phase[m-1] <- NA + } + } + + # at this point we should have exactly nfreq non na values + # we print a message if this is not right, and in that case I suspect some debugging will be needed. + valid_frequencies <- which(!is.na(amp)) + nu <- -nu[valid_frequencies] + amp <- amp[valid_frequencies] + phase <- phase[valid_frequencies] + if (length(valid_frequencies) != nfreq) message (sprintf("something goes wrong : %i valid frequencies, and nfreq = %i", valid_frequencies, nfreq)) -#' Modified Fourier transform for real series + OUT = data.frame(nu=nu, amp=amp, phase=phase) + return(OUT) +} + + +#' Frequency Modified Fourier transform for real series #' #' R-coded version of the Modified Fourier Transform #' with frequency correction, adapted to R. @@ -74,170 +304,103 @@ calc_amp <- function(x, nfreq, Fs){ #' A C-version should be supplied one day. #' #' @importFrom cmna goldsectmax -#' @importFrom pracma gramSchmidt -#' @importFrom gsignal hann #' @param xdata The data provided either as a time series (advised), or as a vector. +#' @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 correction: 0: no frequency correction (equivalent to Laskar); 1 : frequency correction using linear approximation ; 2: frequency correction using sythetic data (both are documented in the Sidlichovsky and Nesvorny paper. #' @param nfreq is the number of frequencies returned, must be smaller that the length of xdata. +#' @param fast (default = TRUE) uses analytical formulations for the crossproducts involving sines and cosines. +#' note: this is not really faster because the bottleneck is actually the goden section search. But more elegant. +#' @return a `mfft_deco` object, based on a data.frame with columns "Freq", "Ampl" and "Phases". #' @author Michel Crucifix #' @references #' \insertRef{sidlichovsky97aa}{gtseries} #' @examples -#' -#' set.seed(12413) -#' t = seq(1024) -#' x_orig = cos(t*0.13423167+0.00) + 1.3 * cos(t*0.119432+2.314) + 0.134994 + 0.4*cos(t*0.653167) + 0.11 * cos(t*0.78913498) + rnorm(1024)*0.12 -#' OUT <- mfft_real(x_orig) -#' print(OUT) +#' +#' data(harmonic_sample) +#' spectrum <- mfft_real(harmonic_sample$data) +#' print(spectrum) #' #' @export mfft_real -mfft_real <- function(xdata, nfreq=5, correction=TRUE){ +mfft_real <- function(xdata, nfreq=5, minfreq=NULL, maxfreq=NULL, correction = 1 , fast=TRUE){ + + N <- length(xdata) + N2 <- N/2. xdata = stats::as.ts(xdata) dt = deltat(xdata) - startx = stats::start(xdata)[1] - N = length(xdata) - times <- (seq(N)-1)*dt+startx - residu = xdata - - # will withold the definitive frequencies - Fs <- rep(-1, nfreq) - k=1 - OUT <- calc_amp(xdata, nfreq, Fs) - - Fs <- OUT$Freq[-1] - - -## correction METHOD ANALYTICAL. CURRENTLY NOT USED BUT SHOULD STILL BE EXPLORED -## -------------------------------------------------------------------------------- -## Freqs <- OUT$Freq -## amps <- OUT$Amp -## Phases <- OUT$Phase -## if (k < (nfreq)){ -## tau = (N-1) / 2. -## Qp <- function(y) { -## (1 / y) * (pi^2/ (pi^2-y^2))*(cos(y) + sin(y) / y * (3. * y^2 - pi^2)/(pi^2 - y^2))} -## Q0 <- {2. / (pi^2) -1./3.} -## -## epsilon = rep(0,nfreq) -## -## # tau = (N-1) / 2. -## -## for (j in seq(k,nfreq-1)){ -## jj = j+1 -## ysj =(Freqs[(jj+1):(nfreq+1)]-Freqs[jj])*tau -## epsilon[j] = sum(amps[(jj + 1):(nfreq+1)]*Qp(ysj)/ -## (amps[jj] * Q0 * tau) * cos(ysj + Phases[(jj+1):(nfreq+1)]-Phases[jj])) -## } -## Epsilon=c(0,epsilon) -## Freqs <- Freqs - Epsilon -## Fs <- Freqs[2:(nfreq+1)] -## } -## - - - # correction (methode 2) - -# } - if (correction){ - xdata_synthetic <- rep(0,N) - Fs <- rep(-1, nfreq) - for (i in seq(nfreq+1)) xdata_synthetic = xdata_synthetic + OUT$Amp[i]*cos(OUT$Freq[i]*seq(N) + OUT$Phase[i]) - OUT2 <- calc_amp(xdata_synthetic, nfreq, Fs) - -# print(OUT2) -# print(OUT$Freq - OUT2$Freq) - # adjust to units - OUT$Freq = OUT$Freq + (OUT$Freq - OUT2$Freq) - OUT$Amp = OUT$Amp + (OUT$Amp - OUT2$Amp) - OUT$Phase = OUT$Phase + (OUT$Phase - OUT2$Phase) - } - OUT$Freq <- OUT$Freq/dt - OUT$Phase <- OUT$Phase - startx*OUT$Freq - - O <- order(OUT$Amp, decreasing=TRUE) - - OUT$Amp = OUT$Amp[O] - OUT$Freq = OUT$Freq[O] - OUT$Phase = OUT$Phase[O] - - - attr(OUT,"class") = "mfft_deco" - attr(OUT,"nfreq") = nfreq - attr(OUT,"xdata") = xdata + my_minfreq <- ifelse (is.null(minfreq), 0, minfreq*dt) + my_maxfreq <- ifelse (is.null(maxfreq), pi, maxfreq*dt) - total_ssq <- sum(xdata^2) + startx = stats::start(xdata)[1] + N <- length(xdata) + OUT <- mfft_analyse(xdata, nfreq, fast, NULL, my_minfreq, my_maxfreq) + + + # correction (methode 2) + if (correction == 2){ + xdata_synthetic <- rep(0,N) + t <- seq(N)-1 + for (i in seq(nfreq)) xdata_synthetic = xdata_synthetic + OUT$amp[i]*cos(OUT$nu[i]*t + OUT$phase[i]) + OUT2 <- mfft_analyse(xdata_synthetic, nfreq, fast, NULL, my_minfreq, my_maxfreq) + OUT$nu = OUT$nu + (OUT$nu - OUT2$nu) + OUT$amp = OUT$amp + (OUT$amp - OUT2$amp) + OUT$phase = OUT$phase + (OUT$phase - OUT2$phase) + } else if (correction == 1){ + + for (j in seq(nfreq)){ + epsilon = OUT$amp[j] * Qprime(-2 * OUT$nu[j] * N2)*cos(2 * OUT$nu[j] * N2 + 2 * OUT$phase[j]) + if ((j+1) <= nfreq) { + for (s in seq(j+1, nfreq)) { + epsilon = epsilon + OUT$amp[s] * + ( + Qprime( (OUT$nu[s] - OUT$nu[j])*N2)*cos((OUT$nu[j] - OUT$nu[s])*N2 + OUT$phase[j] - OUT$phase[s] ) - + Qprime(( OUT$nu[s] + OUT$nu[j])*N2)*cos((OUT$nu[j] + OUT$nu[s])*N2 + OUT$phase[j] + OUT$phase[s] ) + ) + } + } + epsilon = epsilon / Qsecond0 / N2 / OUT$amp[j] + + OUT$nu[j] = OUT$nu[j] - epsilon + } + OUT <- mfft_analyse(xdata, nfreq, fast, nu = OUT$nu, my_minfreq, my_maxfreq) + } - return(OUT) -} + # account for tseries scaling (i.e. uses the vaue of deltat and start encoded + # in the time series object + OUT$nu <- OUT$nu / dt + OUT$phase <- OUT$phase - startx*OUT$nu -#' MFFT reconstruction -#' @rdname mfft_deco -#' @export reconstruct_mfft -reconstruct_mfft <- function(M){ - if (!(attr(M,"class")=="mfft_deco")) stop ("object is not a MFFT decomposition") - xdata <- attr(M,"xdata") - nfreq <- attr(M,"nfreq") - times <- seq(length(xdata))*dt(xdata) + startx(xdata) - reconstructed <- lapply(seq(nfreq), function(i) M$Amp[i]*cos(M$Freq[i]*times + M$Phase[i]) ) -} + # rearrange terms to avoid negative amplitudes and have phases in the right quandrant + # these are cosines so even functions -#' MFFT ANOVA -#' @rdname mfft_deco -#' @export mfft_anova -mfft_anova <- function(M){ - if (!(attr(M,"class")=="mfft_deco")) stop ("object is not a MFFT decomposition") - xdata <- attr(M,"xdata") - nfreq <- attr(M,"nfreq") - N <- length(xdata) - times <- seq(length(xdata))*deltat(xdata) + start(xdata) - reconstructed <- sapply(seq(nfreq), function(i) M$Amp[i]*cos(M$Freq[i]*times + M$Phase[i]) ) - cum_reconstruct <- apply(reconstructed, 1, cumsum) - residual_vars <- apply(apply(cum_reconstruct, 1, function(i) xdata-i) , 2, function(j) sum(j^2)) - var0 <- sum(xdata^2) - master_vars <- c(var0, residual_vars[-length(residual_vars)]) - p2s <- 2*seq(nfreq) - p1s <- c(0, p2s[-length(p2s)]) - F <- (master_vars - residual_vars)/residual_vars * (N - p2s)/(p2s-p1s) -} + to_be_corrected <- which (OUT$amp < 0) + if (length(to_be_corrected)){ + OUT$amp[to_be_corrected] <- - OUT$amp[to_be_corrected] + OUT$phase[to_be_corrected] <- OUT$phase[to_be_corrected] + pi + } + to_be_corrected <- which (OUT$nu < 0) + if (length(to_be_corrected)){ + OUT$nu[to_be_corrected] <- - OUT$nu[to_be_corrected] + OUT$phase[to_be_corrected] <- - OUT$phase[to_be_corrected] + } -#' @rdname mfft_deco -#' @export -as.data.frame.mfft_deco <- function(x) {data.frame(Freq=x$Freq, Amp=x$Amp, Phases=x$Phases)} - - -#' @rdname mfft_deco -#' @export -plot.mfft_deco <- function (M,periods=FALSE,...){ -# O <- order(M$Freq) - plot(abs(M$Freq), abs(M$Amp),'h',ylab="Amplitudes", xlab="", ...) - if (periods) { - frequencies <- pretty(range(M$Freq/(2*pi))) - labels <- as.character(1/frequencies) - if (0 %in% frequencies) labels[which(frequencies == 0)] = "∞" - axis(1, line=3, at=2*pi*frequencies, labels=labels) - mtext("Rate", 1, 2) - mtext("Period", 1, 4) - } else { - mtext("Rate", 1, 3) - } - points(abs(M$Freq), abs(M$Amp),'p',...) -} + # finally, order termes by decreasing amplitude + o <- order(OUT$amp, decreasing = TRUE) + OUT$amp <- OUT$amp[o] + OUT$nu <- OUT$nu[o] + OUT$phase <- OUT$phase[o] -#' @rdname mfft_deco -#' @export -lines.mfft_deco <- function (M,...){ -# O <- order(M$Freq) - lines(abs(M$Freq), abs(M$Amp),'h',...) - points(abs(M$Freq), abs(M$Amp),'p',...) -} + OUT$phase <- (OUT$phase + (2*pi)) %% (2*pi) -#' @rdname mfft_deco -#' @export -print.mfft_deco <- function (M,...){ - print.data.frame(cbind(as.data.frame(M), Period=2*pi/M$Freq)) + # rename for class compatibility + names(OUT) <- c("Freq","Amp","Phases") + attr(OUT, "class") <- "mfft_deco" + attr(OUT, "data") <- xdata + attr(OUT, "nfreq") <- nfreq + return(OUT) } - - diff --git a/R/mfft_real_quatro.R b/R/mfft_real_quatro.R deleted file mode 100644 index bee407e45cef14ff77342cc193cbbbaf3028681c..0000000000000000000000000000000000000000 --- a/R/mfft_real_quatro.R +++ /dev/null @@ -1,417 +0,0 @@ - -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. - -analyse_quatro <- function(xdata, nfreq, fast = TRUE, nu = NULL){ - - - # nu = c(0.0429, 0.0186) - # nu can be provided, in which case the frequencies - # are considered to be given - if (!is.null(nu)) - { - # we will need two versions of nfreq, one for the cosine and one for the cosine - # if if there is a zero frequency we only need it once - nutemp <- unlist ( lapply( nu, function(a) if (a==0) return (a) else return ( c(a,-a))) ) - phase <- unlist ( lapply( nu, function(a) if (a==0) return (0) else return ( c(0, pi/2))) ) - nu <- nutemp - if (length(nu) < 2*nfreq) { - # this will typically occur if one zero frequency was provided - nu[2*nfreq] <- NA - phase[2*nfreq] <- NA - } - } else - { - nu <- rep(NA,2*nfreq) - phase <- rep(NA,2*nfreq) - } - - N <- length(xdata) - T <- N / 2 - hann <- function(N) {return (1-cos(2*pi*seq(0,(N-1))/(N)))}; - hN <- hann(N) - t <- seq(N)-1 - power <- function(x) (Mod(fft(x))^2)[1:floor((N-1)/2)] - - hprod <- function(x,y) sum(x * hN * y)/N - - # B_i = sum[m<= i](A_im) f_m the orthogonal base signals - # x_m what remains of the signal at time step m - - S <- rep(0,2*nfreq) - Prod <- rep(0,2*nfreq) - amp <- rep(NA,2*nfreq) - A <- matrix(0,2*nfreq,2*nfreq) - Q <- matrix(0,2*nfreq,2*nfreq) - f <- list() - x <- list() - B <- list() - # S[m] = hprod (f[m], B[m]) - freqs = 2.*pi*seq(0,(N-1))/N - x[[1]] = xdata - - - for (m in seq(2*nfreq)){ - hx <- hN*x[[m]] - # there is a little tweak here: - # if we have reached 2*nfreq and not nu[m] is na, then - # it means that we have encountered a constant somewhere, and that last step should be skipped. - if ( ! ( m == 2*nfreq && is.na(nu[m]))) { - # ok we are on business - if (is.na(nu[m])){ - # is the "m" frequency already provided ? - # either there were provided from the start and developped in the first lines of this routine - # 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); - 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) - - # fmax = cmna::goldsectmax(function(f) { - # ft <- f*t - # Mod(sum(hx * exp(1i*ft)))}, - # brackets[1], brackets[2], tol=1.e-10, m=9999) - - - - 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 - phase[m] <- 0. - nu[m] <- fmax - phase[m+1] <- pi/2. - nu[m+1] <- -fmax - } else { - # again we are still in thes case where a freqency was not a priori provided - # but the more particular case where the frequency is (with tolerance) considered as zero - # f[[m]] will effectively be constant. - phase[m] <- 0. - nu[m] <- 0. - }} - - f[[m]] <- cos(nu[m]*t + phase[m]) - - if (fast) # we use analytical forms of the products - { - - for (i in seq(m)) - { - - num <- (nu[m] - nu[i])*T - nup <- (nu[m] + nu[i])*T - phim <- (phase[m] - phase[i]) - phip <- (phase[m] + phase[i]) - - Qm <- ifelse(num == 0, 1 , Q(num)) - Qp <- ifelse(nup == 0, 1 , Q(nup)) - cosm <- cos(phim) * cos(num) * Qm - sinm <- sin(phim) * sin(num) * Qm - cosp <- cos(phip) * cos(nup) * Qp - sinp <- sin(phip) * sin(nup) * Qp - - Q[m,i] = 0.5 * (cosm + cosp - sinm - sinp) - } - } else { - for (i in seq(m)) - {Q[m,i] = hprod(f[[m]],f[[i]]) # so remember the convetion - # these are symmetric matrices - # and we only populate the j <= m - } - } - - # - # before normalisation - # B[[m]] = sum_j=1,m A[m,j]*f[j] et - # B[[m]] = ( f[m] - sum_(1,m-1) (f[m] %*% B[i]) * B[i] - # B[[m]] = ( f[m] - sum_(1,m-1) - # (f[m] * sum(j=1,i A_j,i f[i]) ) * - # ( sum (A_j=1, i) A_j,i f[i] ) - - #/ ||B[[m]]|| - - # one can then verify that: - - # B[[m]] %*% B[j]] = - # ( f[m] - sum_(1,m-1) (f[m] %*% B[i]) * B[i] ) %*% B[j] = - # ( f[m] - (f[m] %*% B[j]) * B[j] ) %*% B[j] = - # ( f[m] %*% B[j] - f[m] %*% B[j]) = 0 - - # before normalisation - A[m,] = 0 - A[m, m] = 1. - if (m>1){ - # fmbi = f[m] %*% B[i] - fmbi <- rep(0,(m-1)) - for (j in seq(m-1)) for (i in seq(j)) fmbi[j] = fmbi[j] + A[j,i]*Q[m,i] - for (j in seq(m-1)) for (i in seq(j,(m-1))) A[m,j] = A[m,j] - fmbi[i]*A[i,j] - } - # so, at this stage, we have B[[m]] = sum [A[m,j]] * f[j] - - # B[[2]] = (A[2,1]*f[1] + A[2,2]*f[2]) - # B[[3]] = (A[3,1]*f[1] + A[3,2]*f[2] + A[3,3] * f[3] ) - - - norm = 0 - if (m > 1){ - for (i in seq(m)) { - norm = norm + (A[m,i]*A[m,i]*Q[i,i]) - if (i>1) for (j in seq(i-1)) { norm = norm + 2*A[m,i]*A[m,j]*Q[i,j] - }} - } else { - norm <- A[m,m]*A[m,m]*Q[m,m] - } - - - - A[m,] = A[m,] / sqrt(norm) - - - Prod[m] = hprod(x[[1]],f[[m]]) - # les prods sont corrects - - 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]] - - - #print ('S[m] calcule de deux facons') - #print (m) - #print (S[m]) - #print (hprod(x[[1]], B[[m]])) - #print ('end test') - - # for (j in seq(m)) print (sprintf("B%i%i : %3g", i, j, hprod(B[[j]], B[[m]]))) - # if you are curious, the amplitude of the signal will be - # amp[j] = contribution of the sum( S[m]*B[m]) to the f[[m]] - # amp[j] = sum(m in seq(j)) S(m) * A[m,j] - - # for (i in seq(m)) amp[j] = amp[j]+S[m]*A[m,j] - # and the residual signal - - # x[[m+1]] = x[[m]] - S[m] * B[[m]] - - - 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]] -# print('residu final') -# print(m) -# print (x[[m+1]][seq(20)]) - - } - - } - - #print ('TESTESTEST') - #xtest <- 0 - #xtest2 <- 0 - #coeftests <- rep(0,m) - #phi1 <- 0.4 ; phi2 <- 0.1234 - # coefreals <- c(cos(phi1), -sin(phi1), cos(phi2), -sin(phi2)) - # for (j in seq(m)) xtest <- xtest + coefreals[j] * f[[j]] - # for (j in seq(m)) for (i in seq(j)) xtest2 <- xtest2 + S[j] * A[j,i] * f[[i]] - # for (j in seq(m)) for (i in seq(j)) coeftests[i] <- coeftests[i] + S[j] * A[j,i] - # print ((xtest - x[[1]])[seq(20)]) - # print ('ENDTESTESTET1') - # print ((xtest2 - x[[1]])[seq(20)]) - # print ('ENDTESTESTET2') - # print ("cofreals") - # print (coefreals) - # print ("coftests") - # print (coeftests) - - # donc: j'ai demontre que les deux xtests sont les memes - # coefreals est la bonne solution - # et je sais ussi que B = sum A_mj f_j - # donc je devroais simplement avour les Sj A_mj - - - # at htis point I would like to make a littel check - # make A full - # is A Q t(A) = BtB orthonormal ? - #print ('TEST') - #for (m in seq(2,2*nfreq)) for (j in seq(m-1)) A[j,m]=A[m,j] - #for (m in seq(2,2*nfreq)) for (j in seq(m-1)) Q[j,m]=Q[m,j] - #print (A*t(Q)*t(A)) - print ('END TEST') - - print (A) - print (A %*% S) - # amplitudes - - mmax <- 2*nfreq - if (is.na(nu[mmax])) mmax <- mmax - 1 - - amp[1:mmax] = 0; - for (m in seq(mmax)) for (j in seq(m)) amp[j] <- amp[j] + S[m] * A[m,j] -# print ('ampp...') -# print (amp) - - amp2m <- amp; - - for (m in seq(mmax)){ - # if the perivous "m" was already of the same frequency, we can merge amplitudes and set phase - if ((m > 1) && (nu[m-1] == -nu[m])){ - - phase[m] <- Arg (-1i * amp[m] + amp[m-1]) - amp[m] <- sqrt(amp[m-1]^2 + amp[m]^2) - print ("phase") - print (c(m, phase[m], amp[m], amp[m-1])) - amp[m-1] <- NA - phase[m-1] <- NA - } - } -# print ('amp2m') -# print (amp2m) -# print ('S') -# print (S) -# print ('PRod') -# print (Prod) - - # at this point we should have exactly nfreq non na values - # we print a message if this is not right, and in that case I suspect some debugging will be needed. - - valid_frequencies <- which(!is.na(amp)) - nu <- -nu[valid_frequencies] - amp <- amp[valid_frequencies] - phase <- phase[valid_frequencies] - if (length(valid_frequencies) != nfreq) message (sprintf("something goes wrong : %i valid frequencies, and nfreq = %i", valid_frequencies, nfreq)) -# print ("nu") -# print(nu) -# print ("amp") -# print (amp) - - OUT = data.frame(nu=nu, amp=amp, phase=phase) - return(OUT) -} - - -#' Modified Fourier transform for real series (variant) -#' -#' R-coded version of the Modified Fourier Transform -#' with frequency correction, adapted to R. -#' much slower than mfft (for complex numbers) as the latter is -#' mainly written in C, but is physically -#' more interpretable if signal is real, because -#' it is designed to have no imaginary part in the residual -#' A C-version should be supplied one day. -#' -#' @importFrom cmna goldsectmax -#' @param xdata The data provided either as a time series (advised), or as a vector. -#' @param nfreq is the number of frequencies returned, must be smaller that the length of xdata. -#' @param fast (default = TRUE) uses analytical formulations for the crossproducts involving sines and cosines -#' @param correction: 0: no frequency correction (equivalent to Laskar); 1 : frequency correction using linear approximation ; 2: frequency correction using sythetic data -#' @author Michel Crucifix -#' @references -#' \insertRef{sidlichovsky97aa}{gtseries} -#' @examples -#' -#' set.seed(12413) -#' t = seq(1024) -#' x_orig = cos(t*0.13423167+0.00) + 1.3 * cos(t*0.119432+2.314) + 0.134994 + 0.4*cos(t*0.653167) + 0.11 * cos(t*0.78913498) + rnorm(1024)*0.12 -#' OUT <- mfft_real(x_orig) -#' print(OUT) -#' -#' @export mfft_real_quatro - # will withold the definitive frequencies -mfft_real_quatro <- function(xdata, nfreq=5, correction=1, fast=TRUE){ - N <- length(xdata) - N2 <- N/2. - xdata = stats::as.ts(xdata) - dt = deltat(xdata) - startx = stats::start(xdata)[1] - N <- length(xdata) - OUT <- analyse_quatro(xdata, nfreq, fast) - - - # correction (methode 2) -# } - if (correction == 2){ - xdata_synthetic <- rep(0,N) - t <- seq(N)-1 - for (i in seq(nfreq)) xdata_synthetic = xdata_synthetic + OUT$amp[i]*cos(OUT$nu[i]*t + OUT$phase[i]) - OUT2 <- analyse_quatro(xdata_synthetic, nfreq, fast) - OUT$nu = OUT$nu + (OUT$nu - OUT2$nu) - OUT$amp = OUT$amp + (OUT$amp - OUT2$amp) - OUT$phase = OUT$phase + (OUT$phase - OUT2$phase) - } else if (correction == 1){ - - for (j in seq(nfreq)){ - epsilon = OUT$amp[j] * Qprime(-2 * OUT$nu[j] * N2)*cos(2 * OUT$nu[j] * N2 + 2 * OUT$phase[j]) -# print ('epsilon') -# print (c(epsilon,j, OUT$amp[j], OUT$nu[j], OUT$phase[j])) - if ((j+1) <= nfreq) { for (s in seq(j+1, nfreq)) { - epsilon = epsilon + OUT$amp[s] * - ( - Qprime( (OUT$nu[s] - OUT$nu[j])*N2)*cos((OUT$nu[j] - OUT$nu[s])*N2 + OUT$phase[j] - OUT$phase[s] ) - - Qprime(( OUT$nu[s] + OUT$nu[j])*N2)*cos((OUT$nu[j] + OUT$nu[s])*N2 + OUT$phase[j] + OUT$phase[s] ) ) - }} -# print(sprintf("j= %i, OUT$nu[j]= %.4f, OUT$amp[j] = %.4f", j, OUT$nu[j], OUT$amp[j])) - epsilon = epsilon / Qsecond0 / N2 / OUT$amp[j] - - OUT$nu[j] = OUT$nu[j] - epsilon - } - OUT <- analyse_quatro(xdata, nfreq, fast, nu = OUT$nu) - } - - # account for tseries scaling (i.e. uses the vaue of deltat and start encoded - # in the time series object - OUT$nu <- OUT$nu / dt - OUT$phase <- OUT$phase - startx*OUT$nu - - # rearrange terms to avoid negative amplitudes and have phases in the right quandrant - # these are cosines so even functions - - to_be_corrected <- which (OUT$amp < 0) - if (length(to_be_corrected)){ - OUT$amp[to_be_corrected] <- - OUT$amp[to_be_corrected] - OUT$phase[to_be_corrected] <- OUT$phase[to_be_corrected] + pi - } - - to_be_corrected <- which (OUT$nu < 0) - if (length(to_be_corrected)){ - OUT$nu[to_be_corrected] <- - OUT$nu[to_be_corrected] - OUT$phase[to_be_corrected] <- - OUT$phase[to_be_corrected] - } - - # finally, order termes by decreasing amplitude - o <- order(OUT$amp, decreasing = TRUE) - OUT$amp <- OUT$amp[o] - OUT$nu <- OUT$nu[o] - OUT$phase <- OUT$phase[o] - - OUT$phase <- (OUT$phase + (2*pi)) %% (2*pi) - - - # rename for class compatibility - names(OUT) <- c("Freq","Amp","Phases") - attr(OUT, "class") <- "mfft_deco" - attr(OUT, "data") <- xdata - attr(OUT, "nfreq") <- nfreq - return(OUT) -} - diff --git a/TODO b/TODO index d2b4176b1d5902452dff48ccbe375b885d2f35dc..10d4eb3f763a9571944375cf64612382748c7a5e 100644 --- a/TODO +++ b/TODO @@ -1,8 +1,7 @@ -- make nicer graphs for mfft <- in prgoress -- finalise computations for mfft_real_ter (use notebook) <- seems to be ok forfrequenices -- there still seems to be a bias on amplitudes that is not well corrected with the analytical methods, when there are several frequencies. As if the first components where eating a bit too much power. Not sure why, and not sure how to correct this. - But the bias is much less important with the complex version. So there must be a good mathematical reason that is not - obvious to me (in mfft_complex not quite the same power goes in omega and -omega, so maybe the explanation is there ?) -- use that one only; deprecate mfft_real +- make nicer graphs for mfft <- in progress - work on documentation +- document data and add astronomical data (documentation has to be revised for harmonic and pseudo_log) +- think more about the significance spectrum +- revis the haar bit +- revise the periodogram - and make a paper ! diff --git a/deprecated/mfft_real.R b/deprecated/mfft_real.R new file mode 100644 index 0000000000000000000000000000000000000000..9ca1735d491c0e9dbb4d543daa3501fbeeb2c501 --- /dev/null +++ b/deprecated/mfft_real.R @@ -0,0 +1,243 @@ +calc_amp <- function(x, nfreq, Fs){ + N <- length(x) + t <- seq(N) + xx = matrix(rep(1,N),N,1) + freqs = 2.*pi*seq(0,(N-1))/N + power <- function(x) (Mod(fft(x))^2)[1:floor((N-1)/2)] + residu <- x - mean(x) + for (j in seq(nfreq)) { + +# temps de calcul peut etre gagne en C + + + if (Fs[j] == (-1) ){ + hresidu = gsignal::hann(N)*residu + fbase <- freqs[which.max(power(hresidu))] + brackets <- c(fbase-pi/N, fbase+pi/N); + thresidu <- t(hresidu) + fmax = cmna::goldsectmax(function(f) + (thresidu %*% cos(f*t))^2 + (thresidu %*% sin(f*t))^2, + brackets[1], brackets[2], tol=1.e-10, m=9999) + Fs[j] = fmax + } + + xx <- cbind(xx, cos(Fs[j]*t), sin(Fs[j]*t)) + +# cette partie doit pouvoir etre acceleree en tirant partit du fait +# que l'expression analytitque de xx%*xx est connue +# voir non-working code fmft_real.C +# Il faut utiliser l'integrale connue int_0^2pi cos(x)hann(x) +# on utilise: cos(a)cos(b) = 0.5*(cos(a+b) - cos(a-b)) etc pour toutes les comb. cossin +# et int_0^2T cos(wx) = 1/W(sin[2T]-sin[0]) , ou si on applique le produit scalaire avec fenetre Hann +# et int_0^2T cos(wx)Hann[x] = (pi^2)(pi^2-w^2T^2)(sin(wT)/wT) * (coswT) +# et int_0^2T sin(wx)Hann[x] = (pi^2)(pi^2-w^2T^2)(sin(wT)/wT) * (sin(wT)) + # sin(wT)/(wT) = 1 + # et on se rappelle que la decomposition gramschmidt + # c'est simplement, avec x_i le ie vecteur + # xn_1 = x_1 / ||x_1|| + # xn_(i+1) = x_(i+1) - sum_(j=1,i) x_(i+1)%*%x(i)/||x_(i+1)|| +# en C, cela doit reduire le temps de calcul considerablement + + + B <- pracma::gramSchmidt(xx)$Q + Coefs = t(xx) %*% B + aa <- as.numeric(residu) %*% B + + signal <- B %*% t(aa) + residu <- residu - signal + } + + aa <- as.numeric(x) %*% B + sol = solve(t(Coefs), t(aa)) + sol1 <- matrix(sol[-1],2, nfreq) + + Freqs <- c(0,Fs) + amps <- c(sol[1], sqrt ( apply(sol1^2, 2, sum) )) + Phases <- -c(0, ( apply(sol1, 2, function(i) Arg(i[1] + 1i*i[2]) ))) + + + OUT = data.frame(Freq=Freqs, Amp=amps, Phases=Phases) + attr(OUT,"class") = "mfft_deco" + return(OUT) +} + + + +#' Modified Fourier transform for real series +#' +#' R-coded version of the Modified Fourier Transform +#' with frequency correction, adapted to R. +#' much slower than mfft (for complex numbers) as the latter is +#' mainly written in C, but is physically +#' more interpretable if signal is real, because +#' it is designed to have no imaginary part in the residual +#' A C-version should be supplied one day. +#' +#' @importFrom cmna goldsectmax +#' @importFrom pracma gramSchmidt +#' @importFrom gsignal hann +#' @param xdata The data provided either as a time series (advised), or as a vector. +#' @param nfreq is the number of frequencies returned, must be smaller that the length of xdata. +#' @author Michel Crucifix +#' @references +#' \insertRef{sidlichovsky97aa}{gtseries} +#' @examples +#' +#' set.seed(12413) +#' t = seq(1024) +#' x_orig = cos(t*0.13423167+0.00) + 1.3 * cos(t*0.119432+2.314) + 0.134994 + 0.4*cos(t*0.653167) + 0.11 * cos(t*0.78913498) + rnorm(1024)*0.12 +#' OUT <- mfft_real(x_orig) +#' print(OUT) +#' +#' @export mfft_real +mfft_real <- function(xdata, nfreq=5, correction=TRUE){ + xdata = stats::as.ts(xdata) + dt = deltat(xdata) + startx = stats::start(xdata)[1] + N = length(xdata) + times <- (seq(N)-1)*dt+startx + residu = xdata + + # will withold the definitive frequencies + Fs <- rep(-1, nfreq) + k=1 + OUT <- calc_amp(xdata, nfreq, Fs) + + Fs <- OUT$Freq[-1] + + +## correction METHOD ANALYTICAL. CURRENTLY NOT USED BUT SHOULD STILL BE EXPLORED +## -------------------------------------------------------------------------------- +## Freqs <- OUT$Freq +## amps <- OUT$Amp +## Phases <- OUT$Phase +## if (k < (nfreq)){ +## tau = (N-1) / 2. +## Qp <- function(y) { +## (1 / y) * (pi^2/ (pi^2-y^2))*(cos(y) + sin(y) / y * (3. * y^2 - pi^2)/(pi^2 - y^2))} +## Q0 <- {2. / (pi^2) -1./3.} +## +## epsilon = rep(0,nfreq) +## +## # tau = (N-1) / 2. +## +## for (j in seq(k,nfreq-1)){ +## jj = j+1 +## ysj =(Freqs[(jj+1):(nfreq+1)]-Freqs[jj])*tau +## epsilon[j] = sum(amps[(jj + 1):(nfreq+1)]*Qp(ysj)/ +## (amps[jj] * Q0 * tau) * cos(ysj + Phases[(jj+1):(nfreq+1)]-Phases[jj])) +## } +## Epsilon=c(0,epsilon) +## Freqs <- Freqs - Epsilon +## Fs <- Freqs[2:(nfreq+1)] +## } +## + + + # correction (methode 2) + +# } + if (correction){ + xdata_synthetic <- rep(0,N) + Fs <- rep(-1, nfreq) + for (i in seq(nfreq+1)) xdata_synthetic = xdata_synthetic + OUT$Amp[i]*cos(OUT$Freq[i]*seq(N) + OUT$Phase[i]) + OUT2 <- calc_amp(xdata_synthetic, nfreq, Fs) + +# print(OUT2) +# print(OUT$Freq - OUT2$Freq) + # adjust to units + OUT$Freq = OUT$Freq + (OUT$Freq - OUT2$Freq) + OUT$Amp = OUT$Amp + (OUT$Amp - OUT2$Amp) + OUT$Phase = OUT$Phase + (OUT$Phase - OUT2$Phase) + } + OUT$Freq <- OUT$Freq/dt + OUT$Phase <- OUT$Phase - startx*OUT$Freq + + O <- order(OUT$Amp, decreasing=TRUE) + + OUT$Amp = OUT$Amp[O] + OUT$Freq = OUT$Freq[O] + OUT$Phase = OUT$Phase[O] + + + + attr(OUT,"class") = "mfft_deco" + attr(OUT,"nfreq") = nfreq + attr(OUT,"xdata") = xdata + + total_ssq <- sum(xdata^2) + + return(OUT) +} + +#' MFFT reconstruction +#' @rdname mfft_deco +#' @export reconstruct_mfft +reconstruct_mfft <- function(M){ + if (!(attr(M,"class")=="mfft_deco")) stop ("object is not a MFFT decomposition") + xdata <- attr(M,"xdata") + nfreq <- attr(M,"nfreq") + times <- seq(length(xdata))*dt(xdata) + startx(xdata) + reconstructed <- lapply(seq(nfreq), function(i) M$Amp[i]*cos(M$Freq[i]*times + M$Phase[i]) ) +} + +#' MFFT ANOVA +#' @rdname mfft_deco +#' @export mfft_anova +mfft_anova <- function(M){ + if (!(attr(M,"class")=="mfft_deco")) stop ("object is not a MFFT decomposition") + xdata <- attr(M,"xdata") + nfreq <- attr(M,"nfreq") + N <- length(xdata) + times <- seq(length(xdata))*deltat(xdata) + start(xdata) + reconstructed <- sapply(seq(nfreq), function(i) M$Amp[i]*cos(M$Freq[i]*times + M$Phase[i]) ) + cum_reconstruct <- apply(reconstructed, 1, cumsum) + residual_vars <- apply(apply(cum_reconstruct, 1, function(i) xdata-i) , 2, function(j) sum(j^2)) + var0 <- sum(xdata^2) + master_vars <- c(var0, residual_vars[-length(residual_vars)]) + p2s <- 2*seq(nfreq) + p1s <- c(0, p2s[-length(p2s)]) + F <- (master_vars - residual_vars)/residual_vars * (N - p2s)/(p2s-p1s) +} + + +#' @rdname mfft_deco +#' @export +as.data.frame.mfft_deco <- function(x) {data.frame(Freq=x$Freq, Amp=x$Amp, Phases=x$Phases)} + + +#' @rdname mfft_deco +#' @export +plot.mfft_deco <- function (M,periods=FALSE,...){ +# O <- order(M$Freq) + plot(abs(M$Freq), abs(M$Amp),'h',ylab="Amplitudes", xlab="", ...) + if (periods) { + frequencies <- pretty(range(M$Freq/(2*pi))) + labels <- as.character(1/frequencies) + if (0 %in% frequencies) labels[which(frequencies == 0)] = "∞" + axis(1, line=3, at=2*pi*frequencies, labels=labels) + mtext("Rate", 1, 2) + mtext("Period", 1, 4) + } else { + mtext("Rate", 1, 3) + } + points(abs(M$Freq), abs(M$Amp),'p',...) +} + +#' @rdname mfft_deco +#' @export +lines.mfft_deco <- function (M,...){ +# O <- order(M$Freq) + lines(abs(M$Freq), abs(M$Amp),'h',...) + points(abs(M$Freq), abs(M$Amp),'p',...) +} + + +#' @rdname mfft_deco +#' @export +print.mfft_deco <- function (M,...){ + print.data.frame(cbind(as.data.frame(M), Period=2*pi/M$Freq)) +} + + + diff --git a/R/mfft_real_C.R b/deprecated/mfft_real_C.R similarity index 100% rename from R/mfft_real_C.R rename to deprecated/mfft_real_C.R diff --git a/R/mfft_real_ter.R b/deprecated/mfft_real_ter.R similarity index 100% rename from R/mfft_real_ter.R rename to deprecated/mfft_real_ter.R diff --git a/man/mfft.Rd b/man/mfft.Rd index f547f6e31dbc060a3851ee2daf3e6beb72d77e0d..b50a031fcc6c9afc06b1d6fe67145927992b1cbf 100644 --- a/man/mfft.Rd +++ b/man/mfft.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/mfft.R +% Please edit documentation in R/mfft_complex.R \name{mfft} \alias{mfft} \title{Modified Fourier transform} \usage{ -mfft(xdata, minfreq = NULL, maxfreq = NULL, flag = 1, nfreq = 30) +mfft(xdata, minfreq = NULL, maxfreq = NULL, correction = 1, nfreq = 30) } \arguments{ \item{xdata}{The data provided either as a time series (advised), or as a vector. @@ -12,30 +12,40 @@ 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). +with the time resolution encoded in `xdata` if the latter is a time series. The computed frequencies are in the range given by minfreq and maxfreq.} +\item{correction}{: +The parameter is passed on to `mfft_real` or `mfft_complex` depending on +`xdata` is real or complex, but the meaning of the parameter is not +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{ -STILL NEED TO BE DESCRIBED +a `mfft_deco` object, based on a data.frame with columns "Freq", "Ampl" and "Phases". + note that because of a language glitch (to be fixed), "Freq" actually means "Rate" } \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 +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. +(psi_j) in its decomposition. The algorithm was generalised +by M. Crucifix when the input signal is real. +`mfft` is a wrapper thar redirects the data +to `mfft_real` or `mfft_complex` depending on the nature of +the input data. +TODO: this needs more work because the two routines +do not yet take exactly the same inputs. In particular, the meaning +the `correction` parameter is not the same +TODO: add Laskar reference and say that the correction=0 version + corresponds to Laskar's routine. } \author{ Michel Crucifix for the R code, and David Nesvorny for most of the supporting C code doing the -actual computations + Frequency Modified Fourier transform for Complex Numbers \insertRef{sidlichovsky97aa}{gtseries} } diff --git a/man/mfft_deco.Rd b/man/mfft_deco.Rd index a42683f743a457ae2fead3d1a86ad357b7e3d2c1..0e6f5c9a8c92bfe999f502d085460edcd5a1cff2 100644 --- a/man/mfft_deco.Rd +++ b/man/mfft_deco.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/mfft_real.R +% Please edit documentation in R/mfft_support.R \name{reconstruct_mfft} \alias{reconstruct_mfft} \alias{mfft_anova} @@ -25,4 +25,5 @@ mfft_anova(M) MFFT reconstruction MFFT ANOVA +not ready. do not use. } diff --git a/man/mfft_real.Rd b/man/mfft_real.Rd index 563ad376859d76acb597aa53d0ff77e5cd790d0d..3d21180787c2e2ec4210eb562d8e260718272cd7 100644 --- a/man/mfft_real.Rd +++ b/man/mfft_real.Rd @@ -2,14 +2,33 @@ % Please edit documentation in R/mfft_real.R \name{mfft_real} \alias{mfft_real} -\title{Modified Fourier transform for real series} +\title{Frequency Modified Fourier transform for real series} \usage{ -mfft_real(xdata, nfreq = 5, correction = TRUE) +mfft_real( + xdata, + nfreq = 5, + minfreq = NULL, + maxfreq = NULL, + correction = 1, + fast = TRUE +) } \arguments{ \item{xdata}{The data provided either as a time series (advised), or as a vector.} \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.} + +\item{fast}{(default = TRUE) uses analytical formulations for the crossproducts involving sines and cosines. +note: this is not really faster because the bottleneck is actually the goden section search. But more elegant.} + +\item{correction:}{0: no frequency correction (equivalent to Laskar); 1 : frequency correction using linear approximation ; 2: frequency correction using sythetic data (both are documented in the Sidlichovsky and Nesvorny paper.} +} +\value{ +a `mfft_deco` object, based on a data.frame with columns "Freq", "Ampl" and "Phases". } \description{ R-coded version of the Modified Fourier Transform @@ -22,11 +41,9 @@ A C-version should be supplied one day. } \examples{ -set.seed(12413) -t = seq(1024) -x_orig = cos(t*0.13423167+0.00) + 1.3 * cos(t*0.119432+2.314) + 0.134994 + 0.4*cos(t*0.653167) + 0.11 * cos(t*0.78913498) + rnorm(1024)*0.12 -OUT <- mfft_real(x_orig) -print(OUT) +data(harmonic_sample) +spectrum <- mfft_real(harmonic_sample$data) +print(spectrum) } \references{ diff --git a/man/mfft_real_C.Rd b/man/mfft_real_C.Rd deleted file mode 100644 index 9a0dc5b060e61713107fc98be48913dec79f5a59..0000000000000000000000000000000000000000 --- a/man/mfft_real_C.Rd +++ /dev/null @@ -1,41 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/mfft_real_C.R -\name{mfft_real_C} -\alias{mfft_real_C} -\title{Modified Fourier transform (Real C)} -\usage{ -mfft_real_C(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/mfft_real_quatro.Rd b/man/mfft_real_quatro.Rd deleted file mode 100644 index 1f130018f7b7748e61af5cf59eb633e0bb6e3712..0000000000000000000000000000000000000000 --- a/man/mfft_real_quatro.Rd +++ /dev/null @@ -1,41 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/mfft_real_quatro.R -\name{mfft_real_quatro} -\alias{mfft_real_quatro} -\title{Modified Fourier transform for real series (variant)} -\usage{ -mfft_real_quatro(xdata, nfreq = 5, correction = 1, fast = TRUE) -} -\arguments{ -\item{xdata}{The data provided either as a time series (advised), or as a vector.} - -\item{nfreq}{is the number of frequencies returned, must be smaller that the length of xdata.} - -\item{fast}{(default = TRUE) uses analytical formulations for the crossproducts involving sines and cosines} - -\item{correction:}{0: no frequency correction (equivalent to Laskar); 1 : frequency correction using linear approximation ; 2: frequency correction using sythetic data} -} -\description{ -R-coded version of the Modified Fourier Transform -with frequency correction, adapted to R. -much slower than mfft (for complex numbers) as the latter is -mainly written in C, but is physically -more interpretable if signal is real, because -it is designed to have no imaginary part in the residual -A C-version should be supplied one day. -} -\examples{ - -set.seed(12413) -t = seq(1024) -x_orig = cos(t*0.13423167+0.00) + 1.3 * cos(t*0.119432+2.314) + 0.134994 + 0.4*cos(t*0.653167) + 0.11 * cos(t*0.78913498) + rnorm(1024)*0.12 -OUT <- mfft_real(x_orig) -print(OUT) - -} -\references{ -\insertRef{sidlichovsky97aa}{gtseries} -} -\author{ -Michel Crucifix -} diff --git a/man/mfft_real_ter.Rd b/man/mfft_real_ter.Rd deleted file mode 100644 index 23012375ed4854eb522c366cba970bf0b0ac40f0..0000000000000000000000000000000000000000 --- a/man/mfft_real_ter.Rd +++ /dev/null @@ -1,41 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/mfft_real_ter.R -\name{mfft_real_ter} -\alias{mfft_real_ter} -\title{Modified Fourier transform for real series (variant)} -\usage{ -mfft_real_ter(xdata, nfreq = 5, correction = 1, fast = TRUE) -} -\arguments{ -\item{xdata}{The data provided either as a time series (advised), or as a vector.} - -\item{nfreq}{is the number of frequencies returned, must be smaller that the length of xdata.} - -\item{fast}{(default = TRUE) uses analytical formulations for the crossproducts involving sines and cosines} - -\item{correction:}{0: no frequency correction (equivalent to Laskar); 1 : frequency correction using linear approximation ; 2: frequency correction using sythetic data} -} -\description{ -R-coded version of the Modified Fourier Transform -with frequency correction, adapted to R. -much slower than mfft (for complex numbers) as the latter is -mainly written in C, but is physically -more interpretable if signal is real, because -it is designed to have no imaginary part in the residual -A C-version should be supplied one day. -} -\examples{ - -set.seed(12413) -t = seq(1024) -x_orig = cos(t*0.13423167+0.00) + 1.3 * cos(t*0.119432+2.314) + 0.134994 + 0.4*cos(t*0.653167) + 0.11 * cos(t*0.78913498) + rnorm(1024)*0.12 -OUT <- mfft_real(x_orig) -print(OUT) - -} -\references{ -\insertRef{sidlichovsky97aa}{gtseries} -} -\author{ -Michel Crucifix -}