Luento 25.2: Suodinsuunnittelu

Tunnin aluksi tarkasteltiin demoa, jossa kompleksinen taajuusvaste esitetään painotettujen kompleksisten eksponenttifunktioiden summana. Näin havaittiin mm. ettei pienellä määrällä kertoimia voida muodostaa kovin monimutkaisia amplitudivasteita. Alla olevassa kuvassa suotimessa on 9 kerrointa, jolloin taajuusvaste saadaan yhdeksän kompleksiluvun summana. Amplitudivaste on oikean alakulman käyrä, joka on siis vasemman alakulman käyrän etäisyys origosta.

Tämän jälkeen tutustuttiin suodinsuunnitteluun ikkunamenetelmällä. Suunnittelukriteerit ovat kahtalaiset: suotimen taajuusvasteen määräämiseksi pitää tietää millainen vaihevaste halutaan ja millainen amplitudivaste halutaan.

Vaihevasteen osalta vaaditaan että kaikkien taajuuksien tulee viivästyä yhtä paljon. Tämä toteutuu jos vaihevaste on lineaarinen. Yksinkertaisimmissa tapauksissa vaihevasteen lauseke voi olla siis esimerkiksi muotoa -2w, joka taatusti on lineaarinen. Matlabissa tällainen kuvaaja saadaan esim. komennoilla:

>> [H,W] = freqz([1, 1, 1]);
>> plot(angle(H));
>> grid on

Freqz-funktiosta saa siis ulos taajuusvastefunktion arvoja vektorissa H. Vektorissa on lueteltu taajuusvasteen kompleksiset lukuarvon 512:ssa pisteessä taajuusakselilla. 

Vaihevasteen derivaatasta käyteään nimeä ryhmäviive, ja se ilmaisee suoraan eri taajuuksille tulevan viiveen näytteinä (miinusmerkkisenä). Lopuksi todettiin, että vaihevaste on aina lineaarinen, jos impulssivasteen termit ovat symmetrisesti keskipisteen suhteen.

Amplitudivasteen osalta tavoitteena on saada vaste päästökaistalla ykköseksi ja estokaistalla nollaksi. Käytännössä tämä ei ole mahdollista, vaan suotimelle täytyy antaa hieman toleranssia ja sallia tietty määrä värähtelyä molemmilla kaistoilla. Lisäksi kaistojen väliin täytyy sallia "don't care" -alue, jossa amplitudivaste saa olla mitä vain.

Prujussa ratkaistaan mikä impulssivaste toteuttaisi ideaalisen amplitudivasteen (arvot vain nollaa tai ykköstä). Osoittautuu että impulssivasteen muoto on tuttu sinc-funktio, mutta sen pituus on ääretön. Tämän vuoksi suotimesta ei saataisi ainuttakaan vastearvoa koskaan, vaan laskentaa tarvittaisiin äärettömän paljon.

Tästä ongelmasta päästään katkaisemalla impulssivaste, mutta tämä luonnollisesti vaikuttaa amplitudivasteeseen. Oikealla olevan kuvan mukaisen demottiin, että suoralla katkaisulla ei estokaistan värähtelyä saada millään alle n. 21 desibelin, ja päästökaistallakin suurin heitto on luokkaa 0.7 dB. Ratkaisu tähän on käyttää ikkunointia, eli kertoa katkaistu impulssivaste jollain ikkunafunktiolla. Näin voidaan päästä parempiin vaimennusominaisuuksiin.

Ideaalisen suotimen impulssivasteen pituus on ääretön, eikä sitä voi käytännössä toteuttaa. Näin ollen impulssivaste on katkaistava, mistä seuraa vääristymä amplitudivasteeseen. Matlab-testeillä havaittiin, että tätä ei voi kompensoida esim. kertoimia lisäämällä, vaan on käytettävä ikkunaa, joka pehmentää katkaisun vaikutusta. Ikkunoita on lueteltu esim. sivun 84 taulukossa, ja mitä paremmat vaimennusominaisuudet niillä on, sitä leveämpi siirtymakaistasta tulee. Onneksi tätä voidaan kuitenkin kompensoida kertoimia lisäämällä.

Luento 18.2: Z-muunnoksen laskenta ja sovellukset

Tänään käsiteltiin suotimen analysointi (kappale 4) loppuun.

Menetelmässähän ratkaistaan ensin impulssivaste, sitten siirtofunktio ja lopuksi taajuusvaste. Taajuusvaste on kompleksifunktio, joten sitä ei voida sellaisenaan piirtää 2-ulotteiseen koordinaatistoon. Näin ollen piirretään kaksi kuvaajaa: funktion itseisarvon kuvaaja sekä sen vaihekulman kuvaaja. Näistä edellinen kertoo kuinka paljon eri taajuuksien amplitudit muuttuvat suodatuksessa ja jälkimmäinen paljonko ne viivästyvät suodatuksessa. Amplitudivaste on näistä mielenkiintoisempi, koska sen avulla taajuuksia saadaan esim. poistettua yksinkertaisesti huolehtimalla että amplitudivaste ko. taajuudella on nolla.Vaihevaste puolestaan kertoo paljonko eri taajuudet viivästyvät suodatettaessa.

Amplitudivastetta tarkasteltaessa on kätevämpi käyttää desibeliasteikkoa, joka on logaritminen. Logaritmi tekee kertolaskusta yhteenlaskua, ja korostaa lähellä nollaa olevia eroja, jotka molemmat ovat meille käteviä ominaisuuksia.

Kappaleen 4 kaksi viimeistä lukua käytiin läpi ratkaisemalla kevään 2013 toukokuun tentin tehtävä 4. Tässä tehtävässä on annettu suotimen yhtälö, josta täytyy ratkaista siirtofunktio, piirtää napa-nollakuvio sekä päätellä stabiilisuus.

Lopuksi tutkittiin napa-nollakuvion ja taajuusvasteen suhdetta. Taajuusvastehan saadaan siirtofunktiosta H(z) evaluoimalla se pisteissä z = exp(iw). Geometrisesti tämä tarkoittaa yksikköympyrän reaaliakselin yläpuolella olevia pisteitä. Jokainen napa-nollakuvion nolla laskee taajuusvastetta ja jokainen napa nostaa taajuusvastetta. Tästä nähtiin alla olevan kuvan mukainen demo, jossa hiirellä voidaan sijoitella napoja ja nollia yksikköympyrälle. Alimpaan kuvaan piirtyy jokaisen klikkauksen jälkeen suorimen amplitudi- ja taajuusvasteet.

Luento 11.2: Z-muunnos

Ensimmäisen tunnin aluksi luotiin katsaus Fourier-muunnoksen ja sen yleistysten soveltamiseen koneoppimisessa. Fourier-analyysissähän kysytään kuinka paljon kutakin sinisignaalia on mukana tarkasteltavassa signaalissa. Yleisempi muoto on käyttää jotain muuta signaalikokoelmaa, tai oppia tämä kokoelma datasta. Klassiset menetelmät ovat pääkomponenttianalyysi (PCA) tai Helsingissä kehitetty riippumattomien komponenttien analyysi (ICA), joissa signaalit esitetään sellaisten rakennuspalikoiden avulla että suuria kertoimia tulee mahdollisimman vähän.

Käytimme erästä tällaista hajotelmaa (SPAMS dictionary learning) osallistuessamme Kaggle-alustalla organisoituun linnunlauluntunnistuskilpailuun. Alkuvaiheessa käytimme mm. Fourier-muunnosta, mutta pelkkä taajuuksien analyysi ei tuottanut tulosta. Tämän vuoksi päätimme oppia "sanakirjan" suoraan datasta, ja toivoimme että sanakirjaan päätyisi tyypillisiä eri lintulajien viserryksiä. Näiden viserrysten lukumäärä toimi sitten indikaattorina siitä mitä lintulajeja äänityksessä oli.

Toisena esimerkkinä koneoppimiskilpailusta tarkasteltiin viime kesänä ollutta MEG-aivodatan analyysikilpailua, johon osallistuimme laitoksen kesätyöntekijöiden kanssa.

Toisella tunnilla tarkasteltiin Z-muunnosta ja sen tärkeimpiä ominaisuuksia. Z-muunnoksen avulla voidaan selvittää mm. suotimen stabiilisuus: suodin on stabiili jos kaikki siirtofunktion navat ovat yksikköympyrän sisäpuolella.

Luento 4.2: FFT-algoritmi

Tänään tarkasteltiin Fourier muunnoksen ominaisuuksia, sovelluksia sekä nopeaa toteutusta.

Luennon aluksi esiteltiin alkeellinen menetelmä puheen tunnistukseen. Kirjan Elements of statistical learning kappaleessa 5.2.3 opetetaan tietokone erottelemaan kaksi vokaalia niiden Fourier-muunnosten perusteella. Menetelmä on nimeltään logistinen regressio, joka monimutkaisista kaavoista huolimatta on varsin yksinkertainen toteuttaa: menetelmä etsii kertoimet kullekin Fourier-muunnoksen taajuudelle, ja laskee tulokset yhteen. Jos luku on positiivinen, tulkitaan äänne ä-kirjaimeksi, muutoin a-kirjaimeksi.

Menetelmää demottiin Matlab-toteutuksella, jossa luokittelija opetettiin erottamaan S-äänne muista äänteistä näyttämällä luokittimelle esimerkkejä kahteen luokkaan kuuluvista Fourier-muunnoksista. Menetelmä toimi kohtalaisen hyvin ottaen huomioon opetusaineiston erot testiaineistoon.

Esimerkki kuvaa hyvin tämän päivän signaalinkäsittelyalgoritmia: perusmenetelmiä (Fourier-muunnos, konvoluutio, jne.) käytetään piirregeneraattoreina, jotka tuottavat hieman parempaa raakadataa kuin suora mittaussignaali (esim. taajuustietoa eikä raakaa mittausdataa). Laskettujen piirteiden perusteella sitten nostetaan tiedon abstraktiotasoa edelleen. Esimerkiksi äänteen tunnistuksessa hierarkia on esimerkiksi seuraava:

16000 aikatason näytettä -> 128 taajuustason kerrointa -> 1 bitti, joka kertoo kumpi äänne on kyseessä

Toisena  esimerkkinä oppivasta järjestelmästä mainittiin tamperelaisen Visy Oy:n rekisterikilven tunnistus: tällöinkin luokittelijalle on näytetty kymmeniä tuhansia käsin kerättyjä kirjainmerkkejä, ja algoritmi on oppinut erottelemaan ne toisistaan.

Tämän jälkeen siirryttiin tarkastelemaan Fourier-muunnoksen ominaisuuksia. Ominaisuuksista tutustuttiin lähemmin siirtoon ajassa (esim. laske signaalin x(n+20) muunnos, kun tiedetään x(n):n muunnos) sekä konvoluution muunnokseen (DFT muuntaa konvoluution kertolaskuksi, eli x(n)*y(n) -> X(n)Y(n)). Tämä on perustana mm. dekonvoluutiolle joka on konvoluutiolle käänteinen operaatio. Menetelmää käytettiin mm. Hubble-teleskoopin alkuaikoina, jolloin yhdessä peilissä olleen hiontavirheen vuoksikuvat olivat sumuisia. Kuvantamisprosessia voidaan nimittäin mallintaa (kaksilulotteisella) konvoluutiolla

y(n,m) = h(n,m) * x(n,m),

missä x on todellinen näkymä, y on havaittu sumuinen kuva, ja h on linssin impulssivaste (nk. point spread function; PSF). Yhtälössä y ja h ovat tunnettuja, ja tehtävänä on ratkaista x. Ratkaisu löytyy taajuustasossa, koska

Y(n,m) = H(n,m) X(n,m),

joten (Matlabin syntaksilla ilmaistuna):

x(n,m) = ifft (Y(n,m) ./ H(n,m)).

Dekonvoluutiosta on hyötyä yleisemminkin lineaarisen kanavan aiheuttaman häiriön poistossa. Jos tiedetään signaalin x kulkeneen kanavan h läpi, voidaan vastaanotetusta mittaustuloksesta ypäätellä x, jos meillä on joku käsitys kanavasta h. Esimerkkinä tästä on esim. langattoman tiedonsiirtokanavan estimointi ja sen aiheuttaman vääristymän kompensointi.

Toinen menetelmän tuottama etu on että Fourier-muunnoksen (käytännössä FFT:n) avulla voidaan laskea konvoluutio kaavasta (Matlabin syntaksilla ilmaistuna):

conv(x,y) = ifft(fft(x) .* fft(y))

Lisäksi käsiteltiin nopeaa Fourier-muunnosta eli FFT:tä, joka on vain nopeampi tapa toteuttaa diskreetti Fourier-muunnos (DFT). FFT perustuu signaalin jakamiseen lyhyempiin pätkiin, jotka muunnetaan jakamalla ne edelleen rekursiivisesti kahtia. Rekursio päättyy, kun muunnoksen pituus on 1, jolloin muunnosta ei tarvitse enää tehdä. 1-ulotteisen vektorin tapauksessa muunnosmatriisi on yksinkertaisesti F = [1], joka tarkoittaa pelkkää ykkösellä kertomista eikä sitä tarvitse tehdä. Lyhyemmistä vektoreista saadaan koostettua pidemmät vektorit kaavoilla (3.3) ja (3.4).

Alla on vielä luennon esimerkkikoodi S-kirjaimen tunnistuksesta. Funktio käyttää Stanfordin yliopistossa kehitettyä glmnet-pakettia.

function vokaalin_tunnistus()
%
% Esimerkki vokaalin ja S-kirjaimen erottelusta äänisignaalista.
% heikki.huttunen@tut.fi -- 4.2.2015
%

close all

addpath /home/hehu/Documents/Libraries/glmnet_matlab/

% Ladataan opetusaineisto:

[x, Fs] = audioread('seiska.wav');

[X, H, numFrames] = extractFeatures(x, Fs);
title ('Merkitse S-kirjaimet hiirella');

isConsonant = zeros(numFrames, 1);

while true
    
    [x1, y1] = ginput(1);
    [x2, y2] = ginput(1);
    
    if x1 > x2
        xt = x1;
        x1 = x2;
        x2 = xt;
    end
    
    isConsonant(round(x1 * numFrames) : round(x2 * numFrames)) = 1;
    
    response = questdlg('Jatketaanko annotointia?', ...
        'Kysymys', ...
        'Kyllä', 'Ei', 'Kyllä');
    
    if strcmp(response, 'Ei')
        break
    end
    
end

cvob2 = cvglmnet(X, isConsonant, 'binomial', [], 'class');
yHat = glmnetPredict(cvob2.glmnet_fit, X, cvob2.lambda_min, 'response');

coefficients = H * cvob2.glmnet_fit.beta(:, cvob2.lambda == cvob2.lambda_min);

figure()
subplot(211)
plot(yHat);
ylabel('S-kirjaimen TN')
subplot(212)
stem(coefficients)

response = questdlg('Valmiina tunnistamaan?', ...
        'Tunnistus', ...
        'OK', 'OK');

while true
    
    close all

    myRecObj = audiorecorder(Fs, 16, 1);
    recordblocking(myRecObj, 2);
    y = getaudiodata(myRecObj);

    X = extractFeatures(y, Fs);
    yHat = glmnetPredict(cvob2.glmnet_fit, X, cvob2.lambda_min, 'response');

    figure()
    plot(yHat);

    response = questdlg('Jatketaanko tunnistusta?', ...
        'Kysymys', ...
        'Kyllä', 'Ei', 'Kyllä');
    
    if strcmp(response, 'Ei')
        break
    end
    
end

end

function [F, H, numFrames] = extractFeatures(x, Fs)

    [~,f,t,S] = spectrogram(x, 256, 128, 256, Fs, 'yaxis');
    surf(t, f, 10*log10(abs(S)), 'EdgeColor', 'none');
    axis xy; axis tight; colormap(jet); view(0,90);

    S = log10(S)';
    
    H = [];
    n = (1:size(S, 2))';
    
    for k = 0:3
        H = [H, n.^k];
    end
    
    F = S * H;
    numFrames = size(S, 1);
    
end