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.
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?