An example using Rcpp, inline and FFTW: Convolution

Following a discussion on the Rcpp mailing list I constructed the example below, which illustrates how to use Rcpp.plugin.maker to assure correct linkage of a library – in this case the FFTW library (Fastest Fourier Transform in the West). 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<n;i++) in_1[i]=x[i-n+nx]; fftw_execute(forward); irf_fft=(fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n); for(int i=0;i<n;i++) for(int j=0;j<2;j++) irf_fft[i][j]=out_1[i][j];

for(int i=0;i<ny;i++) in_1[i]=y[i]; for(int i=ny;i<n;i++) in_1[i]=0.0; fftw_execute(forward); y_fft=(fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n); for(int i=0;i<n;i++) for(int j=0;j<2;j++) y_fft[i][j]=out_1[i][j];

for(int i=0;i<n;i++){ in_2[i][0]=irf_fft[i][0]*y_fft[i][0]+irf_fft[i][1]*y_fft[i][1]; 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<n;i++) retPar[i]=out_2[i]/((double)n);

fftw_destroy_plan(forward); fftw_destroy_plan(backward);

return retPar; ‘, plugin=”FFTWconv”)

convolve_fftw=function(x,y){ return(convFFTW(x,rev(y))) }