Recently I tried to compare results for 1D direct convolution and convolution via FFT. I expected to get absolutely the same result, however I faced with a problem that results are different, especially near the left boundary.

enter image description here

enter image description here

FFT convolution should be normalized, however it doesn't change the difference near the left boundary. As I understand this difference appears due to the fact that FFT provides circular convolution, while the direct convolution is linear. I have found that it can be fixed by Overlap–add method, however as I understand it should be used when curve and kernel have the different length, but in my case both of them have equal length.

Mathematica code:

mat = Table[0, {x, 512}, {y, 512}];
SetAttributes[FillMatrix, HoldAll]
FillMatrix[Kep_, Ktrans_, mat_] := 
  Do[Do[mat[[j, i]] = Exp[-Kep (j - i + 1)], {i, 1, j, 1}], {j, 1, 
    Length[mat], 1}];

FillArray[Kep_, Ktrans_, mat_] := 
 Module[{curve = {}}, 
  Do[curve = Append[curve, Exp[-Kep*i]], {i, 1, Length[mat], 1}]; 

AIF = {};
SetAttributes[FillAIF, HoldAll]
FillAIF[AIF_, n_] := 
 Do[AIF = Append[AIF, PDF[GammaDistribution[3, 2], i*0.1]], {i, 1, n, 

FillAIF[AIF, Length[mat]]

FillMatrix[0.005, 1, mat]
conv = mat.AIF;

AIFFourier = Fourier[AIF];
kernelCurve = FillArray[0.005, 1, mat];
kernelCurveFourier = Fourier[kernelCurve];

resSource = InverseFourier[AIFFourier*kernelCurveFourier];
ListPlot[resSource, Joined -> True]
ListPlot[conv, Joined -> True]

Question: Why does this difference appear and how can it be fixed?

  • 1
    $\begingroup$ The Fourier transform implicitly assumes that the function you are transforming is periodic; hence the observed difference. The usual remedy for nonperiodic functions is zero-padding, as Kirill has pointed out. $\endgroup$ – Christian Clason Feb 27 '15 at 14:49

A brute-force approach is to pad the signal and the kernel with enough zeros that there would be no aliasing when doing circular convolution:

AIFFourier = Fourier[PadRight[AIF, 1024]];
kernelCurve = FillArray[0.005, 1, mat];
kernelCurveFourier = Fourier[PadRight[kernelCurve, 1024]];

resSource = InverseFourier[AIFFourier*kernelCurveFourier];
ListPlot[resSource[[1 ;; 512]], Joined -> True]
ListPlot[conv, Joined -> True]
StandardDeviation[resSource[[1 ;; 512]]/conv]

$$ \Longrightarrow 2.37134\times10^{-15} $$

| cite | improve this answer | |
  • $\begingroup$ The term enough is not exact. I believe that the padding should ensure that the signals are padded before computing the FFT so that they become of length at least Signal1 + Signal2 -1 to avoid time aliasing $\endgroup$ – VMMF Feb 21 at 2:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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