# Convolution of two functions using FFTW

I'm trying to perform a discrete convolution of two functions, $$f(x) = 1$$ and $$g(x) = \exp(-x)$$ of length nsize using FFTW. I have followed the procedure for zero-padding my data from here. My resulting values are qualitatively similar, but the magnitudes are way off from the expected result. I have attached the code I have written for reference.

#include <iostream>
#include <cmath>
#include <fstream>
#include "fftw3.h"

int main()
{
/* Size of initial arrays (f(x) and g(x)) arrays */
const int nsize = 100;
/* Size of the convoluted array (size(f(x)) + size(g(x)) - 1) */
const int zero_padded_size = 2*nsize - 1;

/* Declaration of the functions */

/* Initialization of f(x) and g(x), with zero-padding beyond their size */
for(int i = 0; i < nsize; i++)
{
f[i] = 1.0;
g[i] = exp(-2.*M_PI*i/nsize);
}

/* Declaration and memory allocation of complex array to store the fourier transforms of f(x) and g(x) */
fftw_complex *F;
fftw_complex *G;

/* FFTW plans to convert f(x) and g(x) to their fourier domain functions */
fftw_plan realfToFourierF = fftw_plan_dft_r2c_1d(zero_padded_size, f, F, FFTW_ESTIMATE);
fftw_plan realgToFourierG = fftw_plan_dft_r2c_1d(zero_padded_size, g, G, FFTW_ESTIMATE);
fftw_execute(realfToFourierF);
fftw_execute(realgToFourierG);

/* Declaration of complex array to store fourier-space multiplication of F[f(x)] and F[g(x)] */
fftw_complex *H;

/* Point-wise multiplication of F[f(x)] and F[g(x)] */
for(int i = 0; i < zero_padded_size; i++)
{
H[i][0] = (F[i][0] * G[i][0] - F[i][1] * G[i][1]);
H[i][1] = (F[i][0] * G[i][1] + F[i][1] * G[i][0]);
}

/* Array declaration for real-space values after Inverse FT of F[h(x)] = F[f(x)]*F[g(x)] */
/* FFTW plan for IFT of F[h(x)] to h(x) */
fftw_plan fourierHToRealh = fftw_plan_dft_c2r_1d(zero_padded_size, H, h, FFTW_ESTIMATE);
fftw_execute(fourierHToRealh);

/* Normalization factor */
std::ofstream fileout("./convolve.dat");
for(int i = 0; i < zero_padded_size; i++)
{
h[i] *= normFactor;
double t = 2.*M_PI*i/nsize;
fileout<<i<<"\t"<<t<<"\t"<<f[i]<<"\t"<<g[i]<<"\t"<<h[i]<<"\t"<<(1.0 - exp(-t))<<std::endl;
}
fileout.close();
}


I have attached the result obtained from the FFTW convolution against the expected result, and they differ by a constant factor (someConstantFactor) of 32.681202913.

I would be grateful if someone can provide any suggestions as to why this difference occurs and what I'm doing wrong. Thanks!

• Take out all this normalization stuff (which I honestly don't understand and there is not a single comment in your entire code). Before normalization h and g should be identical (within numerical precision). That's the end of the convolution. Whatever happens after this, has nothing to do with the algorithm, Mar 19 at 12:49
• @Hilmar: I have added comments to my code. I hope it is better now. By definition, $h(x) = \int_0^{x} g(t) dt = 1 - \exp(-x)$. They wouldn't be similar in this case. They differ by a constant (Here, $h(x) + g(x) = 1$). But in my case, $h(x)$ and the result of my convolution differ by a factor of 32.681202913. Could you please elaborate on what you mean by identical? Mar 19 at 14:59
• Oops, my bad, I misunderstood your text and code. So you are convolving the a Gaussian with a rectangular impulse response of length 100. You would expect the "flat" part of the amplitude to be the sum of the Gaussian (to be precise the first 100 samples of the Gaussian), i.e. 16.39 or thereabouts. There is no reason for it to be 1. Why would it ? Are you trying to implement an integrator here? Your result is correct but if it doesn't match the expectation you need to explain where your expectation comes from. Mar 19 at 15:31

As Hilmar says in the comments, the sum of $$g$$ is 16.39 or so.

Why would it be 1?

$$\int_0^x g(t) dt = \sinh(x) - \cosh(x) + 1 = 1 \mbox{ for large } x$$

but

$$\sum_{k=0}^{99} e^{-2\pi k/N_{\tt size}} \approx 16.39$$

so

$$\int_0^x g(t) dt \not = \sum_{n=0}^{N_{\tt size}-1} g[n]$$

The assumption is that

$$g(t) = e^{-t}$$

but that

$$g[n] = e^{-2\pi n/N_{\tt size}}$$

## Python code to calculate sum

import matplotlib.pyplot as plt
import numpy as np

nsize = 100
g = np.exp(-2.*np.pi*np.arange(nsize)/nsize)

print(np.sum(g))