RcppArmadillo ifft2 Vs R's native mvfft

我注意到RcppArmadillo支持FFT&二维FFT。不幸的是，根据我的数据，ifft2(RcppArmadillo(和R的原生mvfft(..., inverse = TRUE)之间存在显著差异。这在第零个bin中特别大(这在我的应用程序中非常重要(。差异不是标量倍数。我找不到任何文件或解释这些偏差，尤其是在第零个箱子里。

我已经专门针对ifft(arma::cx_mat input)函数调用调试了这个问题。除非出现意外的内存管理问题，否则这就是罪魁祸首。

示例：ifft2结果(1列前5个条目(：

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

mvfft逆结果(1列前5个条目(：

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

问题

• RcppArmadillo FFT还在开发中吗
• 这是FFT变体中的常见问题吗(FP或DP噪声之外的数值偏差(
• 是否有来自Rcpp或RcppArmadillo的"低级"函数调用来调用R的本地FFT

再现性-下面我尽可能地浓缩了这个问题，并再现了这个问题。更新为最小代码Rcpp代码：

#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(抱歉这太草率了(

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

收益率：

> 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