This is an Interesting problem. Let g be the function and G its integral. If you just want to key in on an improved result, keep the common code before case 0 and go to the code for case 1 (fig 3) below. That plot shows (assuming default colors) the waveform g in yellow, its integral G in blue and a non-fft numerical integral G for comparison purposes in red. The two G’s are basically overlays.
For comparison purposes it’s possible to create an fft method that will agree with G = cumtrapz(g), but it’s a whole lot easier to agree with the approximate integral G = cumsum(g)*delx where delx is the x array spacing. That’s what is done here.
As to why it works, my apologies for what will take quite a few paragraphs. First off, assume N points, N odd, with x array spacing = delx.
For simplicity, assume there are some zeros on each side of the integrand g, as you have done (technically the integrand includes the zeros but by ‘the integrand’ in this context I mean the nonzero part). The elements of g = dG/dx are just first differences of G,
g = [G(1)-G(N), G(2)-G(1), ...
Since an fft array is periodic it wraps around, so G(1) has the end point G(N) subtracted from it. It’s an array of circular differences. Note that if g is summed, every element of G is mentioned twice and you get mean(g) = 0. The fft method necessarily assumes mean(g) = 0.
With the method you used, g is effectively pushed down to have mean 0. That means there is a negative baseline as shown in the red plot in fig 1a. When that gets integrated, it produces a line with negative slope in G. The fft method cleverly creates a sloped line with no discontinuity (i.e. there is no large jump in value between the left and right end points of the (periodic) array), so you end up with the result in case 0 (fig 2). It’s actually an interesting result but of course it’s not what you want.
One way to deal with a function whose mean is nonzero is to insert a negative dirac pulse at the first array point just after the end of the integrand g. This can make the overall mean of g equal to zero and still preserve a baseline of zero. The function g1 in the code includes a dirac pulse.
With Matlab dirac, one has to go look up what the definition is. A dirac pulse of height inf as defined in Matlab is appropriate for a continuous variable. But for a discrete array, a dirac pulse of height 1/delx is more useful.
Fig 1b shows g1, the integrand g plus the appropriate dirac pulse.
Another issue is division by i*k in the frequency domain. Take a look at the derivative of G. From each element of G we subtract the element to its left. A translation by delx in G corresponds to a phase factor of exp(-i*k*delx) in the fourier domain. Including the element itself, the overall factor to create the derivative is
To integrate, you divide by that factor instead of multiplying. This is exact.
When k is small, the factor reduces to i*k*delx and you end up dividing by that. Later in the process the delx is canceled out, so you effectively divide by i*k. But that is just a small k approximation. (It’s exact for the continuous case, which we don’t have). Case 2 (fig 4) shows what happens if you keep the dirac pulse but use the small k approximation. As you can see the result is not too bad but the result picks up some high frequency wiggles.
In both exact and small-k approximation there is the issue of what to do when dividing at k=0. I won’t spend much time on this but the value at k=0 in the k domain is directly related to the mean value of the function in the x domain. These values are set appropriately below.
x = -10:0.05:10;
delx = x(2)-x(1);
N = length(x);
Ncenter = (N+1)/2;
g = zeros(1,N);
xs = x-2;
ind = abs(xs) <=1;
g(ind) = exp(-(abs(sin(xs(ind)))).^2);
plot(x,g,x,g-mean(g)); grid on
Idef = sum(g)*delx;
g1 = g;
f = max(find(abs(g)>0));
g1(f+1) = -Idef/delx;
plot(x,g1); grid on
IIdef = sum(cumsum(g1)*delx)/delx;
Gcumsum = cumsum(g1)*delx;
K = 2*pi/delx;
k = K*(-(N-1)/2:(N-1)/2)/N;
zg = fftshift(fft(g));
zG = zg./(1-exp(-i*k*delx));
zG(Ncenter) = 0;
G = ifft(ifftshift(zG))*delx;
plot(x,G,x,g); grid on
zg1 = fftshift(fft(g1));
zG1 = zg1./(1-exp(-i*k*delx));
zG1(Ncenter) = IIdef;
G1 = ifft(ifftshift(zG1))*delx;
plot(x,G1,x,Gcumsum,'o',x,g); grid on
zg1 = fftshift(fft(g1));
zG2 = zg1./(i*k*delx);
zG2(Ncenter) = IIdef;
G2 = ifft(ifftshift(zG2))*delx;
plot(x,G2,x,G1,x,g); grid on