6

I have adapted the code in Comparing FFT of Function to Analytical FT Solution in Matlab for this question. I am trying to do FFTs and comparing the result with analytical expressions in the Wikipedia tables.

My code is:

a = 1.223;
fs = 1e5; %sampling frequency
dt = 1/fs;
t = 0:dt:30-dt;     %time vector
L = length(t); % no. sample points
t = t - 0.5*max(t); %center around t=0

y = ; % original function in time
Y = dt*fftshift(abs(fft(y))); %numerical soln

freq = (-L/2:L/2-1)*fs/L; %freq vector
w = 2*pi*freq; % angular freq

F = ; %analytical solution

figure; subplot(1,2,1); hold on
plot(w,real(Y),'.')
plot(w,real(F),'-')
xlabel('Frequency, w')
title('real')
legend('numerical','analytic')
xlim([-5,5])
subplot(1,2,2); hold on;
plot(w,imag(Y),'.')
plot(w,imag(F),'-')
xlabel('Frequency, w')
title('imag')
legend('numerical','analytic')
xlim([-5,5])

If I study the Gaussian function and let

y = exp(-a*t.^2); % original function in time

F = exp(-w.^2/(4*a))*sqrt(pi/a); %analytical solution

in the above code, looks like there is good agreement when the real and imaginary parts of the function are plotted:

enter image description here

But if I study a decaying exponential multiplied with a Heaviside function:

H = @(x)1*(x>0); % Heaviside function
y = exp(-a*t).*H(t);

F = 1./(a+1j*w); %analytical solution

then

enter image description here

Why is there a discrepancy? I suspect it's related to the line Y = but I'm not sure why or how.

Edit: I changed the ifftshift to fftshift in Y = dt*fftshift(abs(fft(y)));. Then I also removed the abs. The second graph now looks like:

enter image description here

What is the mathematical reason behind the 'mirrored' graph and how can I remove it?

Medulla Oblongata
  • 3,771
  • 8
  • 36
  • 75
  • The Fourier Transform, although closely related, is not a Discrete Fourier Transform (implemented via the FFT algorithm). So, under some specific conditions you may get very close results, but quite often you will get noticeable differences (just as you observed) – SleuthEye Mar 16 '18 at 11:31
  • Also by taking `abs` you eliminate the imaginary part, so that'll always be 0. I think you need `fftshift` instead of `ìfftshift`. – ViG Mar 16 '18 at 11:37
  • @Zep This is just a MWE - ultimately I need to do FFT and IFFT of much more complicated functions that don't have an analytical form. So I'm doing something like `FFT{ IFFT[F(w)]*IFFT[G(w)] / IFFT[J(w)] }` where I only know the vectors for `F,G,J` in frequency. – Medulla Oblongata Mar 16 '18 at 11:56
  • @Medulla Oblongata I removed my comment because I thought that `ifftshift` would shift AND inverse transform, instead it only shifts. – Zep Mar 16 '18 at 12:00
  • I guess the first experiment, with the Gauss, worked because of the `abs` call? It removes the phase component, which is computed wrong. In the second plot, you see the effect of the `abs`: no imaginary component, purely positive real part. When you removed the `abs` in the 3rd plot, you see the phase is all wrong, but the magnitude is correct. – Cris Luengo Mar 17 '18 at 02:17

1 Answers1

5

The plots at the bottom of the question are not mirrored. If you plot those using lines instead of dots you'll see the numeric results have very high frequencies. The absolute component matches, but the phase doesn't. When this happens, it's almost certainly a case of a shift in the time domain.

And indeed, you define the time domain function with the origin in the middle. The FFT expects the origin to be at the first (leftmost) sample. This is what ifftshift is for:

Y = dt*fftshift(fft(ifftshift(y)));

ifftshift moves the origin to the first sample, in preparation for the fft call, and fftshift moves the origin of the result to the middle, for display.

Edit

Your t does not have a 0:

>> t(L/2+(-1:2))
ans =
  -1.5000e-05  -5.0000e-06   5.0000e-06   1.5000e-05

The sample at t(floor(L/2)+1) needs to be 0. That is the sample that ifftshift moves to the leftmost sample. (I use floor there in case L is odd in size, not the case here.)

To generate a correct t do as follows:

fs = 1e5; % sampling frequency
L = 30 * fs;
t = -floor(L/2):floor((L-1)/2);
t = t / fs;

I first generate an integer t axis of the right length, with 0 at the correct location (t(floor(L/2)+1)==0). Then I convert that to seconds by dividing by the sampling frequency.

With this t, the Y as I suggest above, and the rest of your code as-is, I see this for the Gaussian example:

>> max(abs(F-Y))
ans =    4.5254e-16

For the other function I see larger differences, in the order of 6e-6. This is due to the inability to sample the Heaviside function. You need t=0 in your sampled function, but H doesn't have a value at 0. I noticed that the real component has an offset of similar magnitude, which is caused by the sample at t=0.

Typically, the sampled Heaviside function is set to 0.5 for t=0. If I do that, the offset is removed completely, and max difference for the real component is reduced by 3 orders of magnitude (largest errors happen for values very close to 0, where I see a zig-zag pattern). For the imaginary component, the max error is reduced to 3e-6, still quite large, and is maximal at high frequencies. I attribute these errors to the difference between the ideal and sampled Heaviside functions.

You should probably limit yourself to band-limited functions (or nearly-band-limited ones such as the Gaussian). You might want to try to replace the Heaviside function with an error function (integral of Gaussian) with a small sigma (sigma = 0.8 * fs is the smallest sigma I would consider for proper sampling). Its Fourier transform is known.

Cris Luengo
  • 55,762
  • 10
  • 62
  • 120
  • You're right about the `ifftshift`, and the plot looks ok now. But when I test this line of code with the Gaussian function, the FFT also gives a small imaginary part in the order of 10^(-6). Why is it there and is there a way to remove it? – Medulla Oblongata Mar 18 '18 at 10:04
  • @MedullaOblongata: I would guess it's numerical inaccuracy. Or maybe it's a very small shift in the computation of the input function. Do you have an actual `t==0` in your time axis samples? Calculating `t` using `dt` increments introduces errors, and `t-0.5*max(t)` is probably not precise either. Try making a `t` vector with integers, and dividing by `fs`. Make sure that `t(floor(L/2)+1)==0`. – Cris Luengo Mar 18 '18 at 14:05
  • See here for the difference between the colon operator an `linspace`: https://stackoverflow.com/q/26292695/7328782 – Cris Luengo Mar 18 '18 at 14:15
  • @MedullaOblongata: I've updated the answer with the correct form to construct the time axis `t`. – Cris Luengo Mar 19 '18 at 17:14