The Discrete Fourier Transform (DFT) and lowpass filters

For this lab, you may wish to download /ece431/labs/lab2copy.tgz which contains all of the files needed, including the lab2.htm file you're reading now.

1. Continuous filters

In continuous frequency, one example of an optimal filter is the Gaussian (optimal in the sense that the product of time variance and frequency variance is minimum; if you'd like to learn more on this, see Gabor 1946, Hiesenberg 1926 ):

2. Discrete filters (example: rectangular function) and zero padding in the time domain

Consider a discrete time signal: x[n]=1, 0<=n<=3; x[n]=0, otherwise. 

     x=[1 1 1 1];
        plot(abs(fft(x)));

Explain this plot. 

     y=[1 1 1 1 zeros(1,N-4)];
     plot(abs(fft(y)));

How does this plot relate to the plots in parts (a) and (b). 

     y=[1 1 1 1 zeros(1,4)];
     plot(abs(fft(y)));
 
     How does this plot relate to the plots in parts (a), (b), and (c).

3. Zero padding in the frequency domain

Zero Padding in Frequency: compute the discrete Fourier transform, Y[n]=fft([1 1 1 1 zeros(1,5)]), and zero pad this signal, Y[n], by inserting zeros in the fractional frequency center (the centre of Y[n]). Then convert the signal back to the time domain, using the command yz=ifft(Y); where yz is the version that's been zero padded in frequency. Is yz real? If not, is it close to being real? Should it be real?

4. Resampling of signals by zero padding (e.g. resizing images)

Load and display a matrix (elements of which happen to be pixels from a picture taken by the autofocus camera): 

 
In Octave*,
a=imread("v145.jpg","jpeg");
image(a/4,1);
 
The original color image has been converted to grayscale image already.
*Note, if you have any problem in using imread in octave, try deleting line 52: save_empty_list_elements_ok= empty_list_elements_ok; and line 336 empty_list_elements_ok= save_empty_list_elements_ok; in imread.m provided in Lab2 folder.
 
In Matlab, notice there is imread.m in Lab2 folder, you must either rename it or delete it.
 
a=imread(‘v145.jpg’,’jpeg’);          
 
image(a);axis image;

 

Note that there are 200 such matrices in this directory. These are differently focused pictures of the same subject matter. There are 200 values of camera focus, and matrix number 145, just-so-happens to be one in which the C++ ("C plus plus") can is in focus.

Note the size of the matrix, using the "size" command: 

 size(a)
 
In Octave,
ans =
  240  320
 
In matlab, 
ans =   
  240  320   3

 

Note also that the first axis points down, and the second axis points to the right, because this is how matrices are indexed (e.g. the first index is the row, and the second index is the column). Thus we refer to the matrix as having a size "240 by 320".

There are some zeros around the outside of the matrix, which show up as black areas in the matrix display. This is because the sensor array in the camera does not extend all the way to the edges of the signal reception area. We wish to crop these away: 

In Octave,
a=a(1:235,10:310);
 
In Matlab,
a=a(1:235,10:310,:);

 

Now display "a" again, and note how the black areas are gone.

If you are using Matlab, now convert “a” to grayscale image.

a=rgb2gray(a);
imagesc(a);
colormap (gray);

You can see the gray scale image.

Now compute the Fourier transform of "a": 

In octave,
A=fft(a);
 
In Matlab,
A=fft(double(a));

 

Note that this command computes the Fourier transform of each COLUMN of "a"

Matlab (Octave) commands typically operate down columns of matrices. For example, if you do something silly like try to plot a matrix: 

plot(a)

you will see 301 plots, one for each column of "a".

Let's try something equally silly in the frequency domain: 

 plot(abs(A))

Now you could also try to display the resulting matrix: 

In octave
image(abs(A));
 
In Matlab,
image(abs(A)); colormap gray;

 

which saturates the range of the display, but even if you scale it: 

 
imagesc(abs(A));

 

the peak is so high you won't see too much other than the peak.

In electrical engineering we often use decibels (log units), so try: 

 imagesc(abs((log(A+1))));
 
Note that the base of the logarithm (e.g. log is base e by default, but decibels are defined on 10log10) doesn't matter because we're scaling the result. Thus we see simply the fact that we're on a log scale, irrespective of the base of the logarithm. The log scale tends to make low values of the function quite visible.

Now if you wanted to see just one column, such as the center column, (original signal on a plot together with its Fourier transform magnitude), plotted out, you would try: 

Octave,
subplot(2,1,1); xlabel('TIME'); ylabel('amplitude')
plot(a(:,151))
subplot(2,1,2); xlabel('FREQUENCY'); % notice how ylabel persists from before
plot(abs(log(A(:,151)+1)))
 
Matlab,
subplot(2,1,1); 
plot(a(:,151));
xlabel('TIME'); ylabel('amplitude');
subplot(2,1,2);  
plot(abs(log(A(:,151)+1)));
xlabel('FREQUENCY');ylabel('amplitude');

Now insert some zeros at the point in the array where the fractional frequency is exactly one half. This happens at 118.5, because the first row (row 1) is zero frequency (recall we're using FORTRAN array indexing, e.g. starting at index 1), and thus there are 117 rows 2:118, and then another 117 rows 119:235. Thus if you put some zeros into the matrix A, e.g.:

 z=zeros(100,301);
Az=[A(1:118,:); z; A(119:235,:)];

Now compute the inverse Fourier transform: 

 az=ifft(Az);

Is it real? Should it be real? How real? Explain.

 

In Octave,

 imagesc(real(az),1)
imagesc(imag(az),1)

 

In Matlab,

 imagesc(real(az)); colormap gray; axis image;
figure; imagesc(image(az)); colormap gray; axis image;

In the above example, you have learned how "optimally" to resize an image, using what amounts to a harmonic ("ideal") interpolation.

Now explain how you would change the width of the image, as well as the height, independently. Hint: to get the transpose, At, of a matrix, A, use the ".'" operation, e.g.: 

 At=A.';

5. Frequency content, and smoothness of functions

Each of the 200 matrices (v000.jpg through v199.jpg) in this directory is a differently focused picture of the same subject matter. Display some or all of the matrices (either in Matlab or Octave, or using your favorite image viewer such as xli, xloadimage, or the like). You will find that the low numbered images are those focused near to the camera, and the high numbered images are those focused far away.

Try plotting a row across the matrix (the grayscale image after conversion), selected near the center of the matrix, for some of the 200 matrices: 

In Octave,
a=imread("v000.jpg","jpeg"); a120=a(120,:); plot(a120)
 
In Matlab,
a=imread(‘v000.jpg’,’jpeg’); a=rgb2gray(a);a120=a(120,:); plot(a120)

 

Since the above is a one-liner, you can do this very easily by just using the up-arrow, or the previous command function (CONTROL P) of the Octave shell to edit the command line. Note that the Octave shell command line editing commands are identical to EMACS commands, e.g. CONTROL A goes to beginning of line, CONTROL F goes forward, CONTROL B goes backward, etc.

Comment on the smoothness of the function, in the neighbourhood of the C++ can, for some or all of the 200 matrices.

What can you say about the relationship between smoothness of the functions in the neighbourhood of the can, and how in-focus the can is?

Now try plotting the Fourier transforms of this same row of the matrix for each of various of the matrices: 

In Octave,
a=imread("v000jpg","jpeg"); plot(log(abs(fft(a(120,:)))+1))
a=imread("v100.jpg","jpeg"); plot(log(abs(fft(a(120,:)))+1))
a=imread("v145.jpg","jpeg"); plot(log(abs(fft(a(120,:)))+1))
 
In Matlab,
a=imread(‘v000.jpg’,’jpeg’); a=rgb2gray(a);plot(log(abs(fft(double(a(120,:))))+1))
a=imread(‘v100.jpg’,’jpeg’); a=rgb2gray(a);plot(log(abs(fft(double(a(120,:))))+1))
a=imread(‘v145.jpg’,’jpeg’); a=rgb2gray(a);plot(log(abs(fft(double(a(120,:))))+1))

 

What can you say about the relationship between how in-focus the camera is, and the magnitude of the Fourier transform of this row of the matrix?

6. Image Processing:

In the above, we varied the high frequency content of the images by bringing them in and out of focus. Thus we were able to lowpass filter images (e.g. blur them) by making the camera lens go out-of-focus. In a sense, therefore, the camera lens is an example of a continuous analog lowpass filter. Now we will examine the effects of lowpass filters implemented by digital signal processing:

A camera simulator called "takepicture.m" written in Octave is provided. If you are running it in Matlab, just change all " to '. Try 

 a=takepicture(100);

 

This will return a picture for the focus set to 100. This camera simulator is very simple. To see how it works: 

 
type takepicture

 

and you can see that all it does is load the appropriately numbered picture.

For this part of the lab you must submit a commented octave or matlab ".m" file for each question and a description of how the script works as well as any plots and images. These scripts must be able to run when they are placed into the same directory as the sample images.