1 / 24

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.

floria
Download Presentation

Numerical Methods

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript


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

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

  3. 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…

  4. A labelled plot…

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

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

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

  8. 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)]

  9. 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’

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

  11. 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,…)

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

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

  14. Signal plus noise

  15. 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)')

  16. Frequency content (incl. noise)

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

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

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

  20. 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 =

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

  22. Comparing them Hamming Blackman Hann

  23. 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]

  24. 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…

More Related