Numerical Methods

Numerical Methods Lecture 9 – Frequency Analysis (using FFT) in Matlab Dr Andy Phillips School of Electrical and Electronic Engineering University of Nottingham

Today’s Topics • More Matlab • Frequency Analysis • What it is & why we do it • Windowing • Frequency Analysis in Matlab

More Matlab 1: Labelling Plots • In the last lab we could have provided plot title and x and y axis labels: %excerpt from square wave script file plot(t,sum); title(‘Creation of a Square Wave’) xlabel(‘t’) ylabel(‘x(t)’) • Like comments, labels make your work clearer to the reader…

A labelled plot…

More Matlab 2: abs • To obtain an absolute value or complex magnitude of a scalar or matrix number or variable use abs • abs(-2)2 • if you’ve defined i=sqrt(-1) then abs(1+i)1.4142 • if you’ve defined x=-3 then abs(x)3

More Matlab 3: functions • Matlab has lots of familiar maths functions: • sin, cos, tan • log, log2, log10, exp • log is natural log (i.e. base e) • cosh, sinh, tanh • factorial • And less familiar functions e.g gamma, sinc

Frequency Analysis • One common, and computationally intensive, method we use is Frequency Analysis • The most common task is performing Fourier and Inverse Fourier Transforms • related to Fourier Series, but don’t need signal to be periodic

Fourier Transforms • You’ll do these in Signal Processing next year, so don’t worry about the theory. For info: • Fourier Transform of h (t) is H (f) where: • Inverse Fourier Transform of H (f) is h (t) where: [j=sqrt(-1)]

What does it mean? • Fourier Analysis • This is a technique which takes a data in the time domain and provides its frequency spectrum • The inverse takes us back again • We have a ‘Fourier Transform pair’

DFT and FFT • In engineering functions often represented by finite sets of discrete values e.g. f(t) below could be represented by the set of discrete points (ti ,fi) below

DFT and FFT • Using such discrete data a discrete Fourier Transform (DFT) can be defined: • Computationally intensive (N2) when calculated straightforwardly • Fast Fourier Transform (FFT) computes DFT using very efficient algorithm • uses results of previous calcs to reduce number of operations • can get to Nlog2N operations, so for N=1024 the FFT is approx. 100 times faster than straight DFT • often N required to be an integer power of two (i.e. 2,4,8,16,32,…)

Now… • Having performed the analysis we can • Determine bandwidth • Remove Phase information • Generally play around with the signal • It can be done in C, but it is a lot of work • Use a standard library or a package like Matlab • Your next lab will involve you using Matlab’s fft function

Simple script file using FFT %simple example of FFT use %setting up the signal t=0:0.001:0.6; x=sin(2*pi*40*t)+sin(2*pi*100*t) %40 Hz and 100 Hz sine waves added y=x+2*randn(size(t)); %noise added %randn(size(t)) gives an array of normally distributed random entries % that is the same size as t plot(t(1:50),y(1:50)) %plots first 50 points i.e. to 0.05 s title('40 Hz plus 100 Hz corrupted by noise') xlabel('time (seconds)') pause(2) %waits 2 seconds before continuing

Signal plus noise

Simple script file using FFT pt2 %the fft and output Y=fft(y,512); %performs the fft. Note 512 elements - power of 2 (2^9) Pyy=Y.*conj(Y)/512; %Y will have complex elements. The power spectrum is real %so do an element by element multiplication by complex conjugate f=1000*(0:256)/512; plot(f,Pyy(1:257)) %we are only interested in half the data points %as the same information is provided in the second half of Pyy title('Frequency content of corrupted signal') xlabel('Frequency (Hz)')

Frequency content (incl. noise)

Before FFT a signal however • We often apply a ‘window’ to the data. • This simply means taking the amount we want from the data stream • ie The window is moved along the data; we perform the FFT on this windowed data

However • We extract the data simply by reading data between a start and finish value • in effect we apply a rectangular window • This however causes problems of discontinuities, as shown in the next slide

The Problem… These discontinuities cause errors in our frequency analysis. To avoid this we use a window that ‘tapers’ at the ends

The Hann & Hamming Windows • The most common solution is to use either Hann or Hamming Windows • Are simply cosine bells, scaled to be the width of the window used (in sample points) • They differ only in the choice of  • Hann ->  = 0.5 Hamming ->  = 0.54  - (1-  )cos(2*pi*i/(n-1)) ; 0<i<n 0 ; otherwise xi =

So.. • The filter coefficients w(i) of a window of length n are computed according to the following formulae • Hamming window • w(i) = 0.54 - 0.46*cos(2*pi*i/(n-1)) • Hann window • w(i) = 0.5 - 0.5*cos(2*pi*i/(n-1))

Comparing them Hamming Blackman Hann

Now.. • Since we often apply this window over and over • It might help to put it in a lookup table so we only have to calculate it once • We then simply multiply our array of data by the corresponding element in the window’s array • i.e. • xnew[i] = x[i] * w[i]

And then • We perform our FFT to get the frequency spectrum • As stated before, we could do this in C, or even manually, but it is simpler to use Matlab • (or some standard library e.g. NAG – Numerical Algorithms Group) • More FFTs and windowing in your lab on Wednesday…

Numerical Methods