Time-Frequency Analysis FFT with MATLAB Philippe B. Laval KSU - - PowerPoint PPT Presentation

Time-Frequency Analysis

FFT with MATLAB Philippe B. Laval

KSU

Fall 2015

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 1 / 32

Introduction

We look at MATLAB’s implementation of the Fourier transform. To better understand it, we quickly review the relationship between frequency and period.

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 2 / 32

MATLAB and the FFT

MATLAB implements the Fourier transform with the following functions: fft, ifft, fftshift, ifftshift, fft2, ifft2. We describe them briefly and them illustrate them with examples.

1

• fft. This is the one-dimensional Fourier transform. Assuming a signal

is saved as an array in the variable X, then fft(X) return the Fourier transform of X as a vector of the same size as X. Recall that the values returned by fft(X) are frequencies.

2

• ifft. This is the one-dimensional inverse Fourier transform. Given a

vector or frequencies Y , ifft(Y) will return the signal X corresponding to these frequencies. In theory, if Y=fft(X) then ifft(Y)=X.

3

fft2 and ifft2 are the two dimensional versions of fft and ifft.

4

fftshift rearranges the output of fft, fft2 so the zero-frequency component of the array is at the center. This can be useful to visualize the result of a Fourier transform.

5

ifftshift is the inverse of fftshift.

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 3 / 32

MATLAB and the FFT

We give more explanations about how fft works. Using it is not completely straightforward. As we saw above, the DFT is evaluated for the fundamental frequency of a signal and its multiples. Let us see what that really is. Recall that the fundamental period of a function f is the smallest positive real number T0 such that f (t + T0) = f (t). For example, if a function is periodic over an interval [a, b] then its period will be b − a (the length of the interval over which the function repeats itself. If this is the smallest interval over which the function repeats itself, then b − a is the fundamental period, which we often call the period. The frequency of a signal is the number of periods per second. So, if T is the period of a signal, then its frequency f is: f = 1 T , it is expressed in Hz (Hertz). So, the fundamental frequency f0 is f0 = 1 T0 where T0 is the fundamental period.

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 4 / 32

MATLAB and the FFT

Sometimes, when talking about frequency, we also mean the angular frequency: ω defined by ω = 2πf where f is the frequency. So, the fundamental angular frequency, ω0 is defined to be ω0 = 2πf0 = 2π T0 . It is expressed in rad/s. In what follows, when we say frequency, we really mean angular frequency. From the above, it follows that the fundamental frequency of a periodic function over an interval [a, b] is ω = 2π b − arad/s (or 1 b − aHz). Hence, the multiples, denoted ωn, are ωn = 2πn b − a. Over an interval [−L, L] as it is the case in the examples below, we see that ωn = 2πn 2L . The integer values n are called the wavenumbers.

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 5 / 32

MATLAB and the FFT

To compute the DFT, we sample a signal into a finite number of points, we call N. The DFT, in turn, gives the contribution of the first N multiples of the fundamental frequency to the signal. These fundamental frequencies will also be over an interval centered at 0. So, to summarize, the wavenumbers n will run from −N 2 to N 2 − 1. However, the output of MATLAB’s fft corresponds to wavenumbers

• rdered as
• 0, 1, 2, ..., N

2 , −N 2 + 1, ..., −2, −1

• . This should explain

line 8 of the next two examples. For example, in the next two examples, with numpt = 29 = 512, L = 30, then x (the time) will run from −30 to 30, the wavenumber n will run from −255 to 256 hence the frequencies 2πn 60 = πn 30 will run from (−255) (π) 30 = −26. 704 to (256) (π) 30 = 26. 808.

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 6 / 32

MATLAB and the FFT

For example, in the next two examples, with numpt = 29 = 512, L = 5, then x (the time) will run from −5 to 5, the wavenumber n will run from −256 to 255 hence the frequencies 2πn 10 = πn 5 will run from (−255) (π) 5 = −160. 22 to (256) (π) 5 = 160. 85. A last remark. As noted above, the FFT assumes that the function is periodic over the interval of study. If it is not, the FFT will contain an error at the extremities of the interval of study. This can be seen when starting with a signal, using fft on it, then recovering it with ifft.

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 7 / 32

Example 1

This first example simply illustrates how the Fourier transform find the various frequencies which make up a signal. See the code for this example.

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 8 / 32

Example 1

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 9 / 32

Example 2

This next example defines a signal, adds noise to it, takes the Fourier transform of the noisy signal, removes the frequencies with little significance, recovers the signals from these altered frequencies and compares the two signals. See the code for this example.

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 10 / 32

Example 2

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 11 / 32

Example 2

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 12 / 32

2D FFT

The 1-dimensional DFT and FFT can be generalized to higher dimensions. We will focus on dimension 2 here, and apply the FFT to images. In this section, we just do a very brief overview of the main features

• f the 2-dimensional DFT and FFT with MATLAB. We could spend a

whole semester studying the various techniques which can be developed with the FFT to clean images, enhance images, and many

• ther activities on images.

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 13 / 32

2D FFT - Images Review

Gray images are saved as a 2-dimensional array, that is a matrix. If the image has m × n pixels, then the matrix will be m × n. Each entry in the matrix corresponds to a pixel. The value in that entry is the color (gray) intensity. For 8-bit images, it is an integer between 0 and 255. Color images are saved as a 3-dimensional array. The third dimension indicates how many base colors are used. For example, an m × n image saved in the RGB color scheme will be saved as an m × n × 3 matrix which is in fact 3 m × n matrices, one for each color.

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 14 / 32

2D FFT

An image is 2-dimensional signal. It can be also represented as a sum of sine and cosine waves. But we have sine and cosine waves in each dimension. Assuming we have an image f (x, y) of size M × N, it can be represented in the frequency domain by its two-dimensional DFT which is Df (m, n) =

M−1

• k=0

N−1

• l=0

f (k, l) e−i2π( mk

M + nl N )

(1) You will note that the exponential is in fact a product of two exponentials, e−i2π mk

M e−i2π nl N , one sine and cosine wave for each

dimension. As before, Df (m, n) represents the contribution the sine and cosine waves of frequency e−i2π( mk

M + nl N ) make to the image.

Since for images frequency represents color, we can determine which color (or combination of colors) each pixel has.

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 15 / 32

2D FFT

The inverse Fourier transform is f (k, l) = 1 MN

M−1

• m=0

N−1

• n=0

Df (m, n) ei2π( mk

M + nl N )

It allows to recover the image from the frequencies which make it.

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 16 / 32

2D FFT and MATLAB

Many image manipulations happen in the frequency domain, so the outline to perform these manipulations is as follows: Start with an original image (a signal in the time or spatial domain). Take its Fourier transform, fft2, we are now in the frequency domain. Manipulate the image in the frequency domain. This is often referred to as ’apply a filter’. There are filters for pretty much everything one wants to do such as enhancement, blurring, removing noise, detect edges, ... Take the inverse Fourier transform, ifft2, of the modified image in the previous step. We now have our modified image back in spatial domain.

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 17 / 32

2D FFT and MATLAB

In the following example, we will: load an image, add noise to it, take the Fourier transform of the noisy image

• nly keep the frequencies with high contribution, thus eliminating

most of the noise. take the inverse Fourier transform of the filtered image, this will give us a cleaned image.

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 18 / 32

2D FFT and MATLAB

We will discuss each section of the code. In the process, we will learn the following MATLAB functions (we assume our images are m × n pixels):

• imread. Reads an image from its file. The format is

img=imread(filename); where filename is a string containing the filename of the image, img is the matrix containing the image. Note that if the image is grayscale, img will be an m × n matrix, if the image is a color image (RGB), img will be an m × n × 3 matrix. Both matrices will have integer entries. We can convert it to a grayscale matrix with the next command.

• rgb2gray. Converts a color image to a grayscale image. The format is

newimg=rgb2gray(img); where img is a color image to convert and newimg is the converted image, an m × n matrix with integer entries.

• double. Converts integers to doubles (real numbers). The format is

mat2=double(mat1); where mat1 is a matrix with integer entries, mat2 is a matrix with real entries. Real entries are needed for the FFT.

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 19 / 32

2D FFT and MATLAB

• imshow. Displays an image. The format is imshow(img,[ ]); where

img is a matrix containing the image.

• fft2. Computes the 2D Fourier transform. The format is

fimg=fft2(img) where img is a matrix containing the image and fimg is a matrix containing the Fourier transform.

• fftshift. Shift the zero-frequency component to the center of the
• spectrum. This makes visualizing the Fourier Transform easier. The

format is img1=fftshift(img) where img is a matrix and img1 is the shifted matrix.

• ifft2. Computes the inverse of the 2D Fourier transform. The format

is img=fft2(fimg) where img is a matrix containing the image and fimg is a matrix containing the Fourier transform. The next slides show the MATLAB code for this example as well as the resulting images.

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 20 / 32

2D FFT and MATLAB

%Load the image A=imread(’head.png’); % Convert to grayscale Abw=rgb2gray(A); % Convert to double Abw=double(Abw); % Display the original image (see figure 1) figure(1),imshow(Abw, [ ]);

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 21 / 32

2D FFT and MATLAB

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 22 / 32

2D FFT and MATLAB

% Adding noise % Define the noise a=0; % lower bound of the noise value b=100; % upper bound of the noise value noise=a + (b-a).*rand(nx,ny); % define the noise Abwn=Abw+double(noise); %add noise to the image % Display the image with noise figure(2),imshow(Abwn, [ ]); Note: This code is actually not executed on this image since it already had noise. It is just provided to illustrate how to add noise.

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 23 / 32

2D FFT and MATLAB

% Take the Fourier Transform fftA=fft2(Abwn); % Display the Fourier transform (see figure 2) figure(3),imshow(abs(fftA), [ ]);

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 24 / 32

2D FFT and MATLAB

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 25 / 32

2D FFT and MATLAB

%Try again. In the FFT, high peaks can be so high that they hide the %details. We can reduce the scale with the log function (see figure 3). figure(4),imshow(log(1+abs(fftA)), [ ]);

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 26 / 32

2D FFT and MATLAB

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 27 / 32

2D FFT and MATLAB

%Try again. We use fftshift to put the coefficients with 0 frequency at the center (see figure 4). figure(5),imshow(log(1+abs(fftshift(fftA))), [ ]);

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 28 / 32

2D FFT and MATLAB

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 29 / 32

2D FFT and MATLAB

%remove the noise with a high pass filter Threshold=median(2*median(abs(fftA))); % define what to keep cffA=(abs(fftA)>Threshold).*fftA; % Keep frequencies with contribution > Threshold % Recover the cleaned image using the inverse Fourier transform ifftA=ifft2(cffA); % recover the ’cleaned’ image % Display the cleaned image (see figure ) figure(6),imshow(ifftA, [ ]);

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 30 / 32

2D FFT and MATLAB

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 31 / 32

Exercises

See the problems at the end of the notes on the FFT and MATLAB.

Philippe B. Laval (KSU) MATLAB FFT Fall 2015 32 / 32