Opdracht 1 – Systeemresponsies en signaaltransforma- ties
Deze MATLAB opdracht is onderdeel van de cursus Digitale Systemen (ELEDIS10). Zoals staat beschreven in de cursushandleiding moeten de resultaten van deze opdracht worden opgenomen in een verslag1.
De opdracht moet individueel2 worden uitgevoerd. Hiertoe worden voor elke student andere input geproduceerd (zie hoofdstuk 1). Neem in het verslag (in PDF-formaat) de antwoorden op de vragen voor deze opdracht op in een apart hoofdstuk, eventueel toegelicht met figuren.
Geef duidelijk d.m.v. nummering in het verslag aan op welke vraag je antwoord geeft. Lever het verslag samen met de MATLAB.mat-files in een ZIP-bestand in, waarin variabelen met filtercoëfficiënten en de gefilterde output zijn opgeslagen. Een meer specifieke omschrijving staat in de opdrachten zelf.
Verwerk naam en studentnummer in de naam van de ZIP-file (0999999_Achternaam.zip) en lever deze in via de inlevermap (.../2019-2020/ELEDIS10/) vóór de deadline, welke is vastgesteld op:
zondag 5 april 2020 om 23.59 uur.
Werk dat ná dit moment ingeleverd wordt, zal niet worden beoordeeld en de student in kwestie scoort dan het cijfer 0 voor het verslag van ELEDIS10 en dit dient herkanst te worden conform de beschrijving in de cursushandleiding.
Theoretische vragen over de opdracht dienen te worden gericht aan de eigen docent. Maak een issue3 aan indien er bugs of andere zaken m.b.t. de opdracht of software te rapporteren zijn.
1 Voorbereiding
Voor de vragen worden de getallen van je studentnummer als volgt gehanteerd:
0c1c2c3c4c5c6
Bijvoorbeeld 0987654 geeft:
c1= 9 c2= 8 c3= 7 c4 = 6 c5 = 5 c6= 4
Z.O.Z. VOOR DE OPDRACHTEN!
1 De eisen die gesteld worden aan dit verslag kun je vinden op:https://bitbucket.org/HR_ELEKTRO/dis10/
wiki/Verslag/DIS10_Verslag_1_Eisen.pdf.
2 Vermoedens van plagiaat of fraude worden gemeld bij de Examencommissie volgens het “Reglement tentaminering en examinering” (Artikel 11 van de Hogeschoolgids EAS).
3 Je kunt een issue aanmaken op:https://bitbucket.org/HR_ELEKTRO/dis10/issues.
2 Impuls- en stapresponsie
1.1
Gegeven de volgende recurrente betrekking:y[n] = (0,c1)y[n − 1] − (0,c2)y[n − 2] + (0,c3)y[n − 3] + (0,c4)x[n]
− (0,c5)x[n − 1] − (0,c6)x[n − 2]
A
Bepaal met handberekeningen de eerste 4 waarden van de impulsresponsie h[n]van deze recurrente betrekking waarbij je de getallen van jouw studentnummer hebt ingevuld.
Bijvoorbeeld, als je studentnummer 0987654 is:
y[n] = 0,9y[n−1]−0,8y[n−2]+0,7y[n−3]+0,6x[n]−0,5x[n−1]−0,4x[n−2]
B
Gebruik MATLABom de impulsresponsie te plotten. Allereerst moeten we de coëf- ficiënten bk van x[n − k] en de coëfficiënten ak van y[n − k] worden ingevoerd.Bijvoorbeeld, als je studentnummer 0987654 is:
>> b=[0.6 , -0.5, -0.4];
>> a=[1 , -0.9, 0.8 , -0.7];
Merk op dat het teken van de coëfficiënten van y[n − k] is omgedraaid. Dit is omdat MATLABmet de overdrachtsfunctie werkt:
H(Ω) =Y(Ω) X(Ω) =
PM
m=0bme− jΩm 1−PK
m=1ame− jΩm
Plot met behulp van de onderstaande functie de impulsresponsie:
>> impz(b,a, aantal_gewenste_samples );
Controleer of de eerste waarden van je handberekeningen overeenkomen met de eerste waarden in de plot. Is het systeem stabiel?
C
Plot met behulp van onderstaande functie de stapresponsie van het systeem:>> stepz(b,a, aantal_gewenste_samples );
1.2
Gebruik de onderstaande recurrente betrekking voor de volgende vragen:y[n] =¨(0,5 + 0,c5)y[n − 1] + x[n] als c5 < 5 (−0,4 + 0,c5)y[n − 1] + x[n] als c5≥ 5
A
Bepaal met behulp van een handberekening de eindwaarde van de stapresponsie van deze recurrente betrekking.B
Gebruik MATLABom te verifiëren of je handberekening correct is.3 Voortschrijdend gemiddelde filter
Het blokdiagram van een systeem is te zien in figuur 1.
h[n]
x[n] y[n]
Figuur 1: Het systeem
Gebruik bij de onderstaande vragen het volgende ingangssignaal voor het systeem:
x[n] = sin 2πn
c1c2c3c4c5c6 mod 1000
+ sin
2πn
bc1c2c1003c4c5c6c
1.3
Wat doen demodenflooroperaties precies?1.4
Wat is de periode van de laagste frequentie in dit signaal? Plot precies 2 van deze periodes van x[n]. Stel dat je studentnummer 0987654 is en de periode van de laagste frequentie 1000 samples is (wat in het bovenstaande voorbeeld niet correct is) dan kan de volgende code worden gebruikt:>> n =1:2*1000;
>> x=sin(2*pi*n/(mod(987654 ,1000) )) + ←- ,→ sin(2*pi*n/(floor(987654/100) ));
>> scatter(n,x);
1.5
We willen dat de hoogste frequentie van dit signaal wordt weggefilterd. Dit kan met be- hulp van een voortschrijdend-gemiddelde (moving-average) filter. De impulsresponsie van zo’n filter is als volgt:h[n] =
1
d 0≤ n < d 0 elders
Waarbij d het gewenste aantal samples is waarover het gemiddelde moet worden genomen. Dit kan als volgt in MATLAB worden gezet.
>> b = zeros(1,d)
>> for k = 1:length(b)
>> b(k) = 1/d;
>> end
Bepaal wat de hoogste frequentie van het ingangssignaal is en filter deze weg met behulp van onderstaande code:
>> y=filter(b,1,x);
Hoe komt het dat deze frequentie wordt weggefilterd door dit filter?
1.6
Verifieer visueel of de frequentie is weggefilterd met behulp van het volgende stukje code:>> scatter(n,y)
Geef een verklaring voor het resultaat.
4 Convolutie
Voeg de file wannabe_in_la.wav (te vinden op de wiki) toe aan je working directory van MATLAB. Lees dit bestand in met behulp van het volgende commando:
>> [x,Fs ]=audioread('wannabe_in_la .wav ');
In de variabelex staan nu de samples van dit audiobestand en in Fs de samplefrequentie van dit audiobestand. We gaan dit audiobestand in 2 domeinen filteren. Allereerst in het tijddomein en vervolgens in het frequentiedomein. Na afloop van deze bewerkingen bepalen we de tijd die de CPU nodig heeft gehad om deze bewerking uit te voeren. We bekijken ook de frequentieresponsie van de twee filters en het frequentiespectrum van het gefilterde en ongefilterde audiobestand. In MATLAB is het fragment nu af te spelen doorsoundscte gebruiken:
>> soundsc(x,Fs);
1.7
We beginnen met het audiobestand te filteren in het tijddomein:A
Bepaal de lengte (in seconden) van dit audiobestand op basis van het aantal samples en de samplefrequentie.B
Een eenvoudig laagdoorlaatfilter wordt in het tijddomein gespecificeerd door de volgende impulsresponsie:h[n] =sinπn
p
πn
Hier is p de variabele waarmee de (genormaliseerde) kantelfrequentie van dit filter kan worden ingesteld. Bepaal de waarde van p om een kantelfrequentie
fc= 100 · (c4+ 1) + 200 Hz te krijgen.
p= fs 2 fc
C
Onderstaande MATLABcode genereert 2k+ 1 coëfficiënten van h[n] waarbij in dit voorbeeld p= 20 is gebruikt.>> k = 1e5;
>> n = -k:k;
>> b = sin(pi*n/20) ./(pi*n);
>> b(n==0) = 1/20;
De laatste regel is noodzakelijk omdat h[n] een deling door nul oplevert voor n = 0.
Er zal dus zelf berekend moeten worden wat lim
n→0h[n] is! Geef deze berekening en genereer de coëfficiënten voor de h[n] van jouw filter.
D
Controleer het frequentiespectrum van je filter met behulp vanfreqz. Omdat er geen a coëfficiënten zijn wordt hiervoor simpelweg een 1 ingevuld, waarom wordt hier geen 0 ingevuld?>> freqz(b ,1 ,0: Fs/2,Fs);
E
Varieer met het aantal coëfficiënten van de impulsresponsie in ordes van 10, wat is het resultaat van deze variatie?F
Met behulp van het commandoticwordt een timer gestart en met het commando tocwordt deze timer gestopt en wordt de verstreken tijd in het command window van MATLAB afgedrukt. Filter het audiobestand door de impulsresponsie h[n] te convolueren met de samples van het audiobestand met behulp van onderstaande code waarbijk=1e5 en noteer de tijd dat de CPU hiervoor nodig heeft.>> y=conv(b,x);
G
Luister nogmaals naar het signaal en stel vast of het filter zijn werk goed heeft gedaan met behulp van onderstaande code:>> fft_y =fft(y);
>> % Calculate the frequencies in the array to be plotted as ←- ,→ x-axis
>> df = Fs/length( fft_y );
>> f = (0:length( fft_y ) -1)*df;
>> plot(f, 20*log10(abs( fft_y )));
>> grid;
Deze code transformeert het gefilterde audiosignaal naar het frequentiedomein en plot dit spectrum om te kunnen verifiëren of het filter correct is toegepast. Als het resultaat naar wens is kunnen de coëfficiëntenb en het gefilterde signaal y worden opgeslagen door ze te exporteren naar een.mat-file (pas het studentnummer aan):
>> save('M1deel1_0987654 .mat ','b','y');
1.8
Nu gaan we hetzelfde filter toepassen op het audiobestand maar dan in het frequentie- domein:A
De overdrachtsfunctie van laagdoorlaatfilter met een rechthoekige (brickwall) fre- quentieresponsie in het frequentiedomein is als volgt gespecificeerd:H(Ω) =
¨1 als 0≤ Ω ≤πp 0 als πp ≤ Ω ≤ π
Waarbij p weer de (genormaliseerde) kantelfrequentie bepaalt.
B
Maak een laagdoorlaat filter in het frequentiedomein met de fc en p zoals gevon- den bij opdracht 1.7 deelopdracht B. Controleer met behulp van een plot of de frequentieresponsie van het filter aan de eisen voldoet.C
Filter het audiobestand door de frequentieresponsie H(Ω) te vermenigvuldigen met het frequentiespectrum van het audiobestand en noteer hoe lang de CPU hierover doet. Vergelijk deze tijd dit met het in opdracht 1.7 deelopdracht F gevonden antwoord. Is er een verschil en zo ja wat is de belangrijkste veroorzaker van dit verschil? Het frequentiespectrum van het audiobestand kan je berekenen met de volgende MATLAB code:>> freq_x = fft(x);
D
Stel vast of het filter zijn werk goed heeft gedaan met behulp van onderstaande code:>> %1e plot is het gefilterde frequentiespectrum in dB
>> subplot(2 ,1 ,1)
>> plot(f, 20*log10(abs( freq_y )));
>> %2e plot is het originele frequentiespectrum in dB
>> subplot(2 ,1 ,2)
>> plot(f, 20*log10(abs( freq_x )));
Vergelijk het gefilterde frequentiespectrum met het frequentiespectrum van op- dracht 1.7 deelopdracht G. Wat valt je op?
E
Luister nogmaals naar het signaal door het frequentiespectrum terug te transformeren naar het tijddomein met behulp van de inverse fourier transformatie functieifft().Als het resultaat naar wens is kan de overdrachtsfunctie en het gefilterde signaal in het tijddomeintime_y worden opgeslagen door ze te exporteren naar een .mat-file (pas het studentnummer aan):
>> save('M1deel2_0987654 .mat ','H','time_y ');
