Help with convolving data sets

xraygenfit — Thu, 13 Dec 2007 18:59:21 GMT

Hi guys,

I've been banging my head against the wall trying to convolve two data sets. They are of equal length. One is a gaussian, and the other is a Lorentzian. I haven't been able to get it to work using a regular convolution.

i.e.


for (j = 0; j < 64; j += 1)
{ 
 y= 0.0; 
 for (i = 0; i <= j; i++)
 {
 y += x * h[j - i];
 }
}

and I haven't been able to get it to work using a Fourier Transform (see below). My convolution is orders of magnitude larger than the input functions and a completely wrong shape.

Now, I have been getting all of my information about this from Google, so I must be missing something pretty obvious. Especially since both ways of doing it are failing for me. I know this isn't completely related to the MKL, but I plan on using the MKL wrappers for fftw after i get my test case working. If anyone could give me a hand with either the FT or regular convolution, I would really appreciate it. You can download the project from here

Convolution project

I know I'm using the FFTW correctly, as I can FT and inverse FT a function without a problem. Thanks for any help.



#include "stdafx.h"
#include 
#include 
#include 
#include "fftw3.h"

using namespace std;

int _tmain(int argc, _TCHAR* argv[])
{
 double *x = (double*)fftw_malloc(1000*sizeof(double));
 double *y = (double*)fftw_malloc(1000*sizeof(double));
 double *conv3 = (double*)fftw_malloc(1000*sizeof(double));
 int NumberofPoints = 0;
 
 fftw_plan p;
 fftw_complex* out1 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * 1000);
 fftw_complex* out2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * 1000);
 fftw_complex* final = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)* 1000);

 for(int i = 0; i < NumberofPoints; i++)
 {
  out1[0] = 0;
  out1[1] = 0;
  out2[0] = 0;
  out2[1] = 0;
  final[0] = 0;
  final[1] = 0;
 }

 //Read in our files
 ifstream stream("individgraphs.dat");

 if(stream.is_open() == false)
  cout << "Error" << endl;

 while(stream >> x[NumberofPoints]>>y[NumberofPoints])
 {
  NumberofPoints++;
 }
 stream.close();

 //FFTW

 //FT data file 1
 p = fftw_plan_dft_r2c_1d(NumberofPoint
s, x, out1, FFTW_ESTIMATE);
 fftw_execute(p); 
 fftw_destroy_plan(p);
 //FT data file 2
 p = fftw_plan_dft_r2c_1d(NumberofPoints, y, out2, FFTW_ESTIMATE);
 fftw_execute(p);
 fftw_destroy_plan(p);

 //Convolution - Complex multiplication and scaling
 for(int i = 0; i < NumberofPoints; i++)
 {
   final[0] = (out1[0]*out2[0]-out1[1]*out2[1])*(1.0/(double)NumberofPoints);
   final[1] = (out1[0]*out2[1]+out1[1]*out2[0])*(1.0/(double)NumberofPoints);
 }

 //Invert the product
 p = fftw_plan_dft_c2r_1d(NumberofPoints, final, conv3,FFTW_ESTIMATE);
 fftw_execute(p);
 fftw_destroy_plan(p);

 //Write output file
 ofstream test("outfile.txt");
 for (int i = 0; i < NumberofPoints; i++)
 {
  test << conv3/NumberofPoints << endl;
 }
 test.close();

 //Free memory
 fftw_free(x);
 fftw_free(y);
 fftw_free(out1);
 fftw_free(out2);
 fftw_free(conv3);

 getch();
 return 0;
}

topic Help with convolving data sets in Intel® oneAPI Math Kernel Library

Help with convolving data sets