MATLAB Answers


IFFT integration doesn't match cumtrapz integration

Asked by Nik Arrivederci on 7 Jun 2019
Latest activity Edited by David Goodmanson on 12 Jun 2019
For class, we had to differentiate a piecewise function using diff() and integrate it using cumtrapz() over a defined interval, then do the same through ifft and fft of the originial function using a method given to us by the prof. Here's the entirety of the code relevant to the question
function [] = class4()
X = -10:0.05:10;
N = length(X);
Y = zeros(1,N);
for i = 1:N
if abs(X(i))<=1
Y(i) = exp(-(abs(sin(X(i))))^2);
Y(i) = 0;
x = [n:h:m-h];
Nx = size(x,2);
k = 2*pi/(m-n)*[0:Nx/2-1 0 -Nx/2+1:-1];
IF = ifft(fft(F)./(j*(k+dirac(k))));
plot(X, IY, x, IF);
grid on;
legend('Integral','Fourier Integral');
I've omitted the diff, because it kinda works fineish, but the integral part produces very different results -- IY/cumtrapz gives a steadily rising curve from -1 to 1, while IF/ifft gives either a weirdly shaped curve that dips below y = 0, then grows from -1 to 1 and steadily declines afterwards or a flat line at y = 0. This depends on the inclusion of dirac(k), which prof told us to do, leaving us to figure out where and how to place it . By playing around with coefficients in ifft and fft, I got this exaggerated curve so that you could see what it looks like compared to the cumtrapz integral.
Enough rambling, I guess. Question: Is this normal? If not, what do I do?
I tried just about every placement of dirac(k) I could think of, and prof explained it's mandatory because of problematic integration at k = 0, which is included in the frequency domain.
After adding a linspace(n,m,(m-n)/h+1) to cumtrapz, I got two almost identical graphs, with one of them shifted downwards. Shifting the ifft upwards gave me this.
Now I lean heavily towards this being a property of Fourier transforms, but there's still a possibility I'm doing something wrong, so the question remains the same.


John D'Errico
2019 年 6 月 7 日
Do you really expect a trapezoidal rule integration to be exact?
Nik Arrivederci 2019 年 6 月 8 日
Not really, but I'm much more bothered by the Fourier integral having those weird appendages than by their heights not matching up.
Also, I've modified the code a bit and even after retrying this on a smaller range with a much (0.0001 vs 0.05) step, the appendages are still there.

Sign in to comment.




1 Answer

回答者: David Goodmanson
2019 年 6 月 12 日
編集済み: David Goodmanson
2019 年 6 月 12 日

Hi NIk,
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), ... G(N)-G(N-1)]/delx.
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;
% define the function,
% shift it the the right to provide more generality
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
% insert a negative dirac pulse at the end of g
Idef = sum(g)*delx; % definite integral
g1 = g;
f = max(find(abs(g)>0));
g1(f+1) = -Idef/delx; % dirac ~ 1/delx
plot(x,g1); grid on
% IIdef is related to the second integral of g
IIdef = sum(cumsum(g1)*delx)/delx; % (delx cancels out)
% cumsum integration for comparison purposes
Gcumsum = cumsum(g1)*delx;
% frequency grid with 0 at the center of the array
K = 2*pi/delx;
k = K*(-(N-1)/2:(N-1)/2)/N;
% case 0, no dirac pulse
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
title('case 0')
% case 1, dirac pulse
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
title('case 1')
% case 2, dirac pulse but using 1/i*k factor
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
title('case 2')


Sign in to comment.