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?