1
$\begingroup$

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.

Result from the convolution using FFTW normalized by a constant factor compared with the expected result. Here someConstantFactor = 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!

$\endgroup$
3
  • 1
    $\begingroup$ 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, $\endgroup$
    – Hilmar
    Commented Mar 19, 2023 at 12:49
  • $\begingroup$ @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? $\endgroup$
    – gebegb
    Commented Mar 19, 2023 at 14:59
  • $\begingroup$ 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. $\endgroup$
    – Hilmar
    Commented Mar 19, 2023 at 15:31

1 Answer 1

1
$\begingroup$

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.