# Description:
# Uses the function fft to approximate the Fourier transform.
#
# The Fourier transform ghat of g is defined by ghat(f) =
# integral from - inf to + inf of g(t) * exp( - 2 * pi * 1i0 * t) * dt
clear
#
# total time (so that g(t) is small for |t| > T / 2)
T = 8. 
# number of data points (so dt resolves time integral)
N = 256 
# half of N
N2 = N / 2 
# time between data values
dt = T / N 
# time between frequency values
df = 1 / T 
# time grid
t = (seq(N) - N2 - 1) * dt 
# frequency grid
f = (seq(N) - N2 - 1) * df 
# exponent in definition of g(t)
alpha = 5. 
# g(t) = exp( - alpha * |t|)
g = exp( - alpha * abs(t)) 
# ghat(f) = Fourier transform of g(t)
ghat = 2 * alpha / ((2 * PI * f)^2 + alpha^2) 
# approximation to ghat(f) integral
G = fft(complex(g)) * dt 
#
# just plot the positive frequencies
f = f.row(N2 + 1, N2)
ghat = ghat.row(N2 + 1, N2)
G = G.row(N2 + 1, N2)
#
# set axis to plot positive frequencies
gxaxis("linear", 0., 5, 5, 10) 
# limit as f - > 0 of ghat(f) = .4
gyaxis("linear", 0., .5, 5, 10) 
# plot continuous transform
gplot(f, ghat, "solid") 
# plot real part of approximation
gplot(f, real(G), "cross") 
#
# Label the plot
point = gcoord(3, .4)
gaddtext( {"Curve is true transform", "x is discrete transform"}, point )
gxtitle("frequency")
gytitle("transform")