The resulting R function is named convolve_fftw
The example also illustrates the general use of Rcpp and inline
library(Rcpp)
require(inline)
plug <- Rcpp:::Rcpp.plugin.maker(include.before = "#include
libs = paste("-L/usr/local/lib/R/site-library/Rcpp/lib -lRcpp",
"-Wl,-rpath,/usr/local/lib/R/site-library/Rcpp/lib",
"-lfftw3 -lm -L/usr/lib"))
registerPlugin("FFTWconv", plug )
convFFTW <- cxxfunction(
signature(xIn = "numeric", yIn = "numeric"),
body = '
Rcpp::NumericVector x(xIn);
Rcpp::NumericVector y(yIn);
int nx=x.size();
int ny=y.size();
int n=nx+ny-1;
Rcpp:NumericVector retPar(n,0.0);
double in_1[n];
double out_2[n];
fftw_complex *in_2,*out_1,*irf_fft,*y_fft;
in_2=(fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n);
out_1=(fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n);
fftw_plan forward;
fftw_plan backward;
forward=fftw_plan_dft_r2c_1d(n, in_1, out_1, FFTW_ESTIMATE);
backward=fftw_plan_dft_c2r_1d(n, in_2, out_2, FFTW_ESTIMATE);
for(int i=0;i<(n-nx);i++) in_1[i]=0.0;
for(int i=(n-nx);i
irf_fft=(fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n);
for(int i=0;i
for(int i=0;i
y_fft=(fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n);
for(int i=0;i
for(int i=0;i
in_2[i][1]=irf_fft[i][1]*y_fft[i][0]-irf_fft[i][0]*y_fft[i][1];
}
fftw_execute(backward);
for(int i=0;i
fftw_destroy_plan(forward);
fftw_destroy_plan(backward);
return retPar;
', plugin="FFTWconv")
convolve_fftw=function(x,y){
return(convFFTW(x,rev(y)))
}
0 comments:
Post a Comment