Help with convolution algorithm

[diehardii]

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[j]= 0.0; 
      for (i = 0; i <= j; i++)
      {
         y[j] += x[i] * 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. 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 <fstream>
#include <iostream>
#include <conio.h>
#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(NumberofPoints, 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;
}

[yiannakop]

First of all, when convolving 2 sequences of length M and L the output sequence has length M+L-1. Here is an example of convolving two sequences in C:

for (n = 0; n < L+M; n++)
{
    y[n] = 0;
    for (m = max(0, n-L+1); m <= min(n, M); m++)
        y[n] += h[m] * x[n-m];
}

Anyway, why you think that your result is not correct? Give an example if possible.

Regards,
Theodore