【发布时间】:2019-11-02 15:27:59
【问题描述】:
我正在尝试为FINUFFT 例程编写一个 R 包装器,用于计算不均匀采样系列的 FFT。我几乎没有使用 C/C++ 的经验,所以我正在研究一个将传统傅立叶变换与 NUFFT 进行比较的示例。示例代码如下。
// this is all you must include for the finufft lib...
#include "finufft.h"
#include <complex>
// also needed for this example...
#include <stdio.h>
#include <stdlib.h>
using namespace std;
int main(int argc, char* argv[])
/* Simple example of calling the FINUFFT library from C++, using plain
arrays of C++ complex numbers, with a math test. Barnett 3/10/17
Double-precision version (see example1d1f for single-precision)
Compile with:
g++ -fopenmp example1d1.cpp -I ../src ../lib-static/libfinufft.a -o example1d1 -lfftw3 -lfftw3_omp -lm
or if you have built a single-core version:
g++ example1d1.cpp -I ../src ../lib-static/libfinufft.a -o example1d1 -lfftw3 -lm
Usage: ./example1d1
*/
{
int M = 1e6; // number of nonuniform points
int N = 1e6; // number of modes
double acc = 1e-9; // desired accuracy
nufft_opts opts; finufft_default_opts(&opts);
complex<double> I = complex<double>(0.0,1.0); // the imaginary unit
// generate some random nonuniform points (x) and complex strengths (c):
double *x = (double *)malloc(sizeof(double)*M);
complex<double>* c = (complex<double>*)malloc(sizeof(complex<double>)*M);
for (int j=0; j<M; ++j) {
x[j] = M_PI*(2*((double)rand()/RAND_MAX)-1); // uniform random in [-pi,pi)
c[j] = 2*((double)rand()/RAND_MAX)-1 + I*(2*((double)rand()/RAND_MAX)-1);
}
// allocate output array for the Fourier modes:
complex<double>* F = (complex<double>*)malloc(sizeof(complex<double>)*N);
// call the NUFFT (with iflag=+1): note N and M are typecast to BIGINT
int ier = finufft1d1(M,x,c,+1,acc,N,F,opts);
int n = 142519; // check the answer just for this mode...
complex<double> Ftest = complex<double>(0,0);
for (int j=0; j<M; ++j)
Ftest += c[j] * exp(I*(double)n*x[j]);
int nout = n+N/2; // index in output array for freq mode n
double Fmax = 0.0; // compute inf norm of F
for (int m=0; m<N; ++m) {
double aF = abs(F[m]);
if (aF>Fmax) Fmax=aF;
}
double err = abs(F[nout] - Ftest)/Fmax;
printf("1D type-1 NUFFT done. ier=%d, err in F[%d] rel to max(F) is %.3g\n",ier,n,err);
free(x); free(c); free(F);
return ier;
}
我不需要其中的大部分内容,例如生成测试序列并与传统 FFT 进行比较。此外,我想返回转换的值,而不仅仅是指示成功的错误代码。下面是我的代码。
#include "finufft.h"
#include <complex>
#include <Rcpp.h>
#include <stdlib.h>
using namespace Rcpp;
using namespace std;
// [[Rcpp::export]]
ComplexVector finufft(int M, NumericVector x, ComplexVector c, int N) {
// From example code for finufft, sets precision and default options
double acc = 1e-9;
nufft_opts opts; finufft_default_opts(&opts);
// allocate output array for the finufft routine:
complex<double>* F = (complex<double>*)malloc(sizeof(complex<double>*)*N);
// Change vector inputs from R types to C++ types
double* xd = as< double* >(x);
complex<double>* cd = as< complex<double>* >(c);
// call the NUFFT (with iflag=-1): note N and M are typecast to BIGINT
int ier = finufft1d1(M,xd,cd,-1,acc,N,F,opts);
ComplexVector Fd = as<ComplexVector>(*F);
return Fd;
}
当我尝试在 Rstudio 中获取此信息时,我收到错误“没有匹配函数调用 'as(std::complexFd 的行.我相信错误表明函数'as'没有定义(我知道这是错误的),或者'as'的参数不是正确的类型。 here 的示例包括使用 'as' 转换为 NumericVector 的示例,因此除非复杂值有一些复杂性,否则我不明白为什么这里会出现问题。
我知道使用两个命名空间存在潜在问题,但我认为这不是问题所在。我最好的猜测是,我尝试使用指针的方式存在问题,但我缺乏识别它的经验,我在网上找不到任何类似的例子来指导我。
【问题讨论】:
-
我检查文档的速度非常快,在这里您尝试使用最后一个 as 从 C++ 转到 R,当您尝试获取 ComplexVector 时,对吗?在这种情况下,您应该使用 wrap 。其次,您使用 malloc,因此您丢失了有关长度的信息,我发现很难相信 as 或 wrap 可以得到 F 的长度。我建议您查看此示例以获取 double 向量:stackoverflow.com/a/19866956/4658169
-
谢谢,我将最后一个更改为 wrap 并修复了该错误。我明白你所说的关于 F 的长度,但由于 finufft 例程的运行方式,我已将其设置为这种方式。我认为我无法在不破坏该函数调用的情况下更改携带该对象的方式。 F 中应该有 N 个条目,有没有一种好方法可以让我使用这些信息来构造一个 wrap 可以从指针 F 中处理的对象?我尝试了 Fd.attr("dim") = Dimension(1, N) 之类的方法,但在 Exporter.h 中出现错误:“无法在初始化中将 'SEXP {aka SEXPREC*}' 转换为 'double*'”
-
std::'vector 有一个名为 data() 的成员函数,它返回指向堆分配缓冲区向量使用的原始指针。当我从 c++ 调用 c 函数时,我通常使用 data()