Prev Next dft2d

Two Dimensional Discrete Fourier Tranform
Syntax y = dft2d(x)
Include: include spt\dft2d.oms
See Also fft , ifft , idft2d

ARGUMENTS:  
   INPUTS:
      x = MATRIX, any numerical type. Data to be transformed.
   RETURN: MATRIX, COMPLEX. Forward transformed data.

Description

Perform FORWARD two-dimensional Discrete Fourier Tranform on a matrix.

Argument 'x' may be a matrix of any type. The intrinsic transform 'dft()' is used in this function allowing arbitrary matrix dimensions. Fast fourier transforms are performed if the only prime factors of the number of rows and columns are 2, 3, 5, and 7. Otherwise, the transform will be much slower for large matrices. Compare this to the 'fft2d()' function that uses the zero-centered 'fft()' routine.

Example


include spt\spthead.oms

#################################################################
# EXAMPLE: FIR filter 2-dimensional rectangular pulse and
# display 2D spectrum.
#################################################################

#--- Define x/y dimensions ---

N = 16;
M = 16;

#--- Create Input waveform ---
# Make a 2-D rectangular pulse

d  = 2;      # pulse width in samples, same for x and y direct.
x1 = double(seq(N));
x1 = int( (x1>=(N/2-(d-1))) and (x1<=(N/2+d)) );
x2 = double(seq(M)');
x2 = int( (x2>=(M/2-(d-1))) and (x2<=(M/2+d)) );
x  = x1*x2;  # 2-D input signal, a NXM matrix

#--- Compute input spectrum using 2-D DFT ---

X         = dft2d(x);        # compute DFT
XdB       = db20(X);         # convert to dB for display
XdB       = XdB - max(XdB);  # normalize to 0 dB
XMINDB    = -30d0;           # min clip level
clipem    = (XdB>=XMINDB);   # start clipping
notclipem = not clipem;
XdB       = int(clipem)%XdB + int(notclipem)*XMINDB;

# Swap upper/lower halves
a                      = XdB.blk(1,1,N/2,M);
XdB.blk(1,1,N/2,M)     = XdB.blk(N/2+1,1,N/2,M);
XdB.blk(N/2+1,1,N/2,M) = a;
# Swap left/right halves
a                      = XdB.blk(1,1,N,M/2);
XdB.blk(1,1,N,M/2)     = XdB.blk(1,M/2+1,N,M/2);
XdB.blk(1,M/2+1,N,M/2) = a;

#--- Make a 2-D LOWPASS FILTER   ---
#--- Use linfir() to make filter ---

fs  = 1d0;        # sampling rate
fN1 = .1d0;       # lower cutoff, FIR #1
fN2 = .49d0;      # upper cutoff, FIR #1
fM1 = .1d0;       # lower cutoff, FIR #1
fM2 = .49d0;      # upper cutoff, FIR #1
hN  = linfir(fN1,fN2,fs,N,"LOWPASS","BLACKMAN");
hM  = linfir(fM1,fM2,fs,M,"LOWPASS","BLACKMAN");
h   = hN * hM';   # 2-D filter, a NXM matrix

#--- Compute FIR frequency response using 2-D DFT ---

H = dft2d(h);

#--- Convolve x and h by multiplying transforms ---

Y = H%X;          # spectrum of filtered input

# Swap upper/lower halves
a = Y.blk(1,1,N/2,M);
Y.blk(1,1,N/2,M) = Y.blk(N/2+1,1,N/2,M);
Y.blk(N/2+1,1,N/2,M) = a;
# Swap left/right halves
a = Y.blk(1,1,N,M/2);
Y.blk(1,1,N,M/2) = Y.blk(1,M/2+1,N,M/2);
Y.blk(1,M/2+1,N,M/2) = a;

# Convert output spectrum to dB and nomalize to 0dB MAX.
# also clip data to -100dB MIN

YdB       = db20(Y);
YdB       = YdB - max(YdB);
clipem    = (YdB>=XMINDB);
notclipem = not clipem;
YdB       = int(clipem)%YdB + int(notclipem)*XMINDB;

#--- Compute filtered waveform by inverse 2-D DFT ---

y = double(idft2d(Y));  # for real inputs only
y = y/max(y);       # normalize to unity for display ease

#--- Swap halves to undo FIR delay ---
a                    = y.blk(1,1,N/2,M);
y.blk(1,1,N/2,M)     = y.blk(N/2+1,1,N/2,M);
y.blk(N/2+1,1,N/2,M) = a;
a                    = y.blk(1,1,N,M/2);
y.blk(1,1,N,M/2)     = y.blk(1,M/2+1,N,M/2);
y.blk(1,M/2+1,N,M/2) = a;



Reference

Oppenheim and Schafer, "Digital Signal Processing.", Prentice-Hall.