DSP Labo 3&4: Fourier
24 januari 2005
Inhoudsopgave 1 2 2.1 2.1.1 2.1.2 2.1.3 2.2 2.2.1 2.2.2 2.2.3 2.3 2.4 3 3.1 3.2 3.3 3.4 3.5 3.6 3.7 4 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 5
Inleiding 3 Analyse 3 Fourierreeks 3 Complex 3 Som van sinussen en cosinussen 3 Verband tussen beide vormen 4 Fourierreeks van enkele bekende functies 5 De driehoeksgolf 5 De zaagtand 6 De vierkantgolf 7 Correlatie 8 Verband FFT-DFS-DFT 8 Experimenten 8 Vraag 1 8 Vraag 2 8 Vraag 3 10 Probleem 1 10 Probleem 2 11 Probleem 3 13 Probleem 4 15 Resultaten 15 Vraag 1 15 Vraag 2 16 Vraag 3 19 Voorbeelden 20 Probleem 1 20 Probleem 2 20 Probleem 3 23 Probleem 4 24 Besluit 27
1
Lijst van figuren 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
De driehoeksgolfvorm en de gekozen periode 5 De zaagtand golfvorm en de gekozen periode 6 De vierkant golfvorm en de gekozen periode 7 De eerst vier FS co¨effici¨enten van vraag 2 16 Derde term fourier benadering van 50 Hz driehoeksgolf 17 Zevende term fourier benadering van 50 Hz driehoeksgolf 17 Twaalfde term fourier benadering van 50 Hz driehoeksgolf 18 Eenentwintigste term fourier benadering van 50 Hz driehoeksgolf 18 De golfvorm van vraag 3 met zijn twee harmonischen 19 De synthese van de zaagtand van probleem 1 21 Analyse van de sinusgolf uit probleem 2 21 Analyse van de vierkantgolf uit probleem 2 22 Synthese van de sinusgolf uit probleem 3 23 Synthese van de vierkantgolf uit probleem 3 24 Golfvormen van probleem 4 25 Analyse van probleem 4 26
2
1
Inleiding
In dit labo bestuderen we de Fourier-analyse en synthese zowel in het continue als het discrete domein. Eerst maken we de Fourier synthese van enkele typische golfvormen met de hand in het continue domein. Daarna verifi¨eren we deze bekomen waarden met Matlab. Aangezien Matlab, zoals in vorige labo’s gezien, discreet werkt zullen we hier dus het verband tussen Fourier in het continue en discrete domein moeten leggen. Matlab berekent Fourier aan de hand van de Fast Fourier Transformatie. Deze FFT is een versnelt algoritme voor de berekening van de discrete tijd Fourier transformatie. We onderzoeken in dit labo dan ook hoe we deze FFT kunnen gebruiken om onze zelfberekende Fourier co¨effici¨enten te verifi¨eren. Deze opgedane kennis omtrent Fourier gebruiken we dan om enkele ideale filterbewerkingen uit te voeren op periodieke signalen. De bekomen spectra gebruiken we dan voor enkele analyses.
2
Analyse
2.1
Fourierreeks
Er zijn twee soorten formules voor de fourier reeks analyse en synthese. Enerzijds de complexe vorm en anderzijds de sinus vorm. 2.1.1 Complex De complexe vorm wordt gegeven door: T 2 1 Z Ck = fT (t) · e−j2πkF0 t dx T T
(1)
−2
fT (t) =
+∞ X
Ck · ej2πkF0 t
(2)
k=−∞
2.1.2 Som van sinussen en cosinussen Deze reeks kunnen we ook uitdrukken als een som van sinussen en cosinussen: fT (t) =
∞ ∞ X A0 X + Ak · cos (2πkF0 t) + Bk · sin (2πkF0 t) 2 k=1 k=1
(3)
3
2.1.3 Verband tussen beide vormen Als we nu in (3) volgende substituties doen: ej2πkF0 t + e−j2πkF0 t cos(2πkF0 t) =
2
(4)
ej2πkF0 t − e−j2πkF0 t sin(2πkF0 t) = 2j
, dan krijgen we:
∞ A0 X ej2πkF0 t + e−j2πkF0 t ej2πkF0 t − e−j2πkF0 t fT (t) = Ak + + Bk 2 2 2j k=1 A0 = C0 ⇓ 2 ÃÃ ! Ã ! ! ∞ X Ak Bk j2πkF0 t Ak Bk −j2πkF0 t = C0 + + − e + e 2 2j 2 2j k=1
Ã
⇓ =
Ã
!
Ã
!
!
Ak Bk −j2πkF0 t Ak Bk j2πkF0 t + − e + e = Ck ej2πkF0 t + C−k e−j2πkF0 t 2 2j 2 2j
∞ X
Ck · ej2πkF0 t
(5)
k=−∞
Daaruit kunnen we dus volgende substituties halen: Ak = 2 · ℜ{Ck } Bk = −2 · ℑ{Ck } q 2 2
Ak Bk + = Ck ⇒ |C | = k 2 2j
Ak + Bk
(6)
2 Bk arg(Ck ) = ϕk = arctan Ak
Als we deze in (3) zouden brengen bekomen we: ∞ X A0 fT (t) = |Ck | · cos (2πkF0 t − ϕk ) +2· 2 k=1
We kunnen daarin dan nog Rk = 2 · |Ck | = fT (t) =
q
(7)
A2k + Bk2 dan bekomen we uiteindelijk:
∞ A0 X + Rk · cos (2πkF0 t − ϕk ) 2 k=1
(8)
4
2.2
Fourierreeks van enkele bekende functies
2.2.1 De driehoeksgolf Gegeven:
Figuur 1: De driehoeksgolfvorm en de gekozen periode
¯ ¯ ¯ 2H ¯ ¯ y(t) = ¯ · t¯¯
T
Gevraagd: A0 ;Ak ;Bk
(9)
Oplossing: A0 T
A0
+2 2 Z = y(t) dt T T −2
2 (oppervlakte1periode) T 2 HT = × =H T 2
=
(10)
Ak T
Ak
=
+2 2 Z y(t) cos(kω0 t) dt T T −2
T
=
+2 2H Z 2 t cos(kω0 t) d × T T T −2
T
=
+2 4H 2 Z t d sin(kω0 t) × T2 kω0 0
5
P.I.
=
= = =
T
4H sin(kω0 t)|T /2 − 0 kπT
+2 Z 0
sin(kω0 t)dt
4H 1 [− cos(kω0 t)]T0 /2 0− kπT kω0 ¶ µ 1 4H [cos(kπ) − cos(0)] kπT kω0 4H 2H((−1)k − 1) = − 2 2 (k = oneven) 2 2 k π k π µ
¶
(11)
Bk T
Bk
+2 2 Z = y(t) sin(kω0 t) dt T T −2
⇓
Z
even × oneven = 0
= 0
(12)
y(t) y(t) =
X A +∞ 4H − cos(2π(2k + 1)f0 t) 2 k=0 (2k + 1)2 π 2
(13)
2.2.2 De zaagtand Gegeven:
Figuur 2: De zaagtand golfvorm en de gekozen periode Gevraagd: A0 ;Ak ;Bk Oplossing: De berekeningen doen we op dezelfde manier als bij de driehoeksgolf. Hier zal Ak nul zijn 6
aangezien we een even functie hebben (als we hem symmetrische t.o.v. de x- en y-as leggen). De bekomen waarden: A0 = H (14) Ak = 0 H kπ
(16)
X H H +∞ − sin(2πkf0 t) 2 k=1 kπ
(17)
Bk = − y(t) =
(15)
2.2.3 De vierkantgolf Gegeven:
Figuur 3: De vierkant golfvorm en de gekozen periode Gevraagd: A0 ;Ak ;Bk Oplossing: We voeren opnieuw de berekeningen uit en bekomen: T 2
(18)
Ak = 0
(19)
A0 =
Bk =
T kπ
(k = oneven)
X T +∞ T y(t) = − sin(2π(2k + 1)f0 t) 4 k=0 (2k + 1)π
(20) (21)
7
2.3
Correlatie
De correlatiecofficint r geeft aan wat de kwaliteit van het lineair verband is. Hij wordt gegeven door: r=
N P
i=1
zX en zY worden gegeven door:
zX zY
N −1
− 1 ≤ r ≤ +1
¯ X −X zX =
σX
Y − Y¯ zY =
(22)
(23)
σY
2.4
Verband FFT-DFS-DFT
1 1 Zoals we weten is een DFT spectrum gespreid over [− f2s , + f2s ] of dus [− 2N , + 2N ]. Wanneer we nu in n periode voor N samples zorgen krijgen we het volgende verband tussen T en N :
T = N · Ts ⇒
N fs 1 = ⇒ f0 = f0 fs N
(24)
Aangezien er in de FFT in het fundamentele interval N samples te vinden zijn zal de stap tussen twee samples dus fNs = f0 zijn. Dit is nog eens duidelijk te zien in onderstaande figuur.
3
Experimenten
3.1
Vraag 1
Deze vraag voerden we reeds uit bij de wiskundige analyse. De conclusies volgen later. 3.2
Vraag 2
Fourier-synthese met matrix-wiskunde Aangezien we in Matlab met matrices werken zullen we voor het samenstellen van de golf
8
uit de verschillende co¨effici¨enten beroep doen op matrix wiskunde. We weten uit (3) dat de reeks een som van sinussen is. Als we dit toepassen in Matlab krijgen we voor y(t):
y(t) =
´ A0 ³ + A1 A2 . . . Ak × cos 2
2πf0 t1 2πf0 t2 2π2f0 t1 2π2f0 t2 .. .. . . 2πkf0 t1 2πkf0 t2
· · · 2πf0 tn · · · 2π2f0 tn .. ... . · · · 2πkf0 tn
+ . . . (25)
,dit kunnen we dan herschrijven tot:
´ A0 ³ y(t) = + A1 A2 . . . Ak · cos 2πf0 2
1 2 .. . k
³ · t1
t2 · · · tn
´ + ...
(26)
Op regel 15 van de code passen we dit toe. Vervolgens bepalen we de correlatie van om te controleren vanaf welke benadering de verschillen tussen benadering en origineel verdwijnen. In het laatste stuk van de code maken we een plot van de eerste vier FS co¨effici¨enten.
5
10
15
20
%e n k e l e g r o o t h e d e n d i e ons s i g n a a l i d e n t i f i c e r e n A=1; f 0 =50; T=1/ f 0 ; f s =8000; %de b e r e k e n d e f o u r i e r c o e f i c i e n t e n A0=A; %h e t a a n t a l f o u r i e r c o e f f i c i e n t e n n=1:21; %de m a t r i x met de c o e f f i c i e n t e n An An=2∗A∗ ( ( − 1 ) . ˆ n − 1 ) . / ( n . ˆ 2 ∗ pi ˆ 2 ) ; %de t i j d waarover we ons s i g n a a a l w i l l e n weergeven %r e k e n i n g houden met h e t a a n t a l s a m p l e s d a t we w i l l e n weergeven t =(0:2∗ f s / f 0 ) / f s ; y=A0/2+An∗ cos (2∗ pi ∗ f 0 ∗n ’ ∗ t ) ; y i d e a a l=p u l s t r a n ( t+T/ 2 , 0 :T: 1 , ’ t r i p u l s ’ ,T ) ; plot ( t , [ y i d e a a l ; y ] ) ; xlabel ( ’ t [ s ] ’ ) ; ylabel ( ’ y ’ ) ; t i t l e ( ’ y ( t ) voor k=21 ’ ) ; %t en y b e p e r k e n t o t h e t s t i j g e n d e d e e l van de p e r i o d e t=t ( 1 :T/4∗ f s +1); y=y ( 1 :T/4∗ f s +1); 9
25
30
35
%b e r e k e n de g e m i d d e l d e s en s t a n d a a r d d e v i a t i e s a v g t=mean( t ) ; s t d t=std ( t ) ; avg y=mean( y ) ; s t d y=std ( y ) ; %b e r e k e n de c o r r e l a t i e c o e f f i c i e n t r=sum( ( ( t−a v g t ) / s t d t ) . ∗ ( ( y−avg y ) / s t d y ) ) / ( length ( t ) −1) %c o r r c o e f ( t , y ) ; figure ; f p =0: f 0 : 3 ∗ f 0 ; stem ( fp , abs ( [ A0 , An ( 1 : 3 ) ] ) ) ; xlabel ( ’ f [ Hz ] ’ ) ; ylabel ( ’ A k ’ ) ; title ( ’A k ’ ); Listing 1: vraag2.m
3.3
Vraag 3
Hier passen we (7) toe om het signaal opnieuw samen te stellen.
5
10
%e n k e l e g r o o t h e d e n d i e ons s i g n a a l i d e n t i f i c e r e n f 0 =50; T=1/ f 0 ; f s =8000; %de t i j d waarover we ons s i g n a a a l w i l l e n weergeven %r e k e n i n g houden met h e t a a n t a l s a m p l e s d a t we w i l l e n weergeven t =0:1/ f s :T; y3=2∗cos (2∗ pi ∗3∗ f 0 ∗ t−pi / 4 ) ; y5 =0.5∗ cos (2∗ pi ∗5∗ f 0 ∗ t −3∗pi / 4 ) ; y=2∗( y3+y5 ) ; plot ( t , [ y ; y3 ; y5 ] ) ; t i t l e ( ’ y ( t )=2\ c d ot ( 2 c o s (2\ p i 3 f 0 t −\p i /4)+0.5 c o s (2\ p i 5 f 0 t −3\ p i / 4 ) ) ’ ) ; xlabel ( ’ f [ Hz ] ’ ) ; ylabel ( ’ y ’ ) ; Listing 2: vraag3.m
3.4
Probleem 1
We passen nogmaals de techniek uit vraag 2 toe. %e n k e l e g r o o t h e d e n d i e ons s i g n a a l i d e n t i f i c e r e n %f 0 b e p a l e n we pas b i j de s y n t h e s e 10
5
10
15
20
25
30
35
H=1; f s =8000; f 1 =50; f 2 =100; %de b e r e k e n d e f o u r i e r c o e f i c i e n t e n A0=H; %h e t a a n t a l f o u r i e r c o e f f i c i e n t e n n=1:12; %de m a t r i x met de c o e f f i c i e n t e n An Bn=−H. / ( n∗ pi ) ; %de t i j d waarover we ons s i g n a a a l w i l l e n weergeven %r e k e n i n g houden met h e t a a n t a l s a m p l e s d a t we w i l l e n weergeven t =0:1/ f s : 1 / f 1 ; %y=s y n t h e s e n=12 en f 0 =50 Hz y=A0/2+Bn∗ sin (2∗ pi ∗ f 1 ∗n ’ ∗ t ) ; subplot ( 2 , 1 , 1 ) ; plot ( t , y ) ; xlabel ( ’ t [ s ] ’ ) ; ylabel ( ’ y ’ ) ; t i t l e ( ’ y ( t )= zaagtand ( 5 0 Hz ) ’ ) ; %de t i j d waarover we ons s i g n a a a l w i l l e n weergeven %r e k e n i n g houden met h e t a a n t a l s a m p l e s d a t we w i l l e n weergeven t =0:1/ f s : 1 / f 2 ; %y=s y n t h e s e n=12 en f 0 =100 Hz y=A0/2+Bn∗ sin (2∗ pi ∗ f 2 ∗n ’ ∗ t ) ; subplot ( 2 , 1 , 2 ) ; plot ( t , y ) ; xlabel ( ’ t [ s ] ’ ) ; ylabel ( ’ y ’ ) ; t i t l e ( ’ y ( t )= zaagtand ( 100 Hz ) ’ ) ; Listing 3: probleem1.m
3.5
Probleem 2
We berekenen de FFT van een sinus- en vierkantsgolf, beide 50 Hz. Voor de sinus zien we makkelijk in dat er slecht 1 co¨effici¨ent Bk zal zijn, namelijk deze bij de grondfrequentie met als waarde de amplitude van de sinus. De vierkantsgolf is reeds in (18) tot (21) berekend. Voor het bepalen van de frequentie-as gebruiken we de techniek uit sectie 2.4.
11
5
10
15
20
25
30
35
40
f 0 =50; f s =8000; T=1/ f 0 ; %a a n t a l s a m p l e s i n 1 p e r i o d e = T/Ts N=f s / f 0 ; %t o v e r 1 p e r i o d e = N s a m p l e s %t=n∗Ts t =(0:N−1)/ f s ; y=sin (2∗ pi ∗ f 0 ∗ t ) ; %FFT o v e r N punten = 1 p e r i o d e Y=f f t ( y ,N ) ; f p=−f s / 2 : f s /N: f s / 2 ; f p=f p ( 1 :N ) ; subplot ( 3 , 1 , 1 ) ; plot ( t ( 1 :N) , y ( 1 :N ) ) ; xlabel ( ’ t [ s ] ’ ) ; ylabel ( ’ y ’ ) ; t i t l e ( ’ y ( t )= s i n (2\ p i 5 0 t ) ’ ) ; subplot ( 3 , 1 , 2 ) ; stem ( fp , f f t s h i f t ( abs (Y) ) /N ) ; xlabel ( ’ f [ Hz ] ’ ) ; ylabel ( ’ | C k | ’ ) ; title ( ’ | C k | ’ ); subplot ( 3 , 1 , 3 ) ; %Bk=−2∗IM( Ck ) stem ( fp , −2∗imag( f f t s h i f t (Y) /N ) ) ; axis ( [ 0 f s /2 0 1 ] ) ; xlabel ( ’ f [ Hz ] ’ ) ; ylabel ( ’ B n ’ ) ; title ( ’B n ’ ); figure ; t =(0:N−1)/ f s ; %s q u a r e neemt een symmetrische v i e r k a n t s g o l f dus we moeten %hem een h a l v e a m p l i t u d e naar boven s c h u i v e n y=T/4+T/4∗ s q u a r e (2∗ pi ∗ f 0 ∗ t ) ; subplot ( 3 , 1 , 1 ) ; plot ( t ( 1 :N) , y ( 1 :N ) ) ; xlabel ( ’ t [ s ] ’ ) ; ylabel ( ’ y ’ ) ; t i t l e ( ’ y ( t )=T/4+T/ 4 . s q u a r e (2\ p i 5 0 t ) ’ ) ; Y=f f t ( y ,N ) ; f p=−f s / 2 : f s /N: f s / 2 ; f p=f p ( 1 :N ) ; 12
45
50
55
subplot ( 3 , 1 , 2 ) ; stem ( fp , f f t s h i f t ( abs (Y) /N ) ) ; xlabel ( ’ f [ Hz ] ’ ) ; ylabel ( ’ | C k | ’ ) ; title ( ’ | C k | ’ ); subplot ( 3 , 1 , 3 ) ; stem ( fp , −2∗imag( f f t s h i f t (Y) /N ) ) ; axis ( [ 0 f s /2 0 T/ pi ] ) ; xlabel ( ’ f [ Hz ] ’ ) ; ylabel ( ’ B n ’ ) ; title ( ’B n ’ ); %druk A0 ; B1 ; B3 a f A 0=Y( 1 ) /N B 1=−2∗imag(Y( 2 ) /N) B 2=−2∗imag(Y( 4 ) /N) Listing 4: probleem2.m
3.6
Probleem 3
Bij vraag 2 pasten we reeds een techniek uit de matrix-wiskunde toe om het signaal opnieuw samen te stellen. Deze gebruiken we nu opnieuw, maar dan in functie van Ck . Voor het verwijderen van de imaginaire componenten (veroorzaakt door de afrondingen), gebruiken we REAL. Bij het tweede signaal passen we dan ook eens IFFT toe.
5
10
15
f =50; T=1/ f ; f s =8000; %a a n t a l s a m p l e s i n 1 p e r i o d e = T/Ts N=f s / f ; %t o v e r 1 p e r i o d e = N s a m p l e s %t=n∗Ts t =(0:N−1)/ f s ; y=sin (2∗ pi ∗ f ∗ t ) ; subplot ( 3 , 1 , 1 ) ; plot ( t , y ) ; xlabel ( ’ t [ s ] ’ ) ; ylabel ( ’ y ’ ) ; t i t l e ( ’ y ( t )= s i n (2\ p i 5 0 t ) ’ ) ; %FFT o v e r N punten = 1 p e r i o d e Y=f f t ( y ,N ) ; %g r a f i e k van de c o e f f i c i e n t e n 13
20
25
30
35
subplot ( 3 , 1 , 2 ) ; f p=−f s / 2 : f s /N: f s / 2 ; f p=f p ( 1 :N ) ; stem ( fp , f f t s h i f t ( abs (Y) ) /N ) ; xlabel ( ’ f [ Hz ] ’ ) ; ylabel ( ’Ck ’ ) ; t i t l e ( ’Ck ’ ) ; %h e t o r i g i n e l e s i g n a a l h e r t e k e n e n subplot ( 3 , 1 , 3 ) ; % N−1 %y=1/N x sum ( Ck . exp ( j 2 p i k n /N) ) % k=0 y=r e a l ( 1 /N∗Y∗exp ( i ∗2∗ pi ∗ ( 0 :N− 1 ) ’ ∗ ( 0 :N−1)/N ) ) ; plot ( t , y ) xlabel ( ’ t [ s ] ’ ) ; ylabel ( ’ y ’ ) ; t i t l e ( ’ y ( t )=1/N\ c d ot \Sigma ˆ{N−1} {k=0} C k\ c d ot e ˆ{2 k\ p i n /N} ’ ) ; figure ;
40
45
50
55
60
y=T/4+T/4∗ s q u a r e (2∗ pi ∗ f ∗ t ) ; subplot ( 3 , 1 , 1 ) ; plot ( t , y ) ; xlabel ( ’ t [ s ] ’ ) ; ylabel ( ’ y ’ ) ; t i t l e ( ’ y ( t )=T/4+T/ 4 . s q u a r e (2\ p i 5 0 t ) ’ ) ; %FFT o v e r N punten = 1 p e r i o d e Y=f f t ( y ,N ) ; %g r a f i e k van de c o e f f i c i e n t e n subplot ( 3 , 1 , 2 ) ; f p=−f s / 2 : f s /N: f s / 2 ; f p=f p ( 1 :N ) ; stem ( fp , f f t s h i f t ( abs (Y) ) /N ) ; xlabel ( ’ f [ Hz ] ’ ) ; ylabel ( ’ C k ’ ) ; title ( ’C k ’ ); %h e t o r i g i n e l e s i g n a a l h e r t e k e n e n subplot ( 3 , 1 , 3 ) ; y=i f f t (Y,N ) ; plot ( t , y ) ; xlabel ( ’ t [ s ] ’ ) ; ylabel ( ’ y ’ ) ; t i t l e ( ’ y ( t )=1/N\ c d ot \Sigma ˆ{N−1} {k=0} C k\ c d ot e ˆ{2 k\ p i n /N} ’ ) ; 14
Listing 5: probleem3.m
3.7
Probleem 4
Met het commando load(’Oef2Prob4.mat’) laden we de gegevens in van x en y. Daarvan bepalen we de FFT en halen er de fc en A uit. Dit doen we met behulp van de geplotte grafieken.
5
10
15
load ( ’ Oef2Prob4 . mat ’ ) ; f s =44100; N=length ( x ) ; X=f f t ( x ,N ) ; Y=f f t ( y ,N ) ; f p =(−(1−1/N) / 2 : 1 /N: 1 / 2 ) ∗ f s ; f p=f p ( 1 :N ) ; subplot ( 2 , 1 , 1 ) ; subplot ( 2 , 1 , 1 ) ; stem ( fp , abs ( f f t s h i f t (X) ) /N ) ; xlabel ( ’ f [ Hz ] ’ ) ; ylabel ( ’X ’ ) ; t i t l e ( ’X( f ) ’ ) ; subplot ( 2 , 1 , 2 ) ; stem ( fp , abs ( f f t s h i f t (Y) ) /N ) ; xlabel ( ’ f [ Hz ] ’ ) ; ylabel ( ’Y ’ ) ; t i t l e ( ’Y( f ) ’ ) ; A=abs (Y( 2 ) ) . / abs (X( 2 ) ) Listing 6: probleem4.m
4
Resultaten
4.1
Vraag 1
De co¨effici¨enten berekenden we in de wiskundige analyse. Wanneer we deze golf filteren op 3, 5ω0 zullen dus alle harmonischen boven deze waarde weggefilterd worden. Wat dus wil zeggen dat we enkel de harmonischen tot 3ω0 kunnen weergeven, want k is een geheel getal. De co¨effici¨enten zijn dan A0 tot A3 met dezelfde waarden als de ongefilterde versie.
15
|A | k
1
0.9
0.8
0.7
k
|A |
0.6
0.5
0.4
0.3
0.2
0.1
0
0
50
100
150
f[Hz]
Figuur 4: De eerst vier FS co¨effici¨enten van vraag 2 4.2
Vraag 2
In figuur 4 zijn de eerste vier FS co¨effici¨enten en A0 weergegeven. In figuren 5 tot en met 8 staan de grafieken van enkele benaderingen.
16
y(t) voor k=3 1
0.9
0.8
0.7
y
0.6
0.5
0.4
0.3
0.2
0.1
0
0
0.002
0.004
0.006
0.008
0.01
0.012
0.014
0.016
0.018
0.02 t[s]
Figuur 5: Derde term fourier benadering van 50 Hz driehoeksgolf
y(t) voor k=7 1
0.9
0.8
0.7
y
0.6
0.5
0.4
0.3
0.2
0.1
0
0
0.002
0.004
0.006
0.008
0.01
0.012
0.014
0.016
0.018
0.02 t[s]
Figuur 6: Zevende term fourier benadering van 50 Hz driehoeksgolf
17
y(t) voor k=12 1
0.9
0.8
0.7
y
0.6
0.5
0.4
0.3
0.2
0.1
0
0
0.002
0.004
0.006
0.008
0.01
0.012
0.014
0.016
0.018
0.02 t[s]
Figuur 7: Twaalfde term fourier benadering van 50 Hz driehoeksgolf y(t) voor k=21 1
0.9
0.8
0.7
y
0.6
0.5
0.4
0.3
0.2
0.1
0
0
0.002
0.004
0.006
0.008
0.01 t[s]
0.012
0.014
0.016
0.018
0.02
Figuur 8: Eenentwintigste term fourier benadering van 50 Hz driehoeksgolf Voor deze benaderingen berekenden we volgende waarden voor de correlatieco¨effici¨enten: k r 3 0.9951 7 0.9992 12 0.9997 21 0.9999 We zien dat we dus al bij k = 3 een goede benadering hebben. Als we dit bekijken in de grafiek zien we dat we inderdaad een vrij goede benadering hebben. Het probleem zit 18
hem echter niet in het lineaire deel, dat we al bij lage waarden voor k goed benaderen, maar in de toppen. In de toppen hebben we met het zogenaamde Gibbs-fenomeen te maken. Dit houdt in dat in discontinu¨ıteiten we met overshoot te maken hebben. Bij deze driehoeksgolf levert dit een afgeronde top op. Bij een blokgolf zou dit pieken opleveren in de hoekpunten. Doordat deze top een geringer aantal punten vertegenwoordigt dan het lineaire deel, hebben we dus al vlug een zeer hoge correlatieco¨effici¨ent. Verder zal deze reeks vlugger convergeren (betere benaderingen) dan bijvoorbeeld een zaagtand of blokgolf door het kwadraat in de noemer. Voor k = 21 is de afronding van de toppen verwaarloosbaar geworden en kunnen we dit dus als een zeer goede benadering zien, wat tevens uitgedrukt is in de correlatieco¨effici¨ent. 4.3
Vraag 3
In figuur 9 is samengestelde golfvorm te zien. Hierin zijn ook de twee harmonischen weergegeven. y(t)=2⋅(2cos(2π3f t−π/4)+0.5cos(2π5f t−3π/4)) 0
0
5 y(t) n=3 n=5
4
3
2
y
1
0
−1
−2
−3
−4
−5
0
0.002
0.004
0.006
0.008
0.01 f[Hz]
0.012
0.014
0.016
0.018
0.02
Figuur 9: De golfvorm van vraag 3 met zijn twee harmonischen Indien we nu de vijfde harmonische willen verwijderen met een ideaal laagdoorlaat filter moeten we een afsnijfrequentie nemen die kleiner is dan 5f0 want dan wordt de vijfde 19
harmonische weggefilterd. We moeten ze echter groter nemen dan 3f0 , want anders zou ook de derde harmonische worden weggefilterd. 4.4
Voorbeelden
Als we N laten stijgen neemt de resolutie van onze DFT toe. Dit heeft wel enkele gevolgen. Wanneer we voor N een geheel aantal keer de periode nemen van ons periodiek uitgebreid signaal (zie eerste van voorbeeld 2), krijgen we de juiste FS co¨effici¨enten. We krijgen in ons spectrum dan wel extra nullen tussen de co¨effici¨enten door te toegenomen resolutie. Dit kan dus verwarrend werken. Wanneer we ons signaal aanvullen met nullen, ”zero padding”, komt dit neer op een venster op de functie plaatsen. Daardoor krijgen we wel een hogere resolutie maar merken we sinc-functies rond de twee impulsen. Dit stemt overeen met de in de theorie berekende frequentieresponsie van een venster. Als we de FFT bepalen over een deel van de periode komt dit ook neer op het plaatsen van een venster. De resolutie zal nu zelf dalen. We kunnen zelf de juiste hoogte van de twee impulsen niet meer onderscheiden doordat deze tussen twee samples ligt en dus niet wordt weergegeven. 4.5
Probleem 1
De in de wiskundige analyse berekende co¨effici¨enten zijn onafhankelijk van de frequentie. Ze zijn enkel afhankelijk van de natuurlijke index k. Dit wil dus zeggen dat deze co¨effici¨enten voor elke zaagtand van dezelfde hoogte dezelfde zijn. Waar liggen dan de verschillen? Zoals te zien is in figuur 10 is de periode van de tweede zaagtand logischer wijs de helft. We zien dat we bij de tweede zaagtand echter minder details hebben (minder samples per periode). Dit komt doordat we voor beide signalen dezelfde samplingrate genomen hebben. 4.6
Probleem 2
De gevonden spectra zijn weergegeven in figuren 11 en 12. We gaan ze nu vergelijken met onze berekende waarden.
20
y(t)=zaagtand(50 Hz) 1.2 1 0.8
y
0.6 0.4 0.2 0 −0.2
0
0.002
0.004
0.006
0.008
0.01 t[s]
0.012
0.014
0.016
0.018
0.02
0.007
0.008
0.009
0.01
y(t)=zaagtand(100 Hz) 1.2 1 0.8
y
0.6 0.4 0.2 0 −0.2
0
0.001
0.002
0.003
0.004
0.005 t[s]
0.006
Figuur 10: De synthese van de zaagtand van probleem 1 y(t)=sin(2π50t)
y
1
0
−1
0
0.002
0.004
0.006
0.008
0.01 t[s]
0.012
0.014
0.016
0.018
0.02
|Ck|
|Ck|
1
0.5
0 −4000
−3000
−2000
−1000
0 f[Hz]
1000
2000
3000
4000
2500
3000
3500
4000
Bk
Bk
1
0.5
0
0
500
1000
1500
2000 f[Hz]
Figuur 11: Analyse van de sinusgolf uit probleem 2
21
Voor de sinusgolf is het logisch dat we enkel een B1 zullen hebben met de waarde van de amplitude (de som van sinussen is immers de sinus zelf). Dit zien we dan ook duidelijk in figuur 11. Aangezien Bk = −2 · ℑ{Ck }, zal |Ck | dus de helft zijn van Bk . y(t)=T/4+T/4.square(2π50t)
y
0.01
0.005
0
0
0.002
0.004
0.006
0.008
0.014
0.016
0.018
0.02
x 10
4
k
|C |
0.012
|Ck|
−3
6
0.01 t[s]
2 0 −4000
−3000
−2000
−1000
0 f[Hz]
1000
2000
3000
4000
2500
3000
3500
4000
Bk
−3
x 10 6 Bk
4 2 0
0
500
1000
1500
2000 f[Hz]
Figuur 12: Analyse van de vierkantgolf uit probleem 2 Voor de vierkantgolf berekenden we in (18) waarden: A0 = C0 = 2 T B1 = = kπ T = B3 = kπ
A0 , (20) Bk . We bepalen hiervoor de volgende 0.02 = 0, 0005 4 0.02 = 0, 0064 π 0.02 = 0, 0021 3π
We zien in figuur 12 dat deze waarden overeenstemmen. Als controle lieten we Matlab nog eens deze waarden teruggeven, we kregen: >> A_0 = 0.0050 22
B_1 = 0.0064
B_2 = 0.0021 4.7
Probleem 3
De resultaten zijn te vinden in figuur 13 en 14. y(t)=sin(2π50t)
y
1
0
−1
0
0.002
0.004
0.006
0.008
0.01 t[s]
0.012
0.014
0.016
0.018
0.02
Ck
Ck
1
0.5
0 −4000
−3000
−2000
−1000
0 f[Hz]
1000
2000
3000
4000
y(t)=1/N⋅ΣN−1 C ⋅ e2kπn/N k=0 k
y
1
0
−1
0
0.002
0.004
0.006
0.008
0.01 t[s]
0.012
0.014
0.016
0.018
0.02
Figuur 13: Synthese van de sinusgolf uit probleem 3
23
y(t)=T/4+T/4.square(2π50t)
y
0.01
0.005
0
0
0.002
0.004
0.006
0.008
0.012
0.014
0.016
0.018
0.02
Ck
−3
6
0.01 t[s]
x 10
Ck
4 2 0 −4000
−3000
−2000
−1000
0 f[Hz] N−1
1000
2000
3000
4000
2kπn/N
y(t)=1/N⋅Σk=0 Ck⋅ e 0.02
y
0.01 0 −0.01
0
0.002
0.004
0.006
0.008
0.01 t[s]
0.012
0.014
0.016
0.018
0.02
Figuur 14: Synthese van de vierkantgolf uit probleem 3 Wanneer we ons aantal DFS co¨effici¨enten laten stijgen, laten we dus eigenlijk N stijgen. Onze synthese zal dus meer samples bevatten. Het resulterende signaal zal dus een hogere samplingrate hebben. Bij probleem 1 zal dit een ander effect geven, aangezien we daar de co¨effici¨enten interpreteren als amplitudes en fasen van sinussen. Hoe meer co¨effici¨enten we dan gebruiken hoe lager de rimpel zal zijn, dus hoe beter we het originele signaal benaderen. Dit konden we ook zien bij vraag 2 waar we bepaalden hoever we moesten gaan om een goede benadering te hebben. 4.8
Probleem 4
In figuur 15 zijn ons signaal x en gefilterd signaal y weergegeven.
24
x(t) 1.5 1
x
0.5 0 −0.5 −1 −1.5
0
1
2
3 t[s]
4
5
6 −3
x 10
y(t) 1.5 1
y
0.5 0 −0.5 −1 −1.5
0
1
2
3 t[s]
4
5
6 −3
x 10
Figuur 15: Golfvormen van probleem 4
25
X(f) 0.4
X
0.3 0.2 0.1
−2000
−1000
0 f[Hz] Y(f)
1000
2000
3000
4000
2000
3000
4000
784 Hz 980 Hz
−3000
−980 Hz −784 Hz
0 −4000
0.4
Y
0.3 0.2 0.1 0 −4000
−3000
−2000
−1000
0 f[Hz]
1000
Figuur 16: Analyse van probleem 4 In figuur 16 zijn de fourier-spectra getekend van onze beide signalen. We kunnen daaruit afleiden dat vanaf 980 Hz gefilterd is. Onze afsnijfrequentie moet dus zeker onder deze waarde liggen. Aangezien de vorige harmonische wel nog doorgelaten is, 784 Hz, moet de afsnijfrequentie dus boven deze waarde liggen. Dus: 784 Hz < fc < 980 Hz ,of: 4926, 1 Hz < ωc < 6157, 5 Hz Als versterkingsfactor krijgen we: Als >> A = 1.5000
26
5
Besluit
Dankzij het berekenen van de Fourierreeks van enkele golfvormen met de hand, hebben we ons inzicht in de samenstelling van de Fourierreeks kunnen verbreden. Zo kan je aan de aard van de periodieke golf al uitspraken doen over de te bekomen co¨effici¨enten: • • •
even functie f (t) = f (−t): enkel cos-termen, dus Bk = 0 oneven functie f (t) = −f (−t): enkel sin-termen, dus Ak = 0 gemiddelde is nul: A0 = 0
Hoe meer FS co¨effici¨enten we in onze reeks opnemen, hoe nauwkeuriger onze benadering. Vanzelfsprekend bereiken we in de ”theoretische”oneindigheid de originele functie. In ons voorbeeld kwamen we al uit op een goede benadering bij een zevental co¨effici¨enten. We controleerden dit met de correlatieco¨effici¨ent, alhoewel men zich hier de vraag kan stellen of dit wel de aangewezen manier is. Het probleem zit hem immers in de afgeronde toppen van de driehoeksgolfvorm (Gibb’s fenomeen) en dit wordt niet voldoende doorgedrukt in de correlatieco¨effici¨ent. Een betere manier was misschien geweest om de kwadratische fout te berekenen, wat eenvoudig kan via de identiteit van Parseval. Bij het werken met de FFT dient voldoende aandacht besteed te worden aan het gebruik van de lengte N . Het is zo dat wanneer we deze groter nemen dan de periode of kleiner dan een periode, we eigenlijk een venster plaatsen op onze gesamplede functie (bij grotere N kan dit ook zero padding zijn). Dit geeft tot gevolg dat onze pulsen, Dirac, in ons spectrum omgezet worden in sinc-functies. Bij kleine afwijkingen op de lengte kunnen we nog duidelijk de pulsen onderscheiden aan de toppen van de sinc, maar indien dit teveel afwijkt krijgen we met overlap tussen sinc’s te maken en kunnen twee pulsen bijvoorbeeld als ´e´en sinc zichtbaar worden. Verder dient opgemerkt te worden dat bij de DFT de breedte van het spectrum altijd de samplingfrequentie is. Dit heeft te maken met het Nyquisttheorema. In de praktijk heeft de FFT zijn nut bewezen bij het analyseren van filters. Dankzij deze FFT kunnen we uit het spectrum de afsnijfrequentie en versterkingsfactoren van ideale filters berekenen. Bij het bepalen van de afsnijfrequentie moet er wel opgelet worden. De afsnijfrequentie moet natuurlijk lager zijn dan de laagste harmonische die niet doorgelaten wordt, maar moet ook groter zijn dan de hoogste harmonische die wel doorgelaten wordt. Er is dus een gebied waarbinnen de afsnijfrequentie kan liggen. Tot slot nog iets over probleem 5. We zijn in de labosessies begonnen aan deze opgave, maar doordat we tijd te kort hadden om ze nog af te werken en we geen samplingfrequentie hadden, hebben we dit probleem niet meer ten gronde kunnen oplossen. Een mogelijke oplossing is zoals in de vorige problemen de FFT berekenen. Dit levert ons van het origineel en het gesamplede signaal de reeks co¨effici¨enten. Daarvan kunnen we met abs en angle de verandering in magnitude en fase bepalen. Het samenstellen van het signaal kunnen we dan weer doen zoals in de voorgaande problemen.
27
