c++ - RcppArmadillo ifft2 与 R 的原生 mvfft

Question

我注意到 RcppArmadillo 支持 FFT 和 2-D FFT。不幸的是，我的数据在 ifft2 (RcppArmadillo) 和 R 的原生 mvfft(..., inverse = TRUE) 之间存在显着差异。这在第 0 个 bin 中特别大(这在我的应用程序中非常重要)。差异不是标量倍数。我找不到任何文档或解释这些偏差，尤其是在第 0 个 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

Answer 1

我不是特别喜欢你的例子:

• 因为它还是太复杂了
• 您正在比较的函数中转换数据——通常是个坏主意
• 所以我建议你修正你的输入

这是一个更简单的例子。 help(fft) in R leads of 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

我们可以使用 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> 

并添加逆

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> 

恢复我们在 R 示例中的输入。

据我所知没有错误，而且我没有理由相信这对于 2d 情况有任何不同......

编辑/跟进:错误出在 OP，而不是 Armadillo。这里的主要问题是

• 没有仔细阅读文档
• 不使用最少的示例

这里的主要问题是 Armadillo 的 fft() 可以在 vector 或矩阵上工作，因此(在矩阵情况下)对应于 R 的 mvfft( )。 Armadillo 的 fft2() 只是其他东西，与此处无关。

让我们继续/扩展我们之前的例子。我们重新定义访问器以使用复杂的矩阵值:

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

然后定义一个 5 x 2 的复杂数组，我们将其提供给它:

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> 

这与我们在每一列上单独运行函数得到的输出相同。这也是 R 为 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> 

相同的结果，不同的库/包，据我所知没有错误。