RcppArmadillo ifft2対Rのネイティブmvfft

Question

RcppArmadilloがFFT & 2-D FFTをサポートしています。残念ながら、ifft2（RcppArmadillo）とRのネイティブmvfft(..., inverse = TRUE)には大きな違いがあります。これは、特に私のアプリケーションでは非常に重要なゼロビンで大きくなります。違いはスカラー倍数ではありません。私はこれらの逸脱のためのドキュメンテーションやアカウントを特に第0ビンに見つけることはできません。

特に、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ノイズ外の数値偏差）でよく見られる問題ですか？
• RのネイティブFFTを呼び出すためにRcppまたはRcppArmadilloからの「低レベル」関数呼び出しがありますか？

再現性 - 以下では問題を凝縮して問題を再現しました。

#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

私は特にあなたの例のようにしないでください。 Rcppコード最小限のコード用に更新

• あなたが比較されている機能でデータを変換された静止
• あまりにも複雑であるとして - 一般的に悪い考え
• ので、私はあなたがここにあなたの入力

を修正することをお勧めすることは簡単です例。 Rでhelp(fft)は我々が簡単にRcppArmadilloを使用して再現することができ、この例

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 

とのつながり：

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の例のように、私たちの入力を回復する逆

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> 

を追加します。

私の知る限りないバグ、と私は、これは2Dの場合のために何が違うのであると信じる理由がない...

編集/フォロー：エラーがOPであり、そしてありませんアルマジロとここでの主な問題は

• は、ここでの主な問題はアルマジロのfft()は、ベクトルまたは行列上で動作するので、ないということです最小限の例

で作業していない慎重

ドキュメントを読んではありません（行列の場合）は、Rのmvfft()に対応します。アルマジロのfft2()は単に他のものであり、ここでは関係ありません。

前の例を続行/拡張しましょう。私たちは、複雑なマトリックス値を使用するように私たちのアクセサを再定義：

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

し、我々はそれを養うの寸法5×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> 

同じ結果、異なるライブラリ/パッケージ、私の知る限りがないバグのために何をするかもあります。