Tillämpad matematik TMA 990 KB2 - 2002/2003

Kölista

Datorstudio 2: Att använda diskreta Fouriertransformer

Inledning.

Vi ska i denna övning bekanta oss med hur vi kan använda Matlab för att utföra spektrala analyser, både för att studera frekvensspektrum för icke-periodiska funktioner samt för att identifiera periodiska förlopp i t.ex. mätserier.

Diskreta tidsberoende förlopp.

När vi använder Matlab (eller andra datorverktyg som Excel etc.) dyker det upp en del begränsningar jämfört med teorin. Framför allt handlar det om att vi inte kan hantera oändligheten (som dyker upp i summationen hos serierna och i integrationsgränserna i transformintegralerna) och att vi endast känner funktionsvärden för en funktion vid vissa tidpunkter. En bra bild är att tänka sig en försökssituation där vi läser av och bokför mätvärden med jämna mellanrum. En annan bild kan vara ljudsignalen på en CD som är diskret; man registrerar i det fallet ljudet 44100 gånger per sekund.

Vi börjar nu med att introducera några begrepp som är viktiga att ha i åtanke när vi jobbar med dessa typer av förlopp. Vi talar om en samplad funktion när vi regelbundet bokför funktionsvärdena. Tiden mellan två avläsningar kallas samplingstiden, Ts, och vi kan också prata om samplingsfrekvensen, fs.

När man bara har värden för vissa diskreta tidpunkter kommer detta att påverka vilka funktioner vi kan representera. Skapa i Matlab en vektor med 11 punkter mellan 0 och 2*pi pch plotta sedan cos(t), dvs


>> t = 0: 2*pi/10 :2*pi;
>> plot(t,cos(t))


Resultatet blir som väntat en cosinusfunktion med perioden 2*pi (om än något kantig). Prova nu att plotta cos(5*t). Även denna gång ser funktionen ut som väntat med svängningar mellan -1 och 1.

Prova nu att plotta cos(10*t). Vad händer? Jo, funktionen blir konstant 1 och svänger inte alls som vi är vana vid. Varför blir det så här? Har Matlab räknat fel? Nej, det som händer är att med det avstånd vi har mellan värdena i t hinner den snabba cosinus-svängningen svänga en hel period och komma tillbaka till 1. Vi ser alltså ingen skillnad på den snabba svängningen och den konstanta funktionen!

Vi kan alltså inte representera hur snabbt svängande funktioner som helst. För att vi ska kunna se svängningsförloppet måste vi ha minst två tidpunkter per period. I detta fallet är Ts = 2*pi/10. Periodtiden för den snabbast svängande funktion vi kan representera blir då T = 2*Ts, i detta fallet T = 2*pi/5. Frekvens blir fn = 1/T, denna gränsfrekvens kallas Nyquist-frekvensen, fn.

Vad kan vi säga om funktioner som svänger snabbare än Nyquist-frekvensen? Vi såg att cos(10*t) såg ut som en konstant funktion cos(0*t). Plotta nu cos(9*t). Denna funktion blir identisk med cos(1*t)! Fortsätter vi ser vi att cos(8*t) blir identisk med cos(2*t) osv. Det som händer är att funktioner med hög frekvens f avbildas som funktioner med frekvens "speglad" i Nyquist-frekvensen fn, dvs funktioner med f = fn + ( f - fn ) avbildas med frekvensen fn - (f - fn). Detta kallas alias-effekt.

Anmärkning:
Om det är så att frekvens för den samplade funktionen är högre än dubbla Nyquist-frekvens, dvs den "speglade" frekvens fn - (f - fn) blir negativ, så speglas frekvens en gång till i noll-frekvensen, osv.

Nyquist-frekvensen och alias-effekten måste man alltid ha i åtanke när man jobbar med diskreta tidsförlopp.

Fouriertransformen i Matlab.

Det kommando i Matlab som utför en diskret Fouriertransform heter fft (för "Fast Fourier Transform"). Vi kan tänka oss att det som görs är att numeriskt räkna ut Fourierintegralen, bara det att det är implementerat på ett annat sätt (som vi inte kommer att gå in på just nu). Den inversa transformen får man med ifft.

Transformen returneras med värdet för nollfrekvensen först, sedan följer värden för ökande frekvens upp till Nyquist-frekvensen, därefter kommer värden för negativa Nyquist-frekvensen upp till noll. Antalet punkter är samma som för vektorn som man ger som indata. Värdena vi får tillbaka är i allmänhet komplexa.

I allmänhet vill dock ha nollfrekvensen i mitten av datamängden (för att t.ex. plotta ett amplitudspektrum) och detta kan man få mha kommandot fftshift (med motsvarande "invers" ifftshift).

För att kunna använda den beräknade transformen måste vi också skapa en vektor med frekvensvärden motsvarande de returnerade transformvärdena. Detta gör vi genom att skapa en vektor från -fn till fn och med samma antal punkter som indata.

Uppgift 1: Plotta amplitudspektrum

Vi ska nu jämföra resultatet från fft med analytiskt uträknade Fouriertransformer för några funktioner. Vi skapar därför först en symmetrisk tidsvektor t från t=-T till t=T och med samplingstiden Ts, välj till att börja med T=10 och Ts=0.1. Skapa också en frekvensvektor freq enligt föregående stycke.

Vi börjar med att titta på funktionen f(t) = exp(-t)*theta(t) där theta(t) är Heavisides stegfunktion som är 1 för tider större än noll och 0 för negativa tider. Enligt tipsen i studioövning 1 kunde vi skapa denna funktion genom att skriva


>> f = exp(-t).*(t>=0);

Denna funktion finns behandlad i Exempel 1, p. 4:04, i boken. Fouriertransformen ges av 1/(1+j*w). Notera här att w indikerar vinkelfrekvensen som är 2*pi högre än den vanliga frekvensen.

För att räkna ut Fouriertransformen mha fft och sedan plotta amplituden skriver vi i Matlab


>> Ff = fftshift(fft(f));
>> plot(freq,Ts*abs(Ff))

Faktorn Ts kommer in för att fft inte normerar transformen på samma sätt som vi gör analytiskt. Rita i samma figur också upp absolutbeloppet av den analytiskt uträknade transformen. Tänk på att ta hänsyn till skillnaden mellan frekvens och vinkelfrekvens. Den imaginära enheten i Matlab kallas för i (förutsatt att man inte själv har döpt någon variabel till i). Vi kan skapa en imaginär enhet med namnet j genom att skriva

>> j = sqrt(-1)

Som ni ser stämmer den numeriskt framräknade Fouriertransformen ganska väl med den analytiska. Felet beror på två saker, dels "integrerar" vi inte över oändligheten utan bara från -T till T, dels ser vi alias-effekten som gör att högre frekvenser än Nyquist-frekvens speglas tillbaka. Vi kan förbättra noggrannheten genom att antingen öka T eller minska Ts. Prova!

Testa med några andra funktioner och deras transformer, på sidan 4:61 i boken finns en tabell där ni kan hämta exempel om ni inte hunnit räkna något själva.

Uppgift 2: Transform -> serie

Den här uppgiften går ut på att studera vad som händer med spektrum för en funktion om den blir mer och mer periodisk. Matlabfunktionen fyrpuls.m returnerar ett tåg av fyrkantspulser och tar tiden och antalet pulser i tåget som inargument. För att se hur funktionen ska användas så ta din tidsvektor från första uppgiften och skriv


>> f = fyrpuls(t,1);
>> plot(t,f)
>> f = fyrpuls(t,2);
>> plot(t,f)
>> f = fyrpuls(t,4);
>> plot(t,f)

Börja nu med endast en puls och plotta amplitudspektrum för funktionen på samma sätt som i uppgift 1. Figuren bör likna Figur 6 på sidan 4:06 i boken, fast vi plottar beloppet av funktion.

Upprepa med fler och fler pulser (eventuellt kan du behöva öka T), vad händer?

Varför?

Jämför med amplitudspektrum för serieutvecklingen av fyrkantsvågen.

Uppgift 3: Spektralanalys av mätdata

Här handlar det om att givet en simulerad mätserie identifiera frekvenstoppar mha fft. Mätserierna består av några funktioner med en bestämd frekvens samt lite pålagt brus.

Börja med att spara "Mätserie 1" i ett bibliotek hos dig och läs in filen i matlab (använd kommandot load). Identifiera eventuella frekvenstoppar mha av fft-kommandot enligt metodiken ovan.

För att kontrollera resultatet genomför vi även en mätning med en mindre samplingstid. Ladda hem "Mätserie 2" och genomför en ny spektralanalys. Får du samma resultat som för "Mätserie 1"? Vad beror skillnaden på?

Vad kan man tänka sig händer med nogrannheten när vi ökar T respektive minskar Ts? "Mätserie 3" innehåller data då vi ytterligare minskat samplingstiden och för "Mätserie 4" är istället totala tiden förlängd. Utför en analys även av dessa serier, men försök att gissa dig till vad du kommer att få för förändringar i spektralanalysen innan du gör något.

Fick du ut vad du förväntat dig? Om inte, kan du komma på varför analysen gav ett annorlunda resultat än du trodde?

Mätserie 1, insamlad under 2 sekunder med samplingstid 2.5 millisekunder.
Mätserie 2, insamlad under 2 sekunder med samplingstid 1 millisekunder.
Mätserie 3, insamlad under 2 sekunder med samplingstid 50 mikrosekunder.
Mätserie 4, insamlad under 20 sekunder med samplingstid 1 millisekunder.

Tillbaka till kurshemsidan.


Editor: Rickard Bergström
Last modified: 2002-11-11