Prev Next idft2d

REVERSE Two-Dimensional Discrete Fourier Transform
Syntax y = idft2d(m,ka,fc,pc,fs)
Include: include spt\idft2d.oms
See Also fft , ifft , dft2d

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

Description

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

Argument 'x' may be a matrix of any type. The intrinsic transform 'idft()' is used in this function allowing arbitrary matrix dimensions. Compare this to the 'ifft2d()' function that uses the zero-centered 'ifft()' 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.