What do I not understand about the fft and angle functions?

EDITED VERSION

I am trying to create an example showing how to use the fft function to approximate the Fourier transform. My example is to calculate the Fourier transform of h(t) = 2*pi*exp(-2*pi*t)u(t). The analytic solution is H(f) = 1/(1+jf). Here is my code:

N=1024; % sample count
T=10; % h(t) very small for t > T
t=linspace(0,T,N); % N sample points
h=2*pi*exp(-2*pi*t); % low pass IR 1 Hz cut-off
H=T*fft(h)/N; % FFT calculation of Fourier transform
Hmag=abs(H(1:N/2)); % Magnitude calculation
Hphase=180*angle(H(1:N/2))/pi; % Phase in degrees
f=(0:N/2-1)/T; % frequency bins thru Nyquist freq
subplot(2,1,1);
semilogx(f,20*log10(Hmag))
grid on
subplot(2,1,2);
semilogx(f,Hphase)
grid on

The results are shown in my response to the answer to an earlier version of my question. The magnitude plot is fairly close to the analytic result, but the phase plot goes a bit crazy at higher frequencies. Where am I going wrong?

To repeat, my goal is to show the relationship between the FFT of a discrete time function and the Fourier transform of its continuous time counterpart. I expected a closer relationship.

Answer 1

[Star Strider]

I have some problems following your code. I invite you to take the analytic approach to see if you get the same result. If you have the Symbolic Math Toolbox, try this:

syms w0 t w L real
h=w0*exp(-w0*t);
H1 = int(h *  exp(-j*w*t), t, 0, 2*pi);
H2 = simplify(H1, 'steps', 10);
H3 = rewrite(H1, 'sincos');
ReH3 = simplify(real(H3), 'steps', 10)
ImH3 = simplify(imag(H3), 'steps', 10)
phas = simplify(atan(ImH3/ReH3), 'steps', 10)
fphas = matlabFunction(phas)

The last result is:

fphas = @(w,w0)atan((w.*cos(pi.*w.*2.0)-w.*exp(pi.*w0.*2.0)+w0.*sin(pi.*w.*2.0))./(-w0.*cos(pi.*w.*2.0)+w0.*exp(pi.*w0.*2.0)+w.*sin(pi.*w.*2.0)));

Plug in values for the respective variables and see what you get.

[Warren Koontz]

I am trying to create an example showing how to use fft to calculate the Fourier transform. I have simplified the example to finding the Fourier transform of h(t) = 2*pi*exp(-2*pi*t)u(t), which should be H(f) = 1/(1+j*f). To do this I first create a sampled version of h(t):

N=1024; % sample count
T=10; % time limit
t=linspace(0,T,N); % N sample points
h=2*pi*exp(-2*pi*t); % low pass IR 1 Hz cut-off

I then calculate the Fourier transform using fft (normalizing by T/N) and calculate the magnitude and phase of the first N/2 components (the rest are redundant).

H=T*fft(h)/N; % FFT calculation of Fourier transform
Hmag=abs(H(1:N/2)); % Magnitude calculation
Hphase=180*angle(H(1:N/2))/pi; % Phase in degrees

Finally, I plot the magnitude (converted to dB) and the phase

f=(0:N/2-1)/T; % frequency bins thru Nyquist freq
subplot(2,1,1);
semilogx(f,20*log10(Hmag))
grid on
subplot(2,1,2);
semilogx(f,Hphase)
grid on

Here is the result:

Here is the result of just plotting H(f) = 1/(1+jf) for the same frequency range:

The magnitude plot compares well, at least until you get near the Nyquist frequency. The phase plot, however, gets weird beyond the cut-off frequency.

[Star Strider]

The calculation is not valid beyond the Nyquist frequency. You’re seeing aliasing above the Nyquist frequency, so all bets are off.

[Warren Koontz]

Yes, I have now done some reading about aliasing and leakage. The problem crops up in the phase plot well before the Nyquist frequency, but it seems that inaccuracies can show up early on. Further plots show that both the real and imaginary parts of H get very small even in relatively low frequency bins. Since the phase depends on the ratio of the imaginary to the real, the inaccuracies might be more severe in the phase plot. I guess that's the way it is. Thanks for your help!!

[Star Strider]

My pleasure!

The phase plot depends on the atan2 function, and any calculation involving the division of two uncertain floating-point numbers is going to be even more uncertain. (There are a number of computer science professionals here who may have some thoughts on this problem. I invite their contributions.)