informatie   rooster   hoorcollege   practicum   links

Fouriertransformatie

In dit practicum komt een groot aantal aspecten aan bod van het werken met de discrete Fouriertransformatie (DFT), en de efficiënte implementatie daarvan in de vorm van de FFT.

Er wordt vanuit gegaan dat je je degelijk hebt voorbereid door de teksten over de Fouriertransformatie en het tijd-frequentievlak grondig door te lezen. Het practicum omvat nagenoeg alle onderwerpen die hierin aan bod gekomen zijn en borduurt er hier en daar op voort (in het geval van windows bijvoorbeeld).

De uitvoering van het practicum is eenvoudig. Selecteer alle Matlab-commando's die in fixed-font zijn weergegeven in de volgorde zoals aangegeven. Doorgaans vormen ze aparte paragrafen. Wanneer je een commando vergeet zal Matlab meestal even later klagen over niet-gedefinieerde waarden. Herhaal dan de voorgaande commando's. Zorg ervoor dat je de commando's niet blind uitvoert, maar probeer te begrijpen wat je doet. In veel gevallen geven de Matlab-commando's je al een deel van het antwoord.

Matlab is overigens een computertaal die in de wetenschap en de technologie zeer veel gebruikt wordt om prototypen te bouwen, data te analyseren en systemen te simuleren.

Voor het practicum is een drietal files nodig die in je Matlab "path" moeten komen te staan. Maak bijvoorbeeld ~/matlab/spraak en zet daarin de files uit practFFT.zip:

    fftGram.m
    logsweep.m
    mfcc.m
    winHamming.m
    winHanning.m
    nul123.wav
    

Start Matlab en ga naar je workdirectory (werkbalkje bovenin). Voer het onderstaande commando uit door het via "control-C" en "control-V" naar het Matlab commandovenster te kopiëren.

    format compact                  % Haalt nodeloze lege regels weg
    

onderwerpen

• basisvectoren van de DFT
• DFT
• complexe DFT - de FFT
• impuls-respons van de (I)DFT
• windows (vensters)
• logsweep
• spectrogram
• MFCC (Mell-scaled Cepstral Coefficients)

Basisvectoren van de DFT

De reële discrete Fouriertranformatie (DFT) is besproken in de tekst over de Fouriertransformatie. We gaan eerst een matix maken waarin alle basisvectoren staan. Laten we uitgaan van een N = nfft = 32 dimensionale inputvector. Als basisvectoren hebben we in het tijddomein cos(ω0), cos(ωk) en sin(ωk) voor k = 1:N/2-1, en cos(ωN/2). We kunnen deze 32 basisvectoren in een NxN-matrix stoppen.

Voer onderstaande regels uit in het Matlab commandovenster.

nfft = 32;
B = zeros(nfft,nfft);
k = 0;
B(:,1) = cos(2*pi*k*(0:nfft-1)/nfft)' / sqrt(nfft);
for k = 1:nfft/2-1
  B(:,k*2) = sin(2*pi*k*(0:nfft-1)/nfft)' / sqrt(nfft/2);
  B(:,k*2+1) = cos(2*pi*k*(0:nfft-1)/nfft)' / sqrt(nfft/2);
end
B(:,nfft) = cos(2*pi*nfft/2*(0:nfft-1)/nfft)' / sqrt(nfft);

n = [1 2 3 4 5 nfft]; % basisvector nummer
plot(B(:,n), '-*'); axis([0 nfft+1 [-1.1 1.1]*max(B(:))])
title('Basisvectoren ')

opdracht A

1. Bepaal op basis van het bovenstaande algoritme de default kleur volgorde in Matlab is, dus: {blauw, ... , ... }
2. Waarom staan in bovenstaande algoritme de eerste en de nfft-de basisvector buiten de loop?
3. Geef voor zowel de cosinus- als de sinusbijdragen een formule, waarin het golfgetal k gekoppeld wordt aan kolomnummer (d.w.z. basisvectornummer) n van B.
4. Gegeven een samplefrequentie fs, geef formules voor de frequentie van de basisvectoren, en voor het tijdstip van elk van de samplewaarden in de basisvectoren (start op t = 0 ms). (Zie hier.)
5. Bij een matrixvermenigvuldiging A1*A2 = C geldt dat matrixelementen C(i,j) bepaald worden door de som van het elementsgewijze product van rij i van A1 en kolom j van A2. Geef aan waarom het matlab commando

   imagesc(B*B'); colorbar

   aantoont dat voor B voldaan wordt aan de eisen voor basisvectoren.

DFT

Nu we de basisvectoren hebben kunnen we de FT uitvoeren op nfft-dimensionale inputvectoren. Laten we eens twee pulsen nemen: p1 met de puls op positie 1 en p5 met de puls op positie 5.

nfft = 32;
p1 = zeros(1,nfft); p1(1) = 1; % puls op vectorpositie 1
p5 = zeros(1,nfft); p5(5) = 1; % puls op vectorpositie 5

Zoals uitgelegd zijn de bijdragen van de Fourier-basisvectoren te berekenen als inproduct tussen de input en de basisvectoren, ofwel:

P1_direct = zeros(1,nfft); P2_direct = zeros(1,nfft); % initialisaties om f-domein repr. te alloceren
for ii = 1:nfft
  P1_direct(ii) = sum(p1.*B(:,ii)'); % inproduct op basisvector
  P5_direct(ii) = sum(p5.*B(:,ii)');
end

In Matlab is er een compactere en intuïtieve schrijfwijze die overeen komt met de definitie:

P1 = p1*B; % DFT van p1
P5 = p5*B; % DFT van p5

Check dat beide berekeningswijzen tot identieke resultaten leiden (P1==P1_direct).

opdracht B

1. Voor de terugtransformatie, doorgaans de inverse Fouriertransformatie genoemd, moet gelden dat de inverse transformatie van P1 en P5 identiek moeten zijn aan p1 en p5. Leg uit waarom dit blijkt uit:

sum(abs(P1*B' - p1))
sum(abs(P5*B' - p5))

Vergelijk het resultaat met eps*nfft (eps is het kleinste positieve getal ongelijk 0 dat Matlab kan representeren).
N.B. Dit soort kleine getallen komen regelmatig voor bij het gebruik van de FT en kunnen er bij de complexe FT voor zorgen dat een resultaat niet strict reeel is, terwijl dat wiskundig gezien wel het geval moet zijn.

Het opbouwen van p1 uit cosinus en sinus basisvectoren:

for n = 1:nfft
  stem(P1(1:n)*B(:,1:n)');
  axis([0 nfft+1 -1 1])
  title(['Construction of pulse P1 using n=1:' int2str(n)])
  pause(10./nfft)
end

Het opbouwen van p5 uit cosinus- en sinusbasisvectoren:

for n = 1:nfft
  clf
  stem(P5(1:n)*B(:,1:n)'); hold on
  plot(((P5(1:n)'*(ones(1,nfft))).*B(:,1:n)')'); hold off
  axis([0 nfft+1 -1 1])
  title(['Construction of pulse P5 using n=1:' int2str(n)])
  pause(10./nfft) % pause zonder tijdindicatie in sec wacht op toets
end

opdracht C

1. Voer bovenstaande FOR-lussen nog eens uit, maar selecteer alleen het deel tot het woord pause, d.w.z. laat het deel (10/nfft) weg. Type daarna ; end om de FOR-lus af te maken en de berekening te starten. Doordat de het pause-commando nu geen tijd-indicatie heeft moet je na elke basisvectorbijdrage op een toets drukken. Kijk goed wat er gebeurt.
   Bij de opbouw van p5 worden de gewogen basisvectoren als lijnen geplot. Geef aan de hand van het ontstane patroon aan hoe je met de beschikbare basisvectoren, die geen van allen een bijzonder voorkeur hebben voor een puls op positie 5, toch een scherpe puls op die positie kunt krijgen terwijl alle andere posities nul zijn.

Met behulp van de DFT is elke vector uit sinus- en cosinusbasisvectoren op te bouwen.

x = 2*(rand(1,nfft)-0.5); % random vector met waarden tussen -1 en 1
X = x*B; % FT van randomvector
L = nfft;
for n = 1:nfft
  clf; hold on
  plot(1:L, x, '*r'); axis([0 L+1 -1 1]) % Target
  stem(1:L, X(1:n)*B(:,1:n)'); % Benadering met basisvectoren
  axis([0 L+1 -1 1]); hold off
  title(['Construction of random signal: n = ' int2str(n)])
  pause(10./nfft)
end

opdracht D

1. Voer deze FOR-lus net als bij de vorige vraag stapgewijs uit. Leg uit waarom de benadering van het random signaal als het effect van een (steeds minder) laag-doorlaat filter beschouwd kan worden.
2. Edit (in een editor) bovenstaande FOR-lus zodanig dat het effect van een (steeds minder) hoog-doorlaat filter wordt benaderd.

Complexe DFT - de FFT

De complexe discrete FT, uitgevoerd met het Fast Fourier Tansformatie, is veel efficiënter om uit te rekenen en eigenlijk ook intuïtiever in het gebruik dan de DFT op basis van de matrixberekening zoals tot nu toe. Er is een aantal verschillen met boven: (Kijk hier.)

• Per golfgetal k zitten de cosinus en sinus bijdragen in het zelfde complexe getal, repectievelijk als reëel en imaginair deel.
• Het golfgetal en het indexnummer in het frequentiedomein verhouden zich als k = n+1 voor k < N/2.
• In een strict reëel tijddomeinsignaal geldt voor n > N/2 dat dit de complex geconjugeerden zijn van de rond n = N gespiegelde waarden.
• Er is een andere normalisatie: in het energiedomein is de energie een factor N te groot. Wanneer de energie van x gelijk aan is aan 1 dan is de energie van X = fft(x) gelijk aan sum(X.*conj(X)) = N, maar de energie van ifft(fft(x)) is weer gelijk aan aan 1. De inverse transformatie maakt dit weer goed. (In verband met numerieke efficiëntie worden de sqrt berekening en de delingen voorkomen.)

Voer help fft uit en lees de bijbehorende uitleg.

opdracht E

1. Geef aan hoe je in de definitie van de FFT en de IFFT de inproducten (d.w.z. een som van producten) kunt vinden.
2. Leg uit dat hoewel de normalisering iets anders is dan de eisen voor basisvectoren voorschrijven dat toch geldt dat x = ifft(fft(x)).

de efficiëntie van de FFT
De belangrijkste reden voor de voorkeur voor het toepassen van de FFT boven de directe berekening op basis van het afbeelden op de basisvectoren is dat de FFT veel efficiënter berekend wordt. De FFT is efficiënter doordat optimaal gebruik wordt gemaakt van het feit dat de basisvectoren, als sinusoiden, zoveel symmetrieën vertonen. Daarom wordt maar een minimaal aantal sinus- of cosinuswaarden berekend, de rest volgt uit symmetrieoverwegingen. De FFT is maximaal efficiënt wanneer de inputvector steeds in tweeën kan worden verdeeld. Dit is het meest het geval bij machten van 2. Daarom worden FFT's in de praktijk ook meestal op vectoren die een macht van twee zijn losgelaten.

Het efficiëntie verschil tussen de directe berekening en de FFT kan worden bepaald met de volgende code:

m = 10
nfft = 2.^m % Verzameling nfft-dimensionale basisvectoren
B = zeros(nfft,nfft); k=0;
B(:,1) = cos(2*pi*k*(0:nfft-1)/nfft)' / sqrt(nfft);
for k = 1:nfft/2-1
  B(:,k*2) = sin(2*pi*k*(0:nfft-1)/nfft)' / sqrt(nfft/2);
  B(:,k*2+1) = cos(2*pi*k*(0:nfft-1)/nfft)' / sqrt(nfft/2);
end
B(:,nfft) = cos(2*pi*nfft/2*(0:nfft-1)/nfft)' / sqrt(nfft);

p = rand(1,nfft);
t1 = cputime; for ii = 1:1000; p*B; end; t1 = cputime-t1
t2 = cputime; for ii = 1:1000; fft(p); end; t2 = cputime-t2

opdracht F

1. Welke factor (voor m = 10) is de FFT sneller dan de directe berekening? Voer de berekening meer dan 1 keer uit.
2. Bepaal het verschil in berekeningssnelheid voor een priemgetal in de buurt van 2^10. Welke conclusie trek je over het gebruik van de FFT?

Impulsrespons van de (I)FFT

De impuls-respons van een lineair systeem (in dit geval van de FT) is van groot belang omdat het (bijna) alles zegt over het systeem. Elke input kan immers beschreven worden door een opeenvolging van pulsen, wanneer je dus weet hoe het systeem reageert op een puls, kun je door een optelling van de reactie op elk van een opeenvolging van pulsen de reactie op een willekeurig input berekenen. Laten we een kijken wat er gebeurt wanneer de de puls van positie 1 tot positie nfft opschuiven.

nfft = 32;
for t = 1:nfft
  p = zeros(1,nfft); p(t) = 1; % puls op vectorpositie 1
  P = fft(p);
  Pcos = real(P);
  Psin = imag(P);
  stem(p); hold on
  plot(1:nfft, Pcos, '-r*', 1:nfft, Psin, '-g*'); hold off
  title(['Reele en imaginaire deel van puls op positie ' int2str(t)]);
  axis([0 nfft+1 -1.1 1.1]);
  xlabel('Basisvector nummer');
  pause % Ga door na toetsaanslag
end

opdracht G

1. Zowel de puls zelf als de fouriergetransformeerde ervan zijn in de figuur weergegeven. Leg van beide representaties uit hoe de x-as moet worden geïnterpreteerd.
2. Leg uit dat om een puls op positie 1 te maken, uitsluitend cosinus-bijdragen nodig zij n.
3. Leg uit hoe je uit de figuur kunt afleiden dat voor n > N/2 geldt dat dit de complex geconjugeerden zijn van de rond n = N gespiegelde waarden zijn. Waarom hebben de cosinusbijdragen dit niet?

Windows (vensters)

Windows spelen vooral een rol wanneer je uit een langer signaal, een kleiner deel wilt selecteren om er een FT-analyse op uit te voeren. Laad eerst een stukje spraaksignaal.

[x, fs] = wavread('nul123.wav');
soundsc(x, fs) % Afspelen geluidssignaal x

opdracht H

1. Met welk tijdstip correspondeert sample 2000? Geef de formule.
2. Hoe lang is het totale signaal? (Hint: gebruik de Matlab-functie size.)

Selecteer en window een deel van x vanaf sample 2001 tot 2000 + nfft:

nfft = 256;
y = x(2000+[1:nfft]);
yham = y.*winHamming(nfft);
yhann = y.*winHanning(nfft);

subplot(2,1,1)
plot((1:length(y))/fs*1000,y,(1:length(y))/fs*1000,yham, (1:length(y))/fs*1000,yhann);
axis tight; grid on
xlabel('Time in ms')

opdracht I

1. Bepaal uit dit tijddomeinplaatje de twee meest dominante frequenties. Gebruik de grafiek van het rechthoekig window.
2. Leg aan de hand van de vorm van het Hamming en het Hanning window uit waarom er rond de randen van het window kleine verschillen bestaan tussen de gewindowde signalen.

Bereken en plot de spectra in dB:

subplot(2,1,2)
f_As = 0:fs/nfft:fs/2;
Y = 10*log10(abs(fft(y))); Y = Y(1:nfft/2+1);
Yham = 10*log10(abs(fft(yham))); Yham = Yham(1:nfft/2+1);
Yhann = 10*log10(abs(fft(yhann))); Yhann = Yhann(1:nfft/2+1);
plot(f_As,Y,f_As,Yham,f_As,Yhann);
grid on;
xlabel('Frequentie in Hz'); ylabel('Energie in dB')

Check of je tijddomeinberekening van de twee belangrijkste frequenties klopt.

opdracht J

1. Het verschil tussen de toepassing van het Hamming en het Hanning window enerzijds en het rechthoekige window anderzijds is bijzonder groot. Wat heeft de voorkeur? Waarom?

Logsweep

Als voorbeeldsignaal nemen we nu een logsweep, dat is een signaal dat in frequentie verandert op een wijze die op een logaritmische as lineair is. Het eigenschap van een logsweep (die een lineaire sweep, met de functie linsweep.m, niet heeft) is dat elke octaaf verandering (een factor twee in frequentie) evenlang duurt.

In dit geval gaan we in 3 seconden van 50 Hz naar 3500 Hz. Dat komt overeen met 6.1 octaaf in 3 seconden. Dit komt dus overeen met ongeveer 2 octaaf per seconde (een octaaf is een factor 2).

opdracht K

1. Toon bovenstaande bewering aan met een berekening.

Toets help logsweep en lees de tekst. Het enige bijzondere is FRATE, dit is de inverse (uitgedrukt in Hz) van stepT (in ms) zoals die door fftGram wordt gebruikt. FRATE bepaalt het aantal frequentiewaarden f dat per seconde wordt gegenereerd. Hiermee krijgen we voor elke kolom in het spectrogram een bijbehorende waarde van f.

startf = 50; % Hz
endf = 3500; % Hz
duration = 3; % sec
fs = 8000; % Hz
offsetPhase = 0; % [0,1], 0 is sinus-fase, 0.25 = cosinus-fase
FRATE = 200; % ms, frame rate
[x,f] = logsweep(startf, endf, duration, fs, offsetPhase, FRATE);
t = linspace(0,duration,length(f));

We hebben nu een signaal x (samplefrequentie fs = 8000 Hz) met een frequentieontwikkeling f met bijbehorende tijdstippen t (gesampled op 200 Hz, dus elke 5 ms). Het "geluidssignaal" x is af te spelen met sound(x,fs).
Merk op dat je dit als een continu en gelijkmatig geheel waarneemt. En te plotten met:

clf;
tx = linspace(0,3,length(x));
plot(tx, x);

Dit geeft een slecht interpreteerbaar plaatje.

opdracht L

1. Zoom in op een aantal delen van het signaal, bijvoorbeeld het gebied tussen t = 0.5 s en t = 0.7 s via

   tStart = 0.5; tEnd = 0.7;
   sampStart = fs*tStart;
   sampEnd = fs*tEnd;
   plot(tx(sampStart:sampEnd), x(sampStart:sampEnd));
   axis tight; grid on

   Hoe kan je zien dat de frequentie verandert? Zoom ook in op een interval vlak onder t = 3 s.

Een andere (handiger) manier om in te zoomen is met gebruik van de functie strips (zie help strips) en het gebruik van de muis om in te zoomen: (Dit werk niet zonder Signal Processing Toolbox!) strips(x, .5, fs) % Geeft 6 strips van elk 0.5 sec

opdracht M

1. Na t = 1.5 s zie je waarschijnlijk wat structuur die verdwijnt wanneer je er op inzoomt. In de onderste strip vindt je structuren waar dat niet voor geldt. Geef voor beide verschijnselen een verklaring.

Tot slot van dit onderdeel nog even check van de instantane frequentieontwikkeling zich gedraagt zoals dat zou moeten.

subplot(2,1,1)
plot(t, f); title('Lineaire frequentie-as'); grid
subplot(2,1,2)
semilogy(t, f); title('Logaritmische frequentie-as'); grid

Voor een logaritme geldt dat een vaste afstand op de x-as correspondeert met een bepaalde factor. Controleer dit door een paar punten van de bovenste subplot te nemen.

Spectogram

We gaan ons nu richten op het visualiseren van de frequentieontwikkeling als functie van de tijd in de vorm van spectrogrammen. We gebruiken eerst een rechthoekig window.

stepT = 1000/FRATE; % factor 1000 om ms om te zetten in sec
nfft = 256;
DynRange = 60; % plot met 60 dB dynamische bereik.
clf
[S SInfo] = fftGram(x, fs, nfft, stepT, 'logsweep with rectangular window', DynRange);
colorbar

In ieder geval is het een mooi kleurig plaatje. Het is echter de vraag is of het ook een goed plaatje is. Volgens de tijddomeinanalyse kan het signaal op elk moment beschreven worden met een een enkele frequentie, en zeker niet met steeds in sterkte wisselende bijdragen van vrijwel alle andere frequenties (en dus alle andere FT-basisvectoren). Deze bijdragen worden ook wel "lek" genoemd.

opdracht N

1. Beschrijf de spectrale weergave van de logsweep.
2. Voeg een plot van de frequentieontwikkeling toe met behulp van

   hold on; plot(t,f,'w'); hold off

   Zoom in op de de linker onderkant. Wat concludeer je over de bijdrage van de lek als functie van de instantane frequentie f?
3. Maak een nieuw figuur met de functie figure en gebruik de zoom functie hierboven om een enkel window te visualiseren. Kies een tStart en maak gebruik van tEnd = tStart + nfft/fs; Houdt er rekening mee dat de tijd-as correspondeert met het midden van een window.
   Vergelijk in twee subplots windows met veel en weinig lek. Wat kun je zeggen over grafieken met veel en weinig lek? (Hint: dit probleem wordt gereduceerd door de globale vorm van het Hamming en het Hanning window).

Laten we nu het effect van de al eerder geïntroduceerde windows eens analyseren.

DynRange = 120;
clf;
subplot(2,2,1)
[Srect SInfo] = fftGram(x, fs, nfft, stepT, 'rect', DynRange);
subplot(2,2,2)
[Striang SInfo] = fftGram(x, fs, nfft, stepT, 'triang', DynRange);
subplot(2,2,3)
[Shamming SInfo] = fftGram(x, fs, nfft, stepT, 'hamming', DynRange);
subplot(2,2,4)
[Shanning SInfo] = fftGram(x, fs, nfft, stepT, 'hanning', DynRange);

opdracht O

1. Beschrijf, kwalitatief, de lek van elk window in termen van gemiddeld niveau, variatie met de instantane frequentie, spectrale (vertikale) fijnstructuur en spectrale helling.

We maken een aparte figuur en zoomen in tussen t = 0.6 s en t = 1.2 s en f = 0 Hz tot f = 400 Hz.

figure
subplot(2,2,1)
imagesc(SInfo.T, SInfo.F, Srect./max(Srect(:)));
set(gca, 'ydir', 'normal');
ylabel('Freq in Hz'); title('rect')
hold on; plot(t,f,'w'); hold off
axis([.6 1.2 0 800]);
subplot(2,2,2)
imagesc(SInfo.T, SInfo.F, Striang./max(Striang(:)));
set(gca, 'ydir', 'normal');
ylabel('Freq in Hz'); title('triang')
hold on; plot(t,f,'w'); hold off
axis([.6 1.2 0 800])
subplot(2,2,3)
imagesc(SInfo.T, SInfo.F, Shamming./max(Shamming(:)));
set(gca, 'ydir', 'normal');
xlabel('Time in sec'); ylabel('Freq in Hz'); title('hamming')
hold on; plot(t,f,'w'); hold off
axis([.6 1.2 0 800])
subplot(2,2,4)
imagesc(SInfo.T, SInfo.F, Shanning./max(Shanning(:)));
set(gca, 'ydir', 'normal');
xlabel('Time in sec'); ylabel('Freq in Hz'); title('hanning')
hold on; plot(t,f,'w'); hold off
axis([.6 1.2 0 800])

opdracht P

1. Beschrijf van elke window het gedrag en de vorm van de hoofdpiek in termen van breedte, spectrale vorm en variatie met de instantane frequentie.
2. Welke piek is het scherpst en waarom is dat?
3. Welke is het breedst?
4. Welk window heeft je voorkeur, en waar is deze voorkeur van afhankelijk?

twee componenten
We gaan de situatie nu een beetje realistischer maken door aan de ene geluidscomponent een tweede, genaamd x2, toe te voegen. Dit is ook een logsweep, maar nu precies in tijd omgedraaid (dat was via x2 = fliplr(x); gemakkelijker geweest dan het signaal overnieuw maken zoals hieronder gebeurt). De tweede logsweep x2 loopt van f = 3500 Hz naar f = 50 Hz.

startf = 3500; % Hz
endf = 50; % Hz
duration = 3; % sec
fs = 8000; % Hz
offsetPhase = 0; % [0,1], 0 is sinus-fase, 0.25 = cosinus-fase
FRATE = 200; % ms
[x2,f2] = logsweep(startf, endf, duration, fs, offsetPhase, FRATE);

soundsc(x2, fs) % plays logsweep

Om het nog wat interessanter te maken gaan we x2 een factor 1000 zwakker maken:

soundsc(x+x2/1000,fs) % plays x en x2

opdracht Q

1. Kun je de tweede component horen?
2. Bereken het aantal dB in het energiedomein waarmee een factor 1000 in intensiteit in het amplitudedomein correspondeert.

Van de gewogen optelling gaan we weer een spectrogram maken. Eerst weer een rechthoekig window:

stepT = 1000/FRATE; % ms
nfft = 256;
DynRange = 120;
clf % clear figure
[S SInfo] = fftGram(x+x2/1000, fs, nfft, stepT, 'rect', DynRange);

Het is duidelijk dat er (nagenoeg) niets te zien is van een downsweep in deze spectrale weergave. Laten we nu de andere windows vergelijken.

subplot(2,2,1)
[Srect SInfo] = fftGram(x+x2/1000, fs, nfft, stepT, 'rect', DynRange);
subplot(2,2,2)
[Striang SInfo] = fftGram(x+x2/1000, fs, nfft, stepT, 'triang', DynRange);
subplot(2,2,3)
[Shamming SInfo] = fftGram(x+x2/1000, fs, nfft, stepT, 'hamming', DynRange);
subplot(2,2,4)
[Shanning SInfo] = fftGram(x+x2/1000, fs, nfft, stepT, 'hanning', DynRange);

opdracht R

1. Verklaar de verschillen aan de hand van de windoweigenschappen die je al in kaart hebt gebracht.

Wanneer we inzoomen tussen t = 0.6 s en t = 1.2 s en f = 0 Hz tot f = 800 Hz, kunnen we de verschillen nog wat beter zien.

subplot(2,2,1)
imagesc(SInfo.T, SInfo.F, Srect);
set(gca, 'ydir', 'normal');
xlabel('Time in sec');
ylabel('Freq in Hz'); title('rect')
axis([.5 2.5 0 800])
subplot(2,2,2)
imagesc(SInfo.T, SInfo.F, Striang);
set(gca, 'ydir', 'normal');
xlabel('Time in sec');
ylabel('Freq in Hz'); title('triang')
axis([.5 2.5 0 800])
subplot(2,2,3)
imagesc(SInfo.T, SInfo.F, Shamming);
set(gca, 'ydir', 'normal');
xlabel('Time in sec');
ylabel('Freq in Hz'); title('hamming')
axis([.5 2.5 0 800])
subplot(2,2,4)
imagesc(SInfo.T, SInfo.F, Shanning);
set(gca, 'ydir', 'normal');
xlabel('Time in sec');
ylabel('Freq in Hz'); title('hanning')
axis([.5 2.5 0 800])

opdracht S

1. Welk window zou je voor welk doel willen gebruiken?
   1. het scheiden van verschillende componenten,
   2. het beschrijven van een breedbandig spectrum,
   3. het zo nauwkeurig mogelijk schatten van lokale instantane frequentie van een signaalcomponent.

MFCC (mell-scaled cepstral coefficients

MFCC's zijn de standaard input voor ASR systemen. Kernpunt van de MFCC's is met zo weinig mogelijk getallen de vorm van de spectrale omhullende (waarin de formantinformatie besloten zit) weer te geven. Bij MFCC's wordt niet gestreefd naar het vaststellen van individuele formanten (wat in het algemeen erg moeilijk is), maar naar een representatie die de zelfde informatie vertegenwoordigt, zonder dat je de kans loopt om fouten te maken in de schatting van individuele formanten.

nfft = 128;
stepT = 10;
[x0123, fs] = wavread('nul123.wav');
clf;
[S SInfo] = fftGram(x0123, fs, nfft, stepT, 'hamming', 60);

opdracht T

1. Schat van de 'nul', de 'een' en de 'drie' formantposities en formantbreedtes (in Hz) in een drietal frames.

De MFCC wordt berekend door in elk frame het volgende recept uit te voeren: (in Matlab mfcc.m)

• frameblocking (ophakken in overlappende stukjes van ongeveer 25 ms, elk 10 ms verschoven)
• windowing
• FFT
• kwadrateren of absolute waarde nemen voor "energiespectrum"
• integreren onder 30 tot 40 Mell-scaled triangular filters (de Mell-scale is lineair onder de 1000 Hz en logaritmisch erboven)
• terugtransformeren met een (inverse) cosinustranformatie (dit geeft een zogenaamd cepstrum)
• alleen de 10 tot 14 laagfrequente cepstrale componenten behouden: de MFCC's

Als laatste stap kun je de informatie die je overhoudt weer terugtransformeren met een cosinustransformatie om een visuele representatie van het "spectrum" dat de MFCC's representeren te maken. De frequentie-as van dit spectrum komt overeen met de centrumfrequenties van de driehoekige Mell-filterbank.

Doe nu:

[ceps,freqresp,fb,fbrecon,freqrecon, freqs] = mfcc(x0123, fs, 100);
subplot(2,1,1)
imagesc(ceps); axis xy % De MFCC's
xlabel('Tijd in frames van 10 ms')
ylabel('Cosinus-basisvector #')
subplot(2,1,2);
imagesc(fbrecon); axis xy % Een visuele reconstructie hiervan
xlabel('Tijd in frames van 10 ms')
ylabel('Filterbank #')
freqs(2:end) % Centrum frequenties van de filterbank

De laagste component van de MFCC's (de energie) wijkt vaak af in amplitude. Je krijgt een beter zicht op bijdragen van de cosinusbasisvectoren wanneer je deze component weglaat in de visualisatie.

subplot(2,1,1)
imagesc(3:size(ceps,2), 1:size(ceps,1), ceps(2:end,:)); axis xy
xlabel('Tijd in frames van 10 ms')
ylabel('Cosinus-basisvector #')

opdracht U

1. Bepaal de formant-frequenties uit de reconstructie (gebruik de lijst freqs om van filterbanknummers frequenties te maken) en vergelijk die met de waarden die je zelf had berekend. Beschrijf wat je opvalt.
2. Geef, op basis van de informatie over de cosinus tranformatie, aan waarom de bijdrage van de twee component in de figuur (dus de derde basisvector in werkelijkheid) zo sterk is bij de laatste 4 woorden.