Help with convolving data sets

Help with convolving data sets

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.

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

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.

[CODE]

#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[i][0] = 0;
out1[i][1] = 0;
out2[i][0] = 0;
out2[i][1] = 0;
final[i][0] = 0;
final[i][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[i][0] = (out1[i][0]*out2[i][0]-out1[i][1]*out2[i][1])*(1.0/(double)NumberofPoints);
final[i][1] = (out1[i][0]*out2[i][1]+out1[i][1]*out2[i][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[i]/NumberofPoints << endl;
}
test.close();

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

getch();
return 0;
}
[/CODE]

1 post / 0 new
For more complete information about compiler optimizations, see our Optimization Notice.