
FFT.OMS
Script File:
# 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")
Output:

