3
$\begingroup$

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}]; 
  curve]

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

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?

Improve this question
$\endgroup$
1
  • 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, 2015 at 14:49

1 Answer 1

Reset to default
4
$\begingroup$

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} $$

Improve this answer
$\endgroup$
1
  • $\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, 2020 at 2:17

Your Answer

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

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