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.

Your submission should be printed out and put in the drop box on the 4th floor of galbraith building, outside room 441.

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):

a) find the FT [F(w)] of f(t) = exp(-(t^2)/Z) where Z = 4 hint multiply f(t) by exp(-w^2)*exp(w^2) = 1 Sketch a graph of this function.
try varying the function f(t) by varying Z (try values ranging from Z=0.1 to 100). what do you observe in the time domain? What do you observe in the frequency domain for each of the corresponding time domain functions?

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.

a) Consider a time domain signal, x=[1 1 1 1]; Construct a plot of the magnitude of the Discrete Time Fourier Series (DTFS) of this signal, x. This can be done using the DFT (FFT) function in Octave:
```
     x=[1 1 1 1];
     plot(abs(fft(x)));
     
```
Explain this plot.
b) Consider a discrete time signal: x[n]=1, 0<=n<=3; x[n]=0, otherwise. Determine the Discrete Time Fourier Transform (DTFT) for this signal. Sketch a graph of this function.
c) Zero padding: Consider a discrete time periodic signal, constructed by periodizing the above signal: y[n]=sum(x[n-lN]), for integer l, and integer period N. Suppose N=1000. Construct a plot of the magnitude of the Discrete Time Fourier Series (DTFS). This can be done using the DFT (FFT) function in Octave:
```
     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).
d) Repeat part (c) above, with N=8 (e.g. padding with 4 zeros instead of 996 zeros):
```
     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).
e) In (c) and (d) above, you were padding the signal, in the time domain, with zeros. Explain the effect this time-domain zero padding has in the frequency domain.

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. 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?

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):

a=imread("v145.jpg","jpeg");
image(a/4,1);

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)
ans =
  240  320

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:

a=a(1:235,10:310);

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

Now compute the Fourier transform of "a":

A=fft(a);

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

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:

image(abs(A),1)

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

imagesc(abs(A),1)

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))),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:

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)))

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.

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

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.';

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 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, selected near the center of the matrix, for some of the 200 matrices:

a=imread("v000.jpg","jpeg"); 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:

a=imread("v000.jpg","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))

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?

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) Look at one row of image, lowpass filter it (using a lowpass filter of your choice, such as a rectangular filter, or a Gaussian filter above). Hint: you can multiply by a window filter in the frequency domain, or you can equivalently use the "conv" command in Octave, e.g. conv(f,[1 1 1 1 1 1 1 1]);. Plot the original together with the lowpass filtered version. Describe the effect of lowpass filtering on the smoothness of the function.
b) What is the effect you observe in the frequency domain
c) Given two functions, one being a lowpass filtered version of the other, suggest a way of determining which of these is the one that has been lowpass filtered. Hint: compare the absolute values of the FFTs of the two signals.
d) Apply a measure of the extent to which a signal has been lowpass filtered, based on part (c). Could you use var(fft(X)); var(fftshift(fft(X))); mean(fft(X)); mean(fftshift(fft(X))); mean(fft(X(1:N/2)))? Why is it that none of these really work? Hint: take a look at "extent.m" and either use it, or devise your own measure of extent (e.g. based on first or second moments of the range of Discrete Fourier Transform indices).

Here's the fun part of the lab. The task is to develop an autofocus algorithim for the camera pictured below:

In order to allow people to do this lab at home, the sequence of test images is provided in this directory. The camera has 200 focus positions, and there is one image per focus position. The lower the number of the image, the closer an object has to be to the camera, to be in focus.

A camera simulator called "takepicture.m" is provided. 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 ".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. To judge the algoroithms the TAs will run your script against the sample images as well as in an alternate scene.

a) write a script to focus the camera on the C Plus Plus can. Hint: consider only the center portion of the image, or use a function, such as a sampled version of the Gaussian, or a simple rectangular window, to make your script center-weighted.
b) write a script that focuses on the PCI "interface" poster board in the upper right hand area of the image. Hint: consider only that portion of the image.
c) optional bonus question: for an additional 5% added to your final grade (e.g. approximately the grades of one entire lab), write an algorithm to automatically search through all 200 pictures and combine the sharpest regions of each of these, to obtain a single picture in which everything in the picture is in focus at the same time. For this question, hand in the final picture as well as the code. Your code must work on any arbitrary image sequence from any arbitrary scene, and not just for the specific subject matter depicted in the image sequence. It must also run in a reasonable length of time e.g. certainly not more than an hour on a standard pocket size EyeTap wearable computer (1 GHz P3). What is the ECE431 exam question pool server password written at the top of the ECE431 question server computer keyboard?