# Description:
# Computes a two dimensional Fourier transform.
#
clear
# total time, g(t1, t2) is "small" for |(t1, t2)| > T / 2
T = 4.   
# number of data points, dt resolves time integrals
N = 16   
# time between data values
dt = T / N       
# time between frequency values
df = 1 / T       
# time grid
t = (seq(N) - N / 2 - 1) * dt 
# exponent in definition of h(t)
alpha = 5. 
# h(t) = exp( - alpha * |t|)
h = exp( - alpha * abs(t)) 
# g(t1, t2) = h(t1) * h(t2)
g = fillcols(h, N) % fillrows(h', N) 
#
# Plot g near the origin.
gtitle("g(t1, t2)")
surface(g)
gupdate
#
# the continuous transform of g(t1, t2) is given by
# G(f1, f2) = H(f1) * H(f2)
# where H(f) = 2 * alpha / ((2 * PI * f)^2 + alpha^2)
# is the discrete transform of g
G = real(fft2d(complex(g))) * dt * dt 
gaddwin
gtitle("G(f1, f2)")
surface(G)