R12-kommentteja

Viimeiset Matlab-demot laskennan epätarkkuuksista ja sekvenssin digitaalisesta uudelleennäytteistämisestä.

Lisäkokeiluna (viimeine kuva) miltä kiisseli kuulostaa alasnäytteisettynä, jos ei desimointisuodatusta.

R11 dsp:n estetiikkaa

Mitä 70 minuutissa tulee sanottua? Mitkä on tavoitteet ja miten niihin päästiin? Erityisesti, mitkä ovat olleet implisiittiset esitiedot FFT:n osalta?

Tiistain suomi-englanti-istunnossa aika meni kivasti yhteen tehtävään. Keskiviikkona yhdellä kielellä aikaa meni vajaa tunti eli kielijuttuun menee noin +20% aikaa tason samalla laskiessa.

Kysymyksiä tuli juuri noihin "epäselviin" kohtiin eli perhosyhtälöparin alfaan ja betaan (virtauskaavion rivit) ja muuttujaan l, joka juoksee perhosten mukaan ("yksi taso jäljessä?").

Hahnottamista voisi testata piirrätyttämällä radix-2 DIT FFT jollain isolla N:n arvolla. pari vuotta sitten välikokeissa tällaisia väittämiä olikin. Keskiviikkona pohdittiin hetki, miltä näyttäisi N=16 diagrammi. Tai miten tämä ohjelmoitaisiin - kenties rekurssiivisesti.

Liitutaulukuvissa lisänä kommentteihin on laskennan aikana tulevat tulokset, perhosyhtälöiden aukilasku (tiistai) sekä Psi[2]=X[2] osoittaminen samaksi kuin DFT:n määritelmän mukaan.

Ensin siis sama x[n] sekvenssi DFT:n määritelmän mukaan:



ja sitten radix-2 DIT FFT:



...omaa laskentaa...




FFT on kaunista.

Tehtävä 2 käytiin tiistaina nopeasti (turhaa?) ja keskiviikkona melkein täysin sivuutettiin. Tulee luennolla ensi maanantaina. Lisämatskussa siitä on esitetty ratkaisuna vain taajuuspuolen alasnäytteistyksen kaavan soveltaminen - siis matemaattista silmää vaativa tapa:

R11 - Paperi #6 -- kommentteja

Kurssin loppupuoliskolla käsittellään ilmiöitä, joita syntyy kun suotimen siirretään "matemaattisen täydellisestä, symbolisesta" maailmasta kohti "epätäydellistä, epätarkkaa".

Luku pi on irrationaaliluku, voimme kirjoittaa sen yhdellä symbolilla, joka on matemaattisen tarkka. Mutta jos haluamme kirjoittaa sen murtoluvuksi tai tallettaa tietokoneen muistiin, joudumme käyttämään pyöristystä. Pyöristykset tulevat konkreettisiksi erityisesti käytettäessä pientä bittimäärää.

Tämän viikon laskareissa pääpaino on diskreetin Fourier-muunnosta nopeasti laskevassa algoritmissa, esimerkkinä radix-2 DIT FFT. Tätä kysytään myös välikokeessa.

Toinen esimerkki on signaalin alas- ja ylösnäytteistys digitaalisesti ilman D/A, A/D-muuntimia.

Tehtävä 1, pohjustus.

Katsotaan ensin määritelmän mukaista diskreettiä Fourier-muunnosta ja erityisesti siinä esiintyvää termiä W_N. Nyt jos N=4, niin W_N^kn-pisteet ovat tasavälein yksikköympyrällä.



Nyt voidaan siis vaikkapa käsin laskea DFT jonolle x[n] = {_2_, 3, 5, -1}.
Jono x[n] on neljä merkkiä pitkä, N = 4. Tällöin myös muunnos X[k] on neljä merkkiä pitkä: X[0], X[1], X[2], X[3]. Kukin X[k] lasketaan summana neljästä termistä x[n]W_N^kn. Alla laskettuna X[0] ja X[1]. Laske X[2] ja X[3].



Jos muunnettava jono x[n] on reaaliarvoinen (niinkuin yleensä onkin), niin muunnoksella X[k] on tiettyjä symmetrisiä ominaisuuksia:



DFT:ä tarvitaan siis vaikkapa signaalin taajuusestimoinnissa. Jatkuvasta signaalista voidaan ikkunafunktiolla leikata pala, esim. 256 merkkiä pitkä. Hamming-ikkuna (kuten Hann, Blackman, ...) vaimentavat pätkän alun ja lopun nollaan, jotta spektriin ei tule "säröä" epäjatkuvuudesta johtuen.



Kun tällä lasketaan DFT-256, niin taajuusakseli 0..f_t jaetaan 256 palaan. Taajuusresoluutioksi f_T=8000 Hz tulee noin 40 Hz. Olennaista on siis huomata, että jälleen näemme signaalista vain kaistan 0..f_T/2 eli ensimmäiset 129 DFT:n arvoa.



Kaikki tämä siis voitaisiin laskea kynällä. Lisäksi voidaan osoittaa, että joka kerta lasketaan jo aiemmin laskettuja x[n]W_N-termejä. FFT-algoritmit laskevat DFT:tä yhtä tarkasti mutta "viisaammin", jolloin laskuoperaatioiden määrä putoaa huomattavasti. Tätä esitetään O-notaatiolla. Määritelmän mukaisessa DFT:ssä laskuoperaatoita on O(N^2) kun taas FFT:ssä yleisesti O(N log_2 N).



Esitellään radix-2 DIT FFT, jonka algoritmi voidaan visualioida "virtauskaavioksi". DFT-muunnettava jono x[n] asetaan tietyssä järjestyksessä vasempaan reunaan, lasketaan kerros kerrokselta oikealle päin kunnes oikeassa reunassa on DFT:n arvot X[k].



Tässä algoritmissa x[n] järjestetään bittikäänteiseen järjestykseen. Esim. kymmenjärjestelmän luku 11 on 1011, niin bittikäännetty järjestys on 1101, jota vastaa kymmenkannan 13. Kun N=4 ja bittejä siis b=2, koska 2^b = 2^2 = 4, niin {0, 1, 2, 3} ~ {00, 01, 10, 11} --> {00, 10, 01, 11} ~ {0, 2, 1, 3}. Tämä on "desimointia". Lasketaan kompleksiset kertoimet W_N^s. Tämän jälkeen taso tasolta eteenpäin. Visuaalinen ilme on peräisin "perhosyhtälöistä".



Esimerkkimateriaalissa on annettu esimerkki tilanteesta N=8. Vertaile N=4 ja N=8. Miltä näyttäisi N=16. Tai N=256?



Tehtävä 2, pohjustus

Monen näytteenottotaajuuden (multirate) tekniikoita ja järjestelmiä voidaan tarvita vaikkapa signaalin näytteenottotaajuuden muuttamisessa käyttämättä D/A-A/D-muunnosparia, joka luo aina virhettä tai kun halutaan analysoida tai käsitellä signaalin eri taajuuskaistoja (suodinpankit, filter bank), esim. audiokoodauksessa kriittiset kaistat. Uudelleennäytteistys voi toki tehdä vaikkapa splineillä tai Lagrangen kertoimilla - digitaalinen sämpläys on DSP-juttu. Multirate-järjestelmillä voidaan toteuttaa suotimia tehokkaammin kuin "suoralla toteutuksella".

Alasnäytteistyksessä vähennettään näytteenottotaajuutta jollain kokonaisluvulla. Aikatasossa otetaan vain joka M:s näyte. Taajuuspuolelta nähdään, että uhkana on signaalin vierastuminen.



Ylösnäytteistyksessä näytteenottottaajuutta kasvatetaan kokonaisluvulla L. Aikapuolella jokaisen alkuperäisen näytteen väliin sijoitetaan L-1 kpl nollia. Tämä tuottaa signaaliin korkeita, ylimääräisiä taajuuksia ("image").



Uudelleennäytteistyksessä voidaan luoda L/M-suhde, jonka lisäksi tarvitaan sopiva alipäästösuodatus toisaalta poistamaan "image" ja vierastumisen vaara.



Tehtävä 2 on pieni esimerkki alas- ja ylösnäytteistyksen vaikutuksesta taajuustasossa. Muista, että kuten analogisen signaalin näytteistyksessä kuin myös tässäkin, mitään korkeita taajuuksia ei katoa vaan ne mahdollisesti laskostuvat matalille taajuuksille.

Alla olevassa ratkaisumallissa harjaannutetaan matemaattista silmää käyttämällä alasnäytteistyksen matemaattista kaavaa.



Ylösnäytteistyksessä samoin, mutta tässä tapauksessa hieman helpompi. Huomaa, että kun ensin pudotetiin 1/3-osaan ja sitten nostettiin 3-kertaiseksi, niin lopullinen näytteenottotaajuus on sama kuin alussa. Kaistan pi/3..2pi/3 lisäksi on tullut komponentit 0..pi/3 ja 2pi/3..pi -kaistoille. Nämä voidaan nyt erottaa LP-, BP- ja HP-suotimilla H_0, H_1 ja H_2.

R10 Matlab #5 - kommentteja

Tällä kierroksella kaksi demoa: puhesignaalin piirteistä ja kuvankäsittelystä. Molempien Matlab-koodi valmiina kurssin www-laskarisivulla.

Tehtävä 1. Kannattaa palauttaa mieliin ma 6.4. Mikko Kurimon luento automaattisesta puheentunnistuksesta. Tässä tehtävässä (a) lasketaan pienistä aika-ikkunoista signaalin energiaa, (b) lasketaan ikkunoista spektrit ja piirretään ne spektrogrammina, (c) ryhmitellään aikaikkunoita vastaavat spektrit kolmeen ryhmään automaattisesti k-means-algoritmilla. Matlab-koodin lopussa eri ryhmiä voi kuunnella. HUOM! koodin lopussa on pause-komento, joka odottaa käyttäjältä enter-lyöntiä.

Breakpoint 1: size(x), size(E), E.

Breakpoint 2: size(S): 17 taajuuskaistan energia-arvo eli 17-ulotteinen vektori, näitä 28 kpl.

Breakpoint 3: .. Kokeile eri arvoja k eli ryhmien lkm. /s/ näyttää erottuvan aina.

Tehtävä 2. Ideaalinen alipäästösuodatus valokuvalle.

Tässä luetaan cameraman.tif, joka on 256 pikseliä korkea ja 256 pikseliä leveä harmaasävykuva. class(I) antaa luokaksi uint8 eli "unsigned integer 8" eli etumerkitön kokonaisluku 0 .. (2^8 -1) eli 0 .. 255. Nolla vastaa mustaa ja 1 valkoista. Voi syöttää riville I ja lyödä enter, niin näkee kaikki numerolukuarvot. Katsotaan alla minimi- ja maksimiarvot sekä "leikataan" pieni pala uuteen kuvaajaan!

Kaikki muunnokset, konvoluutiot, spektrit yms. voidaan ajatella 2D-kuvalle yhtälailla kuin 1D-äänisignaalille. Alla olevissa kahdessa kuvassa (taulu ja Matlab-koodin tuottamat kuvaajat 1..5):

Ylärivissä "aikataso" ("spatiaalitaso"), vasemmalla alkuperäinen, keskeltä puuttuu 2D-suotimen impulssivaste (maski) ja oikealla konvoluution avulla saatu ulostulo.

Alarivissä alkuperäisen 2D-signaalin 2D-spektri, suotimen spektri ja oikealla näiden tulo.

Ylä- ja alarivi välillä liikutaan kaksiulotteisen diskreetin Fourier-muunnoksen fft2 avulla. Huomaa, että kun nyt taajuustasossa ideaalinen alipäästösuodin,
niin aikatasossa siitä tulisi jonkun sortin 2D-sinc.

Ideaalisen suodatuksen lopputulos: alipäästösuodatus pehmentää nopeita vaihteluita (reunoja). Kuvankäsittelysoftissa komentoja "blur", "smoothen" ,"unsharpen", "average". Lisäksi syntyy ihmisen silmään ikävän näköistä reunojen aaltoilua. Kun palauttaa mieliin, että konvoluutiossa suotimen 2D-impulssivaste on "sinc", niin sieltähän ne aaltoiluthan ovat peräisin.

Voit kokeilla pehmentää taajuuspuolen suodinta. Tai vaihtoehtoisesti ajattele, että taajuuspuolella LP-suodin olisi sinc-tyyppinen, jolloin aikatason suodin on puhdas NxN-pisteen keskiarvoistava suodin (muista äänisignaalin esimerkki MA-2).

Luento 6.4. Puheentunnistus

Mikko Kurimo puhui automaattisesta puheentunnistuksesta. Paikalla noin 60 opiskelijaa. Materiaali kurssin www-sivulla (ei Nopassa).

Luennolla ei esitetty kysymyksiä. Jos/kun sinulle jäi jotain, jota
olisit voinut kysyä, mutta et ehtinyt, tee se nyt jälkikäteen sähköpostilla
tai kurssin nyyssiryhmässä. Tässäkin blogissa on kommentointimahdollisuus.

R9 / Paperi #5 kommentteja

Tällä kertaa 1,5 tuntia vierähti suodinsuunnittelun merkeissä. Otetaan alkuun perustilanne, jossa x -- h -- y. Otetaan vaikkapa äänisignaali x ja piirretään siitä spektri. Puhe- tai musiikkiäänessä yleensä suurin osa energiasta on matalilla taajuuksilla. Jos havaitaan korkeataajuista häiriötä, halutaan suunnitella LTI-suodin, joka on alipäästösuodin sellaisella rajataajuudella, että häiriö vaimentuu pois (ylempi kuva). Jos valitsemme rajataajuuden väärin (vrt demotehtävän kuva 1c), niin häiriö ei suodatukaan oikein (alempi kuva).

Tehtävä 1. Suotimen vaatimusmäärittely (speksit) kirjoitetaan rajataajuuksien ja värähtely/vaimenemisominaisuuksien suhteen. Tehtäväannosta piirretään seuraavat visuaaliset määrittelyt.

Mallivastauksissa speksit ja toteutunut |H(z)| on piiretty kuvaajiin. Jos |H(z)| käppyrä menee speksien väliin, niin ok. Jos ei, niin jotain on tehtävä, jotta speksit täyttyy: muokattava rajataajuutta, värähtelyspeksejä tai lisättävä astelukua, mikä vaikuttaa suotimen jyrkkyyteen ja värähtelymääriin.

Havaitaan, että (a) 4. asteen elliptinen täyttää annetut speksit, (b) 10. asteen Cheb II täyttää, mutta on liiankin tiukka etenkin estokaistalla, (c) rajataajuuden suhteen ei täytä speksejä.

Liian tiukka suodin == turhan suuri asteluku vaatimuksiin suhteutettuna, turhaa laskentaa, turhia muuttujia, turhaa prosessorin ja ilmakehän lämpenemistä jne...

VK2-pistelaskarit tehtävä 2: kerroin K skaalaa suotimen magnitudivasteen maksimin ykköseksi. Skaalauksilla ei merkitystä symbolisessa, ideaalissa laskennassa, mutta on merkitystä, jos esitystarkkuus on rajattu.

Tiedetään siis, että 8. asteen kaistanpäästösuotimen maksimi on kohdassa omega = pi/2. Koska z = e^jw, niin z= e^(j pi/2) = j, ja maksimiarvo voidaan "helposti" laskea.

Tässä tuli sattumoisin esimerkki äärellisestä laskentatarkkuudesta. Tehtäväpaperissa oli neljä desimaalia, laskin tiistaina taululla katkaisulla yhteen ja keskiviikkona pyöristyksellä yhteen desimaaliin. Katkaisulla pääsin tilanteeseen, jossa jakaja oli nolla, pyöristyksellä maksimiarvoksi tuli 7, kun se oikeasti lienee 17,.. tai jotain. Tässä nyt erityisesti vaikutti se, että nimittäjä on pieni luku, joten laskenta on herkkää.

K (gain) voi ajatella olevan "voimakkuuden säädin". Suotimien sarjaankytkennässä sisäinen laskenta skaalataan sopivaksi, jotta laskentatarkkuus olisi järkevää.

Voi ajatella myös signaalin näytteistyksessä, esim. puheen taltioinnissa. Jos puhutaan mikkiin hiljaa tai kaukaa, niin signaali saa koko dynaamisen skaalan -1..+1 sijasta arvoja vaikkapa välillä -0.01..+0.01. Nyt jos lyödään K=100, niin päästään takaisin tuohon skaalaan -1..+1, mutta myös kvantisointikohina on yhtä lailla vahvistunut!

Kotistereoissa nupit kaakossa lienee suurin piirtein 0 dB ja sitten hiljaisimmillaan kuiskaustasolla noin -70 dB. Mites se menikään: puhe 60 dB, suihkukone 130 dB tai jotain... - sama suhde, lieneekö oikealla hehtaarilla?

Tehtävä 2. Digitaalisen IIR-suotimen suunnittelu bilineaarimenetelmällä. Teoreettinen tehtävä: Matlabissa tehtäisiin kahdella komennolla IIR-suodin. Voisi ajatella neljä vaihetta:
(I) speksit kohdasta 2b
(II) prewarppaus kohdasta 2c (s-taso -> z-taso)
(III) analoginen suodin H(s) kohdasta 2a
(IV) bilineaarimuunnos 2c (z-taso -> s-taso)

(I) digispeksit, nyt tiedetään vain, että rajataajuudella 100 Hz vaimennus -3 dB

(II) taajuuden esikorjaus tarvitaan kompensoimaan bilineaarimuunnoksen (kohta IV) taajuusvääristymä. Digimaailmassa nähdään taajuudet 0 .. fT/2, tässä 0 .. 500 Hz. Jatkuvassa maailmassa taasen nähdään taajuudet 0 .. oo Hz. Eli vääristymä eestaas tulee, kun kohta 500 Hz mäpätään äärettömäksi taajuudeksi. Kohta Omega_pc s-tasossa vastaa taajuutta 104 Hz (ravistettu hihasta tulkkaamalla vakio k = 2f_T).

Vakiosta k ei tarvitse tietää mitään, koska se tulee automaagisesti katoamaan kohdan IV sijoituksessa.

(III) Nyt siis meillä on _Butterworth_-tyyppinen 1. asteen stabiili H(s)-suodin. Siitä tullaan saamaan kohdassa (IV) 1. asteen stabiili H(z)-suodin. Muistetaan, että stabiililla analogiasuotimilla navat vasemmassa puolitasossa, kun taas stabiiliilla digisuotimella navat yks.ympyrän sisällä.

(IV) Bilineaarimuunnos sijoituksilla s = ... ja Omega_pc = ...

Tehtävä 3. Digitaalisen FIR-suotimen suunnittelu. Erilaisia vaihtoehtoja. Matlabin komentoja fir1, fir2, firls, firpm.

Katsotaan tässä ikkunamenetelmä (fir1).

Ikkunafunktioita käytetään erityisesti spektriestimoinnissa katkaisemaan signaali x pieniin paloihin. Siksi termit Hamming, Hann, Blackman voivat ponnahtaa esiin.

Pääperiaatetta voi demota myös kurssin www-sivuilta löytyvällä Matlab-ohjelmalla demoFIRwindow.m

Lasketaan siis ideaalinen, äärettömän pitkä h_d[n], ikkunoidaan se w[n]:llä:

h_FIR[n] = h_d[n] . w[n]

Kuten ylläolevasta kuvasta kävi ilmi, kaikkia vaiheita voi tutkia niin aika- kuin taajuustasossa.

R8 Matlab #4 laskarien kommentteja

Kurssin alkupuolella suotimen analysointia. Nyt vastakkainen operaatio eli määrittelyistä lasketaan suotimen H(z) kertoimet.

HUOM! Tehtävien 2 ja 3 koodipohja on saatavilla kurssin kotisivulta, kuten myös muutamat muut apufunktiot speksit.m, speksitFIR.m, tf2latex.m.

[10-15 min] Tehtävässä 1 käydään läpi suodinvaatimusten määrittelyä Matlabia varten. Tässä erityisesti IIR-suotimien määrittelyt: päästökaistalla maksimivaihtelu Rp (dB), estokaistalla minimivaimennus Rs (dB), päästökaistan rajataajuus Wp (skaalattuna, näytteenottotaajuuden puolikas vastaa ykköstä, joten aina 0 < Wp < 1) ja estokaistan rajataajuus Ws.

Syötetään arvot Matlabiin. speksit.m on kurssin oma funktio.

[15-20 min] Tehtävässä 2 käydään esimerkki IIR-suotimen suunnittelusta. Digitaalisen IIR-suotimen suunnittelussa käytetään hyväksi analogisen IIR-suotimen H(s) suunnittelualgoritmeja (Butterworth, ...), jotka bilineaarimuunnoksella muutetaan digitaaliksi H(z).

Ladataan synteettinen signaalin M4002.wav, joka koostuu 100 Hz, 300 Hz, 500 Hz, 700 H, ..., kosinikomponenteista.

Tehtävänannossa sivulla 2 ("Task") annetaan suotimen vaatimukset: Rp, Rs, Wp, Ws, ja lisäksi, että halutaan elliptinen IIR-suodin (ellipord, ellip).

IIR-suotimien toteutus Matlabilla on kaksi riviä koodia: suotimen asteluvun arviointi (-ord.m) ja itse suotimen kertoimien laskeminen.

LUE esim. "help ellipord" JA "help ellip", jotta tiedät, mitä arvoja annetaan funktiolle ja mitä se palauttaa!!!

Täältä pitäisi tulla siis nätti 3. asteen alipäästösuodin. Piirretään suotimen magnitudivaste samaan kuvaajaan vaatimusten kera ja todetaan, että vaatimukset toteutuvat ("meet specfications"). Suodatuksen jälkeen voi katsoa tulosta signaalin spektrogrammista.

[15-20 min] Tehtävässä 3 tehdään FIR-suodin Parks-McClellan-algoritmilla (firpmord, firpm).

Jos tulee huomanneeksi, että speksit eivät täyty tarkalleen, lue "help ellipord" :n lopusta "CAUTION 1":

Saadaan sitten lineaarivaiheisen suotimen nätit kuvaajat asteluvulla 40 (alkuperäinen) tai 42 (korjattu).

Kokeile muuttaa päästö- ja estokaistan väliä pienemmäksi. Miten asteluvulle käy?