• Matrix representation of DFT
• Orthonormal nature of M
• How FFT works
• Significance of phase
• Gaussian filter

The Discrete Fourier Transform (DFT) can also be summarized as a matrix operation. The DFT formula:

Translates into the matrix operation MX, where X is the discrete sample data length N and the square matrix (NxN) M is defined by the Matlab function:
 
 

This operation is verified by getting identical results to Matlab's FFT
X = [ 1 2 3 4 3 2 1 0]
 
 
 

Matlab can be used to demonstrate that M is orthonormal, which requires showing (1) at least one dot product of different vectors from the M vector space is zero; and (2) the length of eacah vector is 1 (dot product of vector with itself). Meeting the second requirement required normalizing M; I believe we can normalize here because there is no definition for normalization. The only requirement is that 1/N is multiplied from fft() to ifft() . The following line is added the last line of the MDFT function:

M = M./sqrt(N)

Here are examples that work for all cases. The results are within the tolerance to be zero.
 

More evidence of the orthogonality of M can further be verified by graphing the real and imaginary components. The real component of M is an even function, which means the function is symmetrical about the Y-axis (see the corners of the center star). This is a result of the cos() function (even function) which creates the real component. The imaginary component is an odd function, which means it is symmetrical about x=. This result is from the sin() function (odd function) that creates the imaginary component (compare the same corners of the star as the real component).

Figure: REAL component M

Figure: IMAG component M
 
 

Resource: Dirty Pixels by Jim Blinn
 

How the FFT works

As the name Fast Fourier Transform indicates, the FFT is a faster algorithm to find the Discrete Fourier Transform (DFT) of a sampled signal. This is a high level overview how the FFT works. I do not claim to understand the all the subtleties that accomplish the final FFT.

The FFT algorithm takes the "divide and conquer" approach to solving the program similar to how Merge sort works. The divide phase divides the samples into one-sample long signals. However, there is a unique reordering of the samples: The position (index) of every sample is swapped. The new location is the bit reversal of the sample index; for instance sample 3 (0x011), would swap places with 6 (0x110). This arrangement setups the samples for the combination phase. It is amazing that by the end of the division process, the individual time samples are now frequency values.

The division process is similar the merge process of Merge-sort except that the combine phase is where the computations of the FFT take place. Samples are combined by interlacing the data values of two groups that are merged. The interlacing of the samples corresponds to complex multiplication. The interlacing is followed by a "butterfly" calculation. I believe the butterfly operation calculates the complex conjugate, which gives the Fourier Transforms' graph its symmetry.

The FFT growth function is O(n lgn) verses the O(n^2) of DFT. The FFT growth is just like Merge-sort, divide and conquer. The n^2 nature of DFT can be seen in the matrix operation in the previous section.

Resources:

• DSP Guide Ch 12 
• FFT Tutorial 

This section demonstrates that the defining details of a picture are in the phase component of a 2D-Fourier Transform of an image. Here two 256x256 images are loaded and the magnitude and phase of the images is swapped.

Concluding from the results, the phase information defines the picture. This could be explained by considering the final image is constructed by superposition of sinusoidal waves. Changing the magnitude will only efect magnitude and intensity of the picture. However, offsetting the waves will create interference patterns. It is within these patterns that the images are defined.
 
 

Original Image
» t1 = imread('test1.jpg');
Image Fourier Transform magnitude plot
» t1 = imread('test1.jpg'); » colormap hsv » f1 = abs(fftshift(fft2(t1))); » imagesc(f1)
Swapping magnitude and phase
» f1 = fft2(t1); » f2 = fft2(t2); » a1 = angle(f1); » a2 = angle(f2); » m1 = abs(f1); » m2 = abs(f2); » imagesc(abs(ifft2(complex(m2.*cos(a1), m2.*sin(a1)  ))))	» imagesc(abs(ifft2(complex(m1.*cos(a2), m1.*sin(a2)  ))))
Resulting images

Gaussian

Designing a Gaussian filter that should suppress 1 cycle in 4 pixels to 5%. The following picture defines the relationship of the spatial and frequency domain. The problem is to find the std for the Gaussian that is at 5% when the frequency is at 0.25Hz. The frequency domain equation is easier to solve because I has one less std. So translating 0.25Hz in the spatial domain to frequency has to be in rad/sec as can be seen by "w" (thanks to Miguel for pointing this out). So the  frequency is 2*pi*0.25. Solving for std (spatial) in frequency domain Gaussian at 1.57 equal to 5%: sqrt(-2*ln(0.05)/w^2) = 1.56 std.
 

Gaussian surface

» x= -7.5:1:7.5; » y=x; » [X Y] = meshgrid(x,y); » Z = (1./sqrt(2*pi*1.56)).*exp(-( Y.^2 + X.^2)./(2.*1.56.^2)); » surf(Z)

Gaussian Picture

» imagesc(Z);

Image

» imagesc(t1)

Filtered image

» imagesc(filter2(Z,t1))