caseyk caseyk - 3 months ago 10
C++ Question

RcppArmadillo ifft2 Vs R's native mvfft

I have noticed that

RcppArmadillo
supports FFT & 2-D FFT. Unfortunately there is a significant difference between
ifft2
(
RcppArmadillo
) and R's native
mvfft(..., inverse = TRUE)
with my data. This is especially large in the zeroth bin (which is incredibly important in my application). The difference is not a scalar multiple. I cannot find any documentation or account for these deviations especially in the zeroth bin.

I have debugged the issue specifically to the
ifft(arma::cx_mat input)
function call. Unless perhaps there is an unforeseen memory management issue, this is the culprit.

Example:
ifft2
result(1 column first 5 entries):

[1] 0.513297156-0.423498014i -0.129250939+0.300225299i
0.039722228-0.093052563i -0.007956237+0.018643534i 0.001181177-0.002768473i


mvfft
inverse result (1 column first 5 entries):

[1] 0.278131988-0.633838170i -0.195699114+0.445980950i
0.060070320-0.136894940i -0.011924932+0.027175865i 0.001754788-0.003999007i


Questions


  • Is the
    RcppArmadillo
    FFT still in development?

  • Is this a common issue across FFT variants(numerical deviations outside of FP or DP noise)?

  • Is there a 'low-level' function call from Rcpp or RcppArmadillo to call R's native FFT?



Reproducibility - Below I condensed the problem as much as I could and reproduced the issue.
Updated for minimal code
Rcpp code:

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

using namespace Rcpp;
// [[Rcpp::export]]
//profile is the dependent variable of a given variable x,
//q is a vector containing complex valued information for a single column after a tcrossprod
//Size is a scalar value which the FFT depends upon.
arma::cx_mat DebugLmnCPP( arma::cx_vec Profile, arma::cx_vec q) {
std::complex<double> oneeye (0,1);//Cmplx number (0 + 1i)
arma::cx_mat qFFT = ifft2( exp( oneeye * (Profile * q.st() ) ) );
return(qFFT );
}
// [[Rcpp::export]]
//For pedagogical purposes
arma::cx_mat DebugIFFTRCPP( arma::cx_mat input) {
arma::cx_mat qFFT = ifft2( input );
return( qFFT );
}


RCode (sorry this is sloppy)

library(Rcpp)
library(RcppArmadillo)
sourceCpp("/home/FILE.cpp")

#Use C++ function
qt <- c(6.0+0i, 5.95+0i, 0.10+0i)
prof <- 0.25* sin( (1:512)*(2*3.1415)/512 ) + 0.25#Offset Sine wave
Debug1 <- DebugLmnCPP( Profile = prof, q = qt )

#Use R function
FFTSize <- 2^9
DebugLmnR <- function(Profile, q) {
g <- (0+1i)*(as.matrix(Profile ) %*% t(q))
qFFT <- mvfft( exp(g) , inverse = TRUE) / FFTSize
return( qFFT )
}
#Call function
Debug2 <- DebugLmnR( Profile = prof, q = qt )

#Use R and C++
DebugLmnRC <- function(Profile, q) {
g <- (0+1i)*(as.matrix(Profile ) %*% t(q))
qFFT <- DebugIFFTRCPP(exp(g))
return( qFFT )
}
#Call function
Debug3 <- DebugLmnRC( Profile = prof, q = qt )
#Compare Results
Debug1[1:5,1] #CPP
Debug2[1:5,1] #R
Debug3[1:5,1] #R and CPP


yields :

> Debug1[1:5,1]
[1] 0.359632774+0.35083419i -0.037254305-0.36995074i 0.015576046+0.15288379i -0.004552119-0.03992962i
[5] 0.000967252+0.00765564i
> Debug2[1:5,1]
[1] 0.03620451+0.51053116i -0.04624384-0.55604273i 0.02204910+0.23101589i -0.00653108-0.06061692i
[5] 0.00140213+0.01167389i
> Debug3[1:5,1]
[1] 0.359632774+0.35083419i -0.037254305-0.36995074i 0.015576046+0.15288379i -0.004552119-0.03992962i
[5] 0.000967252+0.00765564i

Answer

I don't particularly like your example:

  • as it is still way too complex
  • you are transforming data in the functions you are comparing -- generally a bad idea
  • so I would suggest you fix your inputs

Here is a simpler example. help(fft) in R leads of with this example

fftR> x <- 1:4

fftR> fft(x)
[1] 10+0i -2+2i -2+0i -2-2i

fftR> fft(fft(x), inverse = TRUE)/length(x)
[1] 1+0i 2+0i 3+0i 4+0i

which we can easily reproduce using RcppArmadillo:

R> cppFunction("arma::cx_mat armafft(arma::vec x) { return fft(x); }", 
+              depends="RcppArmadillo")
R> armafft(1:4)
      [,1]
[1,] 10+0i
[2,] -2+2i
[3,] -2+0i
[4,] -2-2i
R> 

and adding the inverse

R> cppFunction("arma::cx_mat armaifft(arma::cx_mat x) { return ifft(x); }", 
+              depends="RcppArmadillo")
R> armaifft(armafft(1:4))
     [,1]
[1,] 1+0i
[2,] 2+0i
[3,] 3+0i
[4,] 4+0i
R> 

recovering our input as in the R example.

No bug as far as I can tell, and I have no reason to believe this is any different for the 2d case...

Edit/Followup: The error is with the OP, and not with Armadillo. The primary issues here are

  • not reading the documentation carefully
  • not working with a minimal examples

The main issue here is that Armadillo's fft() can work on vectors or matrices and does hence (in the matrix case) correspond to R's mvfft(). Armadillo's fft2() is simply something else and not relevant here.

Let us continue / extend our previous example. We redefine our accessor to use complex matrix values:

R> cppFunction("arma::cx_mat armafft(arma::cx_mat x) { return fft(x); }",
+              depends="RcppArmadillo")
R>

and then define a complex array of dimension 5 x 2 which we feed to it:

R> z <- array(1:10 + 1i, dim=c(5,2))
R> z
     [,1]  [,2]
[1,] 1+1i  6+1i
[2,] 2+1i  7+1i
[3,] 3+1i  8+1i
[4,] 4+1i  9+1i
[5,] 5+1i 10+1i
R> 
R> armafft(z)
              [,1]          [,2]
[1,] 15.0+5.00000i 40.0+5.00000i
[2,] -2.5+3.44095i -2.5+3.44095i
[3,] -2.5+0.81230i -2.5+0.81230i
[4,] -2.5-0.81230i -2.5-0.81230i
[5,] -2.5-3.44095i -2.5-3.44095i
R> 

This is the same output we would get from running the function separately on each column. And that is also what R does for mvfft() (cf help(fft))

R> mvfft(z)
              [,1]          [,2]
[1,] 15.0+5.00000i 40.0+5.00000i
[2,] -2.5+3.44095i -2.5+3.44095i
[3,] -2.5+0.81230i -2.5+0.81230i
[4,] -2.5-0.81230i -2.5-0.81230i
[5,] -2.5-3.44095i -2.5-3.44095i
R> 

Same result, different libraries / packages, no bug as far as I can see.

Comments