Numerical methods
Download
1 / 24

Numerical Methods - PowerPoint PPT Presentation


  • 130 Views
  • Uploaded on

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.

loader
I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.
capcha
Download Presentation

PowerPoint Slideshow about ' Numerical Methods' - floria


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.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.


- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript
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
Today’s Topics

  • More Matlab

  • Frequency Analysis

    • What it is & why we do it

    • Windowing

  • Frequency Analysis in Matlab


More matlab 1 labelling plots
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…



More matlab 2 abs
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
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
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
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
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
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 fft1
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 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



Simple script file using fft pt2
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)')



Before fft a signal however
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
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
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 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
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
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…


ad