Convolution of two functions using FFTW

Question

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 */
    double f[zero_padded_size] = {};
    double g[zero_padded_size] = {};

    /* 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;
    F = (fftw_complex *) fftw_malloc(sizeof(fftw_complex)*zero_padded_size);
    G = (fftw_complex *) fftw_malloc(sizeof(fftw_complex)*zero_padded_size);

    /* 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;
    H = (fftw_complex *) fftw_malloc(sizeof(fftw_complex)*zero_padded_size);

    /* 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)] */
    double h[zero_padded_size] = {};
    /* 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 */
    double normFactor = 1./zero_padded_size;
    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!

Answer 1

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

Why would it be 1?

In fact,

$$\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))