Often you need to convolve a particular vector with a lot of other vectors of equal length, then the program below can be used. It illustrates the use of Rcpp/inline and FFTW used from within R.

### The program below can be used, when a vector is convolved several times

### with vectors of equal length

### In the first step the FFT of x is stored and FFTW plans are created,

### in subsequent steps the FFT of x is reused and FFTW plans recalled.

library(Rcpp)

require(inline)

### Definition of plugin, makes linkage to FFTW possible

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 )

### Convolution function used after initiation

convFFTW <- cxxfunction(

signature(x_fftIn = "numeric", yIn = "numeric"),

body = '

Rcpp::NumericMatrix x_fft(x_fftIn);

Rcpp::NumericVector y(yIn);

int ny=y.size();

int n=x_fft.ncol();

Rcpp:NumericVector retPar(n);

double in_1[n];

double out_2[n];

FILE *inputfile;

fftw_complex *in_2,*out_1,*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;

inputfile = fopen("stored.wisdom", "r");

fftw_import_wisdom_from_file(inputfile);

fclose(inputfile);

forward=fftw_plan_dft_r2c_1d(n, in_1, out_1, FFTW_MEASURE);

backward=fftw_plan_dft_c2r_1d(n, in_2, out_2, FFTW_MEASURE);

for(int i=0;i
for(int i=ny;i
fftw_execute(forward);

for(int i=0;i
in_2[i][0]=x_fft(0,i)*out_1[i][0]+x_fft(1,i)*out_1[i][1];

in_2[i][1]=x_fft(1,i)*out_1[i][0]-x_fft(0,i)*out_1[i][1];

}

fftw_execute(backward);

for(int i=0;i

fftw_destroy_plan(forward);

fftw_destroy_plan(backward);

return retPar;

', plugin="FFTWconv")

### Convolution function used for initiation / creation of FFTW plans

convFFTWstart <- 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:NumericMatrix retPar(3,n);

double in_1[n];

double out_2[n];

FILE *inputfile;

fftw_complex *in_2,*out_1,*x_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_MEASURE);

backward=fftw_plan_dft_c2r_1d(n, in_2, out_2, FFTW_MEASURE);

inputfile = fopen("stored.wisdom", "w");

fftw_export_wisdom_to_file(inputfile);

fclose(inputfile);

for(int i=0;i<(n-nx);i++) in_1[i]=0.0;

for(int i=(n-nx);i
fftw_execute(forward);

x_fft=(fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n);

for(int i=0;i
for(int j=0;j<2;j++) x_fft[i][j]=out_1[i][j];

for(int i=0;i
for(int i=ny;i
fftw_execute(forward);

for(int i=0;i
in_2[i][0]=x_fft[i][0]*out_1[i][0]+x_fft[i][1]*out_1[i][1];

in_2[i][1]=x_fft[i][1]*out_1[i][0]-x_fft[i][0]*out_1[i][1];

}

fftw_execute(backward);

for(int i=0;i
retPar(0,i)=out_2[i]/((double)n);

retPar(1,i)=x_fft[i][0];

retPar(2,i)=x_fft[i][1];

}

fftw_destroy_plan(forward);

fftw_destroy_plan(backward);

return retPar;

', plugin="FFTWconv")

convolve_fftw=function(x,y){

return(convFFTW(x,rev(y)))

}

###R function which handles Rcpp/inline functions above

conv.iteration <<- 0

convolve_fftw=function(x,y){

if(conv.iteration==0){

conv.start <<- convFFTWstart(x,rev(y))

conv.iteration <<- 1

return(conv.start[1,])

}

else return(convFFTW(conv.start[2:3,],rev(y)))

}

### The program below can be used, when a vector is convolved several times

### with vectors of equal length

### In the first step the FFT of x is stored and FFTW plans are created,

### in subsequent steps the FFT of x is reused and FFTW plans recalled.

library(Rcpp)

require(inline)

### Definition of plugin, makes linkage to FFTW possible

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 )

### Convolution function used after initiation

convFFTW <- cxxfunction(

signature(x_fftIn = "numeric", yIn = "numeric"),

body = '

Rcpp::NumericMatrix x_fft(x_fftIn);

Rcpp::NumericVector y(yIn);

int ny=y.size();

int n=x_fft.ncol();

Rcpp:NumericVector retPar(n);

double in_1[n];

double out_2[n];

FILE *inputfile;

fftw_complex *in_2,*out_1,*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;

inputfile = fopen("stored.wisdom", "r");

fftw_import_wisdom_from_file(inputfile);

fclose(inputfile);

forward=fftw_plan_dft_r2c_1d(n, in_1, out_1, FFTW_MEASURE);

backward=fftw_plan_dft_c2r_1d(n, in_2, out_2, FFTW_MEASURE);

for(int i=0;i

for(int i=0;i

in_2[i][1]=x_fft(1,i)*out_1[i][0]-x_fft(0,i)*out_1[i][1];

}

fftw_execute(backward);

for(int i=0;i

fftw_destroy_plan(forward);

fftw_destroy_plan(backward);

return retPar;

', plugin="FFTWconv")

### Convolution function used for initiation / creation of FFTW plans

convFFTWstart <- 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:NumericMatrix retPar(3,n);

double in_1[n];

double out_2[n];

FILE *inputfile;

fftw_complex *in_2,*out_1,*x_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_MEASURE);

backward=fftw_plan_dft_c2r_1d(n, in_2, out_2, FFTW_MEASURE);

inputfile = fopen("stored.wisdom", "w");

fftw_export_wisdom_to_file(inputfile);

fclose(inputfile);

for(int i=0;i<(n-nx);i++) in_1[i]=0.0;

for(int i=(n-nx);i

x_fft=(fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n);

for(int i=0;i

for(int i=0;i

for(int i=0;i

in_2[i][1]=x_fft[i][1]*out_1[i][0]-x_fft[i][0]*out_1[i][1];

}

fftw_execute(backward);

for(int i=0;i

retPar(1,i)=x_fft[i][0];

retPar(2,i)=x_fft[i][1];

}

fftw_destroy_plan(forward);

fftw_destroy_plan(backward);

return retPar;

', plugin="FFTWconv")

convolve_fftw=function(x,y){

return(convFFTW(x,rev(y)))

}

###R function which handles Rcpp/inline functions above

conv.iteration <<- 0

convolve_fftw=function(x,y){

if(conv.iteration==0){

conv.start <<- convFFTWstart(x,rev(y))

conv.iteration <<- 1

return(conv.start[1,])

}

else return(convFFTW(conv.start[2:3,],rev(y)))

}

## Comments