Fourier methods in signal processing

Kristian Sandberg
Dept. of Applied Mathematics
University of Colorado at Boulder

Instructions

This lab is due Monday March 6 in lecture. You are encouraged to discuss the lab with each other, but your solutions must be independently written using you own words and formulations.

Goal

The goals with this lab are:

1.
Learn how to analyze a one-dimensional signal by using the discrete Fourier transform (DFT).
2.
Develop intuition for how spectrum behaves for different types of signals.
3.
Apply Fourier analysis to de-noising and compressing signals.

4.
Compare different methods to perform a convolution.
5.
Learn how to use Fourier methods for image processing (compression and low-pass filtering).

One dimensional signals

Background

The first examples of a linear space (or a vector space) we usually encounter when we study mathematics, are the Euclidean spaces $\mathbb{R} ^2$ and $\mathbb{R} ^3$. We have a good intuition for these spaces in terms of basis vectors, lengths and angles. However, there are other useful linear spaces. Another important class of examples are function spaces which is an important concept in applied mathematics. The functions in a function space correspond to what we think of as vectors in $\mathbb{R} ^2$ and $\mathbb{R} ^3$. We define scalar multiplication and addition for functions in a function space in a linear manner and we can even define norms on this space which corresponds to our notion of ``length'' in Euclidean geometry. We can find a correspondence to the notion of angles by introducing an inner product in our function space. We can introduce a basis on our space.

In this lab we will study the function space spanned by the trigonometric functions (or functions of the type eikx if we use complex numbers). We will try to represent a signal in terms of these trigonometric functions and then work with the coefficients that our signal is represented by in this basis.

Throughout this lab, the elements in the DFT vector of a signal will be referred to as the spectrum of the signal. In other words, the coefficients which represent the signal in the Fourier basis is the spectrum of the signal. Compare this terminology to the content of the spectral theorem in linear algebra (Theorem 1.109) in Frazier [1].

One dimensional signals

In this lab we will study one dimensional signals (first part) and two dimensional signals (section 4). An image is a two dimensional signal but we will see that once we understand Fourier analysis of one dimensional signals, we will almost automatically understand how we can use Fourier analysis for images. However, one dimensional signals (sometimes referred to as time series) have plenty of applications!

The smoothness of a signal

A signal or function is referred to as being ``smooth'' if the function and ``sufficiently many'' of its derivatives are continuous. The more derivatives that are continuous, the smoother is our signal.

Since Fourier analysis means representing a signal with respect to a basis of periodic functions, our signal should also be periodic for an efficient representation in this basis. This is of course not true for all signals. Since we will study finite signals in this lab, we therefore think of our finite signal to be periodically extended as described in the beginning of Chapter 2.1 in [1]. Once we have periodically extended our signal in this manner, there is a useful ``rule of thumb'' that relates the structure of our signal to its spectrum:

This rule of thumb can of course be stated (and proved) more rigorously. However, for this lab the above statement should be sufficient to explain some of the spectra you will encounter.

The fundamental principle for signal processing using Fourier analysis

A fundamental principle when processing a signal using Fourier analysis is to manipulate the spectrum of a signal rather that manipulating the signal itself. Hence, the usual procedure is to find the DFT of a signal, manipulate the DFT vector by, for example, letting some of its elements equal zero to get rid of unwanted frequencies, and then transform back our signal using the inverse DFT.

Compression

When we compress a signal using Fourier analysis we neglect frequencies that have zero or almost zero magnitude in the spectrum. In many cases this means that we only need a small fraction of the spectrum to represent a signal. Instead of storing the whole signal, we just store the coefficients for the largest frequencies in the spectrum.

Let us take an example. Consider a signal with 512 samples, that is, 512 data points. To store the whole signal requires 512 pieces of information to be stored. If we take the DFT (fft() in Matlab) of our signal, we get a new vector with 512 elements. Usually a great portion of this vector is ``almost'' zero. It may very well be that only 50 of the elements are ``large enough'' to significantly contribute to the signal. Therefore, we store only these 50 elements of the DFT vector. Once we need our signal, we just take an inverse transform of our DFT vector to reconstruct a good approximation of our original signal (ifft() in Matlab). For a definition of how to measure compression ratio, please see The Singular Value Decomposition of an Image.

The l2-norm of a signal

When measuring the error in signal processing, we usually use the l2-norm. Let our (finite) signal be given by $z=(z(0), z(1),\ldots, z(N-1))$ and a processed version of the original signal be given by $\tilde{z}=(\tilde{z}(0), \tilde{z}(1),\ldots, \tilde{z}(N-1))$. We define the relative l2-error as

\begin{displaymath}E_{rel}=\sqrt{\frac{\sum_{k=0}^{N-1}\vert z(k)-\tilde{z} (k)\vert^2}{\sum_{k=0}^{N-1}\vert z(k)\vert^2}} \ . \end{displaymath} (1)

Are Fourier methods still competitive in signal processing?

Even though wavelets outperforms Fourier methods in many situations, there are still some signals which are better represented using Fourier methods. Also, even if we decide to use wavelets for an application, understanding wavelet analysis often requires solid knowledge in Fourier analysis.

Some useful Matlab commands

The following table gives a few commands that will be useful for this lab.

Operation: Matlab command
The DFT of a signal (vector) z.  
(non-normalized) fft(z)
The inverse DFT of a vector fz.  
(Normed with factor 1/N, where  
N is the length of the vector) ifft(fz)
Switch position of first and second half  
of vector fz. fftshift(fz)
Start a ``stop watch'' tic
Find the time elapsed since the stop watch started. time=toc
Generate a symmetric Toeplitz matrix where  
the first column is given by the vector c toeplitz(c)
Generate a Toeplitz matrix where  
the first column is given by the vector c and the first row  
is given by the vector r toeplitz(c,r)
The real part of a scalar (or vector) z. real(z)
The imaginary part of a scalar (or vector) z. imag(z)
The magnitude part of a scalar (or vector) z. abs(z)

Note that the last three commands in the table act element-wise on vectors. For example, if z=(-1,3+4i,-2i) then abs(z) will give z=(1,5,2) as output. When asked to find the magnitude of a signal in the exercises below, this means that you should use the command abs(). (As opposed to finding the norm of a vector for which you want to use the command norm().)

Warning! When multiplying a matrix A with a vector x to form Ax, the vector has to be a column vector! You can form the transpose in Matlab with the ' operator.

Displaying the spectrum

When you display the DFT of a signal in Matlab, the default setting is that the low frequencies are displayed at the right and left part of the plot with high frequencies in the center. This is how the spectra in Chapter 2.1 in [1] are shown. However, in many texts, one displays the spectrum with the low frequencies in the center, and the high frequencies at the left and right edge of the plot. Which way you choose is up to you. You can accomplish the latter alternative by using the command fftshift().

Exercises

Exercise 1: An introductory example

For each of the signals $z_1(t)=cos(8\pi t)$ and $z_2(t)=sin(16\pi t)$ plot: Reconstruct the signal from its DFT and verify that you indeed get back the original signal.

Compare the spectra of the two signals. What are the differences, similarities? Do the spectra look like you expected? Discuss why/why not.

Hint. An easy way to generate the first signal in Matlab using 128 samples is to first type grid=0:127; followed by z1=cos(pi*grid/16);. Plot the signal by typing plot(grid,z1).

Exercise 2: Manipulating a signal

For the signal $z(t)=cos(8\pi t)+sin(16\pi t)$ plot: Now manipulate the spectrum of the signal so that the reconstruction based on the manipulated spectrum only contains the sin-component of the signal. Plot the reconstruction.

Exercise 3: The spectrum of z(t)=t2

Plot the signal z(t)=t2 on the interval [0,1]. Plot the spectrum of the signal. Can you see any major qualitative difference between this spectrum and the spectra for the signals in the previous exercises? Can you explain the difference?

Exercise 4: The spectrum of a discontinuous vs. a continuous signal

Down load the two signals square wave and saw tooth. Plot the signals along with their spectra. For each signal, set the last 90% of the elements in the DFT vector equal zero. Plot the spectrum again along with a reconstruction of the processed signal. Which reconstruction is most faithful to the original one? Compute the relative l2-error for both the processed signals and comment on the result. Relate your conclusion to the continuity of the original signals.

Exercise 5: The uncertainty principle

a) Plot the signal z(t)=e-a(t-0.5)2 on the interval [0,1] along with its spectrum for a=1, 100, 10000 and 106.

b) What can you say about the signal and its spectrum as a increases? Can you state a general (qualitative) principle that relates the ``width'' of a signal with the width of its spectrum?

Note: The principle you should find is related to the famous ``Heisenberg's uncertainty principle'' in quantum mechanics which says that we cannot determine a particle's position and velocity exactly simultaneously. In a certain sense, the momentum of a particle in the quantum mechanical picture is the Fourier transform of the position of the particle.

Exercise 6: Comparison of two signals with identical energy

Consider the functions z1(t)=t and z2(t)=1-2|t-1/2| on the interval [0,1]. By computing the energy of these functions on [0,1] (by evaluating $\sqrt{\int_0^1 z_1(t)^2\ dt}$ and $\sqrt{\int_0^1 z_2(t)^2\ dt}$) one finds that these signals have the same energy. Now plot the spectra of these two functions and compare the decay carefully. Which decays faster? Can you explain why?

Exercise 7: The global property of the DFT

Plot the function

\begin{displaymath}z(t)=\left\{ \begin{array}{ll} cos(8\pi t) & t\in[0,1/2) \\ cos(32\pi t) & t\in[1/2,1] \end{array} \right. \end{displaymath}

on the interval [0,1). Plot the spectrum of the signal

It is very tempting to say that ``this is a superposition of two pure frequencies'' but this is not right! If this were a superposition of two pure frequencies, what would the spectrum look like? Why does the spectrum of our signal differ from a spectrum of two pure frequencies?

Relate this example to Lemma 2.13 in [1] and the following discussion. Take some time and think about what this lemma means and why the spectrum of the signal above does not only consist of two frequencies. This concept will be very important later on when we study wavelets!

Exercise 8: De-noising a signal 1

Down-load the signal z here. Plot the signal and its spectrum. Remove the noise from the signal!

(Hint. Noise can be created by a random process which can be modeled in Matlab using the command rand(). To get an extra hint for how to de-noise, plot a random signal and its spectrum in Matlab by using rand(1,N) where N is the length of the signal.)

Exercise 9: De-noising a signal 2

You are working for a record company that wants to release an old vinyl record on a CD. The problem is that the only analog recording they have left has a scratch at one place which listeners will find very annoying. Your task it to down-load a digitized version of the analog recording and get rid of the scratch by using your knowledge in Fourier analysis. (Do the de-noising in the Fourier domian. In this case, it would in principle be possible to ``clean'' this signal in the space domain, but in more general cases, this would not be practically possible.)

Down-load the signal here.

Exercise 10: Signal compression

Plot the signal $z(t)=e^{-4t}sin(16\pi t)$ on the interval [0,1] and its spectrum. Now compress the signal by letting appropriate elements in the DFT vector equal zero. See how many elements in the DFT vector you need to keep in order for the compressed signal to look roughly the same as the original signal. What is the relative l2-error after the compression? What is the compression rate?

Bonus question: Can you describe a physical situation where this kind of signal occurs?

Exercise 11: Convolution

Down load the following signal z. Define a vector w that samples the function $cos(32\pi x)$ at 1024 points on the interval between 0 and 1. Compute the convolution w*z in the following ways:
1.
In the space domain. One way to do this is to construct the matrix representing the translation invariant operator associated with w and then multiply the resulting matrix with z.
2.
Fast convolution. Transform w*z into Fourier space, perform the convolution and then transform back to the space domain.

Verify that both methods give the same result (hand in a plot of the convolved signal). Now time the both methods and comment on the difference!

Even though there are convolution commands in the Image processing toolbox in Matlab, try to solve this exercise without using these. However, you can of course use the built in fft() commands!

  
Fourier methods applied to image processing

Background

The theory introduced for one dimensional signals above carries over to two dimensional signals with minor changes. Our basis functions now depend on two variables (one in the x-direction and one in the y-direction) and also two frequencies, one in each direction. See Exercises 2.15-2.18 in [1] for more details. The corresponding l2-norm for a two dimensional signal now becomes

\begin{displaymath}\vert\vert A\vert\vert _2 = \sqrt{\sum_{i=1}^M \sum_{j=1}^N \vert a_{ij}\vert^2} \end{displaymath} (2)

where aij are the elements in the $M\times N$ matrix representing the two dimensional signal. It is computed in Matlab using the Frobenius norm.

Some useful Matlab commands

The following table gives some commands that will be useful for this part of the lab.

Operation: Matlab command
The DFT of a 2D signal (matrix) A.  
(non-normalized) fft2(A)
The inverse DFT of an $M\times N$ matrix fA.  
(Normalized with a factor 1/MN.) ifft2(fA)
Switch position of first and third quadrant and  
second and third quadrant of matrix fA. fftshift(fA)

 

The commands real, imag and abs work on matrices elementwise just as in the one dimensional case.

Displaying the spectrum

When you display the DFT of a two dimensional signal in Matlab, the default setting is that the low frequencies are displayed towards the edges of the plot and high frequencies in the center. However, in many situations one displays the spectrum with the low frequencies in the center, and the high frequencies near the edges of the plot. Which way you choose is up to you. You can accomplish the latter alternative by using the command fftshift().

There is a very useful trick to enhance the visualization of the spectrum of a two dimensional signal. If you take the logarithm of the gray-scale, this usually give a better plot of the frequency distribution. In case you want to display the spectrum fA, I recommend to type imshow(log(abs(fA))) for a nice visualization of the frequency distribution. You may also want to use fftshift as described above, but that is more a matter of taste.

Exercises

Exercise 12: The two dimensional Fourier basis

Define the basis function Fm,n for an $M\times N$ matrix by

\begin{displaymath}F_{m,n}(k,l) = \frac{1}{\sqrt{MN}}e^{i2\pi(\frac{mk}{M}+\frac{nl}{N})},\ k=1,2,\ldots,M,\ l=1,2,\ldots,N \end{displaymath} (3)

To get a feeling for what these functions look like, plot the real and imaginary part of F1,1, F2,0, F0,2 and F4,4. Do this by first evaluating these functions on a square grid (128 by 128) and then display the resulting matrix using the command imshow.

Finally, form a sum of two (real) basis functions such that the resulting imshow-plot ``resembles'' a chess board.

You may have to re-scale the matrix elements before you display the matrix using imshow. You can accomplish this using the mat2gray() command.

Exercise 13: Thresholding

To compress an image, set elements in the Fourier space of a signal with a magnitude less than some number $\epsilon$ to zero. (A matrix with many zeros can be stored very efficiently.)

Write a function threshold that takes a matrix and a scalar representing $\epsilon$ as input. The function should loop through all elements in the matrix and put all elements with a magnitude less than epsilon to zero and compute and print (on the screen) the compression ratio (use definition in previous lab and assume that every non-zero element is ``one piece of information'' that needs to be stored). Finally, the function should return the resulting matrix.

Exercise 14: Image compression

Using the function threshold, compress (by performing thresholding in the Fourier domain) the two images pic-home.jpg and basket.jpg you worked with in the previous lab. Discuss the compression in terms of visual quality, compression rate and the l2-error. Is the performance different for the two images and if so why? How does this compression method compare to the SVD-compression for the previous lab?

Exercise 15: Low pass filtering

A ``naive'' way to low pass filter an image (that is, remove high frequencies of the image) is to Fourier transform the image, put all elements representing high frequencies of the Fourier transformed signal equal to zero, and then take an inverse transform of the signal.

Design a low pass filter that removes the highest frequencies in both the x-direction and in the y-direction. Experiment by removing more and more high frequencies and see how this changes the the filtered image.

The Image processing toolbox in Matlab has some built-in filtering commands, but try to use your own commands, with the exception of the built-in fft2() commands.

This filtering is nothing else than a convolution. Could we do this without using the Fourier transform? If so, why is that not a desirable approach?

As mentioned above, this is a ``naive'' way of low pass filtering. It amounts to filter the frequencies with a ``rectangular'' function, and this may not be the best approach in all cases.

Format of your lab report

Write down your comments to all of the exercises above (word processed). Include a representative collection of the code(s) that you used as an appendix. Make sure to include comments in your code explaining what it does. Include any plots or images you need to justify your conclusions.

Questions?

In case you have questions regarding the material in this lab, do not hesitate to contact me at kristian.sandberg@colorado.edu or visit me during my lab hours Mondays 4:30-6 pm and Thursdays 4-5:30 pm in ECCR 143.

Lycka till!

References

[1] Frazier, M.W., An introduction to wavelets through linear algebra, Springer-Verlag, 1999.

About this document ...

Fourier methods in signal processing

This document was generated using the LaTeX2HTML translator Version 98.1 release (February 19th, 1998)

Copyright © 1993, 1994, 1995, 1996, 1997, Nikos Drakos, Computer Based Learning Unit, University of Leeds.

The command line arguments were:
latex2html -split 0 -no_navigation fa1.tex.

The translation was initiated by Kristian Sandberg on 2000-02-17

Kristian Sandberg
2000-02-17