i
i
“Legisa-vesti” — 2017/6/30 — 9:01 — page 68 — #1 i
i
i
i
i
i
ZANIMIVOSTI
STRNJENO VZORČENJE
Uvod
Za razumevanje nove metode strnjenega vzorčenja (snemanja, zajemanja),
ki ji v angleščini pravijo compressed sensing, compressive sampling, si
najprej oglejmo nekaj dejstev o digitalizaciji slik, zvoka . . . Pri tem bomo
spoznali matematična orodja, ki nastopajo tudi v novi metodi zajemanja
signalov.
Stiskanje
Digitalne fotografije v surovi obliki vsebujejo informacijo o intenziteti treh
osnovnih barv za vsak piksel posebej. To pomeni množico podatkov. Zato
pred posredovanjem fotografij že kamera ali pa naknadno uporabnik foto-
grafijo stisne, navadno v zelo uspešni in razširjeni JPEG format, ki tudi
nima problemov s patenti. Večinoma lahko stisnemo za faktor 4 praktično
brez izgube podrobnosti. Tudi stisk za faktor 10 daje navadno zelo do-
bre rezultate. Na slikah imamo namreč pogosto večje površine (kot modro
nebo), kjer ni podrobnosti. Tako lahko na sliki najdemo pravokotnike, ki
so približno enakomerno obarvani. Tak del slike lahko podamo z mnogo
manj podatki. Izvedba tega stiskanja je nekoliko zapletena in matematično
zanimiva.
Zaradi enostavnosti privzemimo, da imamo črno-belo sliko. Denimo, da
kamera daje 8 bitov informacije o osvetljenosti piksla. Torej imamo 28 = 256
odtenkov sivine, od povsem črne (0) do povsem bele (255). Razdelimo sliko
na pravokotnike velikosti 8× 8 pikslov. Tak del slike podamo najprej z od-
stopanji svetlosti pikslov od vrednosti 27 = 128. Dobimo matriko velikosti
8×8, ki jo imamo lahko tudi za vektor v R64. Na matriki opravimo dvorazse-
žno diskretno kosinusno transformacijo, DCT. Gre za razvoj po ortogonalni
bazi. V baznih matrikah nastopajo kosinusi, zelo podobno kot pri klasičnem
razvoju v Fourierovo vrsto. Enorazsežna DCT deluje na prostoru RN tako,
da vektorju x = (x(0),x(1), . . . ,x(N−1)) priredi vektor X v istem prostoru
po formuli
X(k) =
N−1∑
j=0
x(j) cos
( π
N
(
(j +
1
2
)
k
)
.
68 Obzornik mat. fiz. 64 (2017) 2
i
i
“Legisa-vesti” — 2017/6/30 — 9:01 — page 69 — #2 i
i
i
i
i
i
Strnjeno vzorčenje
Formula za dvorazsežno transformacijo je seveda obsežneǰsa in si jo lahko
ogledate na Wikipediji.
DCT je v tesnem sorodu z diskretno Fourierovo transformacijo (DFT),
le da so koeficienti pri razvoju po DCT realni.
DFT vsakemu vektorju x ∈ CN priredi vektor x̂ v istem prostoru, tako
da je
x̂(k) =
N−1∑
j=0
x(j) exp
(−2πkj
N
)
.
Za izvedbo DFT pri velikih N imamo odličen algoritem, imenovan hitra
Fourierova transformacija – FFT, ki jo mnogi dobro poznate.
Vrnimo se k naši 8 × 8 matriki osvetljenosti pikslov. Dvorazsežna dis-
kretna kosinusna transformacija da novo 8 × 8 DCT matriko iz DCT ko-
eficientov. Koeficient z indeksom (0, 0) je povprečje vrednosti elementov
originalne matrike. Na tej stopnji lahko iz DCT matrike rekonstruiramo
originalno matriko.
Ampak koeficienti z majhnima indeksoma so najpomembneǰsi, ker po-
dajajo dele slike, ki najbolj padejo v oči. Na tipični sliki so elementi DCT
matrike z vǐsjimi indeksi majhni v primerjavi z elementi zgoraj levo, saj ni
veliko visokofrekvenčnih elementov na sliki. Težko recimo pričakujemo, da
bomo na sliki imeli močno komponento z indeksom (6, 5) v obliki nekakšne
modificirane šahovnice s 7×6 polji, kot si lahko ogledate v ilustracijah baze
za DCT na [1, str. 15]. Na 8× 8 pikslov velikem delu modrega neba bo po
absolutni vrednosti daleč največji element z indeksom (0, 0).
Zdaj naredimo kvantizacijo. To je enostaven in grob postopek, ki odseka
najmanj pomembne informacije in večinoma zmanǰsa natančnost preostalih
podatkov. Obstajajo posebne kvantizacijske matrike, dobljene izkustveno.
Elementi kvantizacijske matrike so naravna števila in večinoma naraščajo,
ko se gibljemo desno in navzdol po matriki.
V literaturi najdemo standardno kvantizacijsko matriko Q. Matrika Q
ima vrednost 16 na mestu (0, 0), 121 na mestu (6, 5) in 99 na mestu (7, 7)
– glej [2, str. 4] ali [1, str. 20].
Elemente DCT matrike delimo z istoležnimi elementi kvantizacijske ma-
trike, odvržemo decimalke. Prav tako vse vrednosti nad 255 porežemo na
255. Podobno naredimo z vrednostmi pod −255. (Angleški izraz za rezanje
je clip.) Gre za agresivno zmanǰsanje in pogosto eliminacijo že tako majhnih
visokofrekvenčnih komponent ter zmanǰsanje velikosti elementov kvantizi-
rane matrike. Elementi v kvantizirani matriki so do predznaka predstavljeni
z največ 8 biti. To je bistven korak na poti do stisnjene JPEG slike. Za-
Obzornik mat. fiz. 64 (2017) 2 69
i
i
“Legisa-vesti” — 2017/6/30 — 9:01 — page 70 — #3 i
i
i
i
i
i
Zanimivosti
vrgli smo del informacije na originalni sliki (večinoma nezanimivi del), se
zadovoljili s približki nekaterih drugih podatkov in originala ne moremo več
natančno rekonstruirati. Na tipični sliki ima kvantizirana matrika mnogo
ničel, predvsem v desnem spodnjem delu. Zgoraj omenjena matrika Q obi-
čajno sliko stisne za faktor približno 7. Matriki, v kateri je večina elementov
ničelnih, preostali pa nimajo posebne strukture, rečemo razpršena matrika.
Kvantizirana matrika je torej praviloma razpršena.
V fotoaparatu z velikim senzorjem (APS-C ipd.) nastavitev na fino (an-
gleško fine) da kvantizacijsko matriko z bistveno manǰsimi elementi, velikosti
recimo od 1 do 6. To pomeni nižjo kompresijo, nekako za faktor 2. Pri malih
tipalih z diagonalo pod 8 mm bo kvantizacijska matrika v načinu fine imela
elemente recimo od 1 do 15, saj ustrezne optike običajno nimajo zelo dobre
ločljivosti.
Pri dekodiranju pomnožimo matriko nazaj z istoležnimi elementi kvan-
tizacijske matrike in opravimo inverzno transformacijo k DCT. Dobimo pri-
bližek prvotne slike našega kvadrata. Na slikah z mehkimi prehodi med
svetlimi in temnimi deli slike to deluje sijajno. Algoritem za JPEG stiska-
nje je računsko nezahteven, hiter in robusten. Manǰsi problem se pojavi pri
slikah z ogromno podrobnostmi, kot so trava, krzno. Bolǰse kamere tako
sliko prepoznajo in bistveno manj stisnejo, se pravi uporabijo drugo kvan-
tizacijsko matriko kot sicer. (Slika z ogromno šuma je tudi težava, vendar
pa tu nismo zainteresirani za podrobno reprodukcijo.) Pri poceni kamerah
pa lahko kombinacija nekakovostnega zoom objektiva in neprilagodljivega
stiskanja travnik spremeni v zeleno plundro. JPEG tudi ni najbolǰsi za re-
produkcijo grafičnih podrobnosti. Za manǰse risbe in grafike profesionalci
raje uporabljajo format PNG.
Za zvok je nastal na podlagi JPEG priljubljeni, za zdaj še patentirani
format MP3. Omogoča stiskanje v različnih kakovostih. Zelo dober za
kompresijo zvoka je tudi prostokodni format (Ogg) Vorbis.
Vzorčenje in digitalizacija
Nekateri študenti na izpitih rǐsejo grafe funkcij »po točkah«. Večinoma se to
ne obnese. Vendar pa je mogoče velik razred funkcij popolnoma rekonstru-
irati iz njihovih vrednosti na diskretni množici točk. V knjigi [3] najdemo
na str. 373 izrek, ki ga ni težko dokazati:
Izrek 1. Naj bo f ∈ L2(R) zvezna in naj bo njena Fourierova transformi-
ranka f̂ enaka 0 zunaj intervala [−L,L], kjer je L > 0. Potem je f določena
70 Obzornik mat. fiz. 64 (2017) 2
i
i
“Legisa-vesti” — 2017/6/30 — 9:01 — page 71 — #4 i
i
i
i
i
i
Strnjeno vzorčenje
z vrednostmi (= vzorci) v točkah{ n
2L
;n ∈ Z
}
takole:
f(t) =
∞∑
n=−∞
f
( n
2L
)sin(2Lπt− nπ)
2Lπt− nπ
. (1)
Našo funkcijo f torej lahko rekonstruiramo iz njenih vrednosti = vzor-
cev v točkah, medsebojno oddaljenih za intervalček ∆t = 1/(2L), ki mu
pravimo Nyquistov interval. V eni sekundi torej vzorčimo 2L-krat, in to
je Nyquistova frekvenca.
Knjiga pravi, da izrek iz leta 1915 pripada E. T. Whittakerju.
Kakšna je Fourierova transformiranka dušenega kosinusnega vala f(t) =
exp(−bt2) cos(2πνt) s krožno frekvenco ω = 2πν? Če je b majhen (tako da
je val opazen dalj časa), je lahko izračunati, da ima f̂ dva ozka vrhova v
±ν in je bolj ali manj koncentrirana okrog teh dveh točk. Čim manǰsi je b,
tem bolj izraziti sta ti konici. Če je ν nekoliko manǰsi od L, je f̂ praktično
enaka nič zunaj intervala [−L,L]. Ker je f ∈ L2(R), približno zadošča
pogojem izreka. Tako laže razumemo znameniti Shannonov izrek (1948), ki
poenostavljeno pravi: za vzorčenje signala, ki vsebuje le frekvence pod ν,
moramo vzorčiti vsaj s frekvenco 2ν. Podobne rezultate so že prej dokazali
drugi matematiki: Kotelnikov leta 1933 . . .
Ljudje navadno ne slǐsimo zvokov nad 20 kHz. Torej moramo za dobro
digitalizacijo glasbe najprej odfiltrirati vse zvoke s frekvenco nad 20 kHz
in nato vzorčiti s frekvenco več kot 40 kHz. Za CD kakovost vzorčimo
s frekvenco 44,1 kHz. To nam omogoča, da iz vzorčenega signala verno
rekonstruiramo (sfiltrirani) original. Pri digitalni telefoniji vzorčimo le z 8
kHz, saj so telefoni konstruirani za frekvenčno območje od 300 Hz do 3,4
kHz. To je bolj zasilna rešitev, saj človeški glas sega nekako od 60 Hz do
14 kHz. Gosteǰse vzorčenje in bolǰso kakovost zvoka že zdaj ponuja Skype.
Uvajajo se tudi nove rešitve kot HD Voice (vzorčenje s 16 kHz).
Aliasing
Če vzorčimo s premajhno frekvenco, lahko pride do problemov, ki jih opi-
šemo z angleško besedo aliasing. Pri tem se lahko sinusoida v originalu
reproducira po digitalizaciji kot sinusoida s povsem drugo frekvenco. Ori-
ginalna sinusoida dobi torej svoje drugo življenje – svoj alias, ki pa nam
Obzornik mat. fiz. 64 (2017) 2 71
i
i
“Legisa-vesti” — 2017/6/30 — 9:01 — page 72 — #5 i
i
i
i
i
i
Zanimivosti
Slika 1. Funkciji f(x) = sin(2πx) (prekinjena črta) in g(x) = − sin((2/3)πx) (polna črta)
se ujemata v točkah
{
3k
4
; k ∈ Z
}
.
dela le težave. Pomislite recimo na funkcijo x 7→ sin(2πx). Tu je ν = 1.
Denimo, da jo vzorčimo s korakom 0,5. Če začnemo v x = 0, so vsi vzorci
enaki 0. Rekonstrukcija signala iz teh vzorcev nam da napačen rezultat –
funkcijo, ki je identično enaka 0. Če našo sinusoido vzorčimo s korakom
3/4 in začnemo v 0, vzorci ustrezajo tudi funkciji x 7→ − sin((2/3)πx), kot
vidimo na sliki 1. Pri rekonstrukciji vzorčenega signala bomo v tem primeru
dobili alias prvotnega signala v obliki funkcije z eno tretjino frekvence origi-
nala in z drugačno fazo, se pravi spet povsem napačen rezultat. Za pravilno
rekonstrukcijo moramo vzorčiti s korakom, manǰsim od 0,5.
Kaj pa, če ne bi vzorčili v enakomernih presledkih, ampak v slučajno
izbranih točkah? V tem primeru tudi pri vzorčenju s korakom, dalǰsim od
0,5, skoraj gotovo ne bi prǐsli do gornjih napačnih rezultatov.
V digitalni fotografiji se problemom z aliasingom poskušajo izogniti z
anti-aliasing (AA) filtri pred senzorjem. Ti zameglijo (sfiltrirajo) preveč
drobne detajle, a obenem malce zmanǰsajo kakovost slike. Efekte aliasinga
v fotografiji navadno označimo s francosko besedo moiré, ki jo izgovorimo
kot muaré.
Strnjeno vzorčenje
Vzorčenje s slučajnim korakom je ena od idej za t. i. strnjeno vzorčenje, an-
gleško compressed sensing, compressed sampling. Začetek tega sega
72 Obzornik mat. fiz. 64 (2017) 2
i
i
“Legisa-vesti” — 2017/6/30 — 9:01 — page 73 — #6 i
i
i
i
i
i
Strnjeno vzorčenje
v sedemdeseta leta preǰsnjega stoletja, kot lepo opisuje poljudno napisani
vir [4]. Geologi, ki so iskali nafto, so eksperimentalno ugotovili, da lahko
dobijo želeno sliko zemeljskih slojev z mnogo manj seizmičnimi meritvami,
kot bi jih bilo treba po teoriji. To je bila seveda super novica, saj je pov-
zročanje mini potresov z eksplozivi moteče. Razložili so si zadevo s tem, da
so posamezni sloji precej homogeni. Zanimale so jih le meje med sloji, ki
so na instrumentih dajale špice, in pa prisotnost sloja, ki bi lahko vseboval
nafto. Skratka, zanimalo jih je le majhno število pomembnih informacij.
V devetdesetih letih so metode geologov začeli preučevati matematiki, med
njimi David Donoho z Univerze Stanford.
Magnetna resonanca (MR) je odlična, ampak draga metoda za slikanje
mehkih tkiv. Francoski matematik Emmanuel Candès s Caltecha in njegov
mentor pri doktoratu Donoho sta oba začela ugotavljati, da lahko dobita
dobre MR slike vzorcev tudi s precej manj posnetki, kot jih predvideva
teorija, če vzorčita nekoliko drugače. Kot pri običajnih fotografijah tudi pri
MR navadno ni vse na sliki zanimivo.
Kot smo videli zgoraj, lahko digitalno sliko ali zvok skoraj zmeraj močno
stisnemo z minimalno izgubo kakovosti. Ampak najprej zajamemo čisto vse
podatke, da se nato odločimo, kako jih bomo stisnili. Podobno dobri slikarji
lahko hitro prepoznajo bistvo portretiranca ali naravne scene in ga podajo
z nekaj potezami.
Pri strnjenem vzorčenju želimo že vnaprej zajeti manj podatkov, kot
jih zahteva Shannonov pogoj. Metoda deluje le takrat, ko ǐsčemo majhno
število informacij.
Spomnimo se na naslednji problem. Imamo 9 kovancev, med katerimi jih
ima 8 enako maso, eden pa je lažji. Le z dvema tehtanjema s staromodno
ravnovesno tehtnico, brez uteži, lahko izločimo lažji kovanec, če se tega
lotimo pravilno.
Candès pripoveduje, da je imel na konferencah ob predstavitvi svojih
rezultatov težave, ker so nekateri mislili, da malce prireja rezultate. Poskušal
je stvari strogo dokazati, pa se mu je po prvih uspehih zataknilo. Naključje
je hotelo, da je Candès imel otroka v isti mali šoli kot Terence Tao z UCLA.
(Tao je nekaj let kasneje, 2006, dobil Fieldsovo medaljo). Francoz je Tau
zaupal svoje probleme, a mu sogovornik sprva ni verjel in mu je obljubil
nasprotni primer za njegove trditve. Tega pa ni mogel najti, zato je stvar
vzel resneje, začel sodelovati s Candèsom in prispeval manjkajoči dokaz. Z
njima je delal Justin Romberg z Georgia Tech. Preprint njihovih rezultatov
je izšel leta 2004, obenem pa je Donoho oznanil podobna dognanja.
Obzornik mat. fiz. 64 (2017) 2 73
i
i
“Legisa-vesti” — 2017/6/30 — 9:01 — page 74 — #7 i
i
i
i
i
i
Zanimivosti
Na kratko bom povzel Taovo predavanje [5] iz leta 2009. Rekonstruirati
želimo signal x dolžine N , kjer je N ogromen, denimo reda velikosti 106. Za
x vemo poleg tega, da je s-razpršen, kar pomeni, da ima največ s koordinat
od 0 različnih. Tu je s mnogo manǰsi kot N , kar pǐsemo s  N . Imamo
poddoločen sistem m linearnih enačb za x v obliki Ax = b. Tu je b vektor
meritev in m < N . Vemo, da ima tak sistem poleg našega signala neskončno
rešitev.
Trditev 2. Privzemimo, da je vsakih 2s stolpcev matrike A linearno neod-
visnih. Če je signal x s-razpršen in zadošča sistemu Ax = b, je ta signal
edina s-razpršena rešitev tega sistema.
Dokaz. Denimo, da ima sistem še eno s-razpršeno rešitev x′. Vzemimo
vektor w = x−x′. Ta vektor ima največ 2s koordinat od 0 različnih. Velja
Aw = 0. Toda produkt na levi je enak w(1) krat prvi stolpec matrike
A+ . . ., kjer sumacija poteka po neničelnih koordinatah za w, torej ima ta
linearna kombinacija največ 2s sumandov. Iz predpostavke sledi, da so vsi
w(i) enaki 0 in tako w = 0.
Predpostavka trditve pomeni, da je število m vrstic v matriki A vsaj 2s,
se pravi, da je meritev vsaj 2s.
Posledica 3. Signal x je tista rešitev sistema Ax = b, ki ima najmanj
neničelnih koordinat.
To je zelo lep rezultat, žal pa je iskanje najbolj razpršene rešitve nume-
rično izredno zahtevno. V splošnem je to NP-težak problem. Izkazalo pa
se je, da namesto tega lahko ǐsčemo rešitev, za katero bo vsota absolutnih
vrednosti koordinat (se pravi `1–norma) vektorja x minimalna. Tej me-
todi iskanja rešitve so nadeli angleško ime basis pursuit. Gre za konveksno
optimizacijo, ki se je lahko lotimo z metodami linearnega programiranja.
Rezultati teorije so tesno povezani z verjetnostjo. Naslednji izrek so Can-
dès, Romberg in Tao objavili v IEEE Transactions on Information Theory
leta 2006:
Izrek 4. Izberimo naključno števila k1, k2, . . . , km iz množice {0, 1, . . . ,
N − 1}. Obstaja taka konstanta C, da lahko pri pogoju m > Cs lnN
z verjetnostjo blizu 1 vsak s-razpršen vektor x ∈ CN rekonstruiramo
iz vrednosti diskretnih Fourierovih koeficientov x̂(ki), kjer i teče od 1 do m.
74 Obzornik mat. fiz. 64 (2017) 2
i
i
“Legisa-vesti” — 2017/6/30 — 9:01 — page 75 — #8 i
i
i
i
i
i
Strnjeno vzorčenje
Normalno bi za rekonstrukcijo potrebovali vseh N Fourierovih koeficien-
tov. Numerični poskusi sugerirajo, da lahko večinoma rešitev dobimo točno
celo pri pogoju, da je meritev vsaj 4s.
Metoda deluje pri precej širokem razredu matrik, ne samo pri matriki
za diskretno Fourierovo transformacijo. Tu močno pomaga dejstvo, da je N
ogromen.
Metodo lahko razširimo tudi na primer, da je signal samo približno raz-
pršen; v tem primeru dobimo približno rešitev. Če je meritve pokvaril šum,
se pravi, da imamo enačbo Ax = b + z, kjer je z šum, bo naša metoda dala
približek rešitve.
Teorija je zelo zanimiva. Manǰsi problem je, da je metoda računsko
precej zahtevna.
Konec novembra 2016 je Siemens oznanil, da začenja prodajo opreme
za magnetno resonančno (MR) slikanje, ki sloni na strnjenem vzorčenju.
Namenjena je slikanju delovanja srca. Po reklami v 25 sekundah naredi
tisto, za kar so potrebovali prej do 6 minut. Pacient mora samo enkrat
zadrževati dih, namesto desetkrat do štirinajstkrat pri preǰsnjem postopku.
Obdelava podatkov poteka v sami napravi. Skrčitev morda ne gre samo na
račun strnjenega vzorčenja, ampak tudi drugih algoritemskih in inženirskih
izbolǰsav. Prve verzije te inovacije so v uporabi že vsaj leto dni v uglednih
bolnǐsnicah, tako da stvar deluje. Siemens obljublja, da bo kmalu ponudil
podobne naprave tudi v druge namene.
Strnjeno vzorčenje je lahko uporabno marsikje, kjer so senzorji zelo
dragi. Verjetno se bo strnjeno vzorčenje pojavilo tudi pri računalnǐski to-
mografiji (CT) in zmanǰsalo velike doze sevanja.
LITERATURA
[1] A. Bruna, JPEG, Catania, 2008, dostopno na www.dmi.unict.it/~battiato/EI_
MOBILE0708/JPEG\%20(Bruna).pdf, ogled: 6. 4. 2017.
[2] K. Cabeen in P. Gent, Image compression and the discrete cosine transform, Math
45, College of the Redwoods, dostopno na www.math.cuhk.edu.hk/~lmlui/dct.pdf,
ogled: 6. 4. 2017.
[3] G. Bachmann, L. Narici in E. Beckenstein, Fourier and wavelet analysis, Universitext,
Springer-Verlag, New York, 2000.
[4] D. Mackenzie, Compressed sensing makes every pixel count, What’s Happening in
Mathematical Sciences 7 (2009), 114–127, dostopno na www.ams.org/samplings/
math-history/hap7-pixel.pdf, ogled: 6. 4. 2017.
[5] T. Tao, Compressed sensing or: the equation Ax = b revisited, Mahler
lecture 2009, dostopno na www.math.hkbu.edu.hk/~ttang/UsefulCollections/
compressed-sensing1.pdf, ogled: 6. 4. 2017.
Peter Legǐsa
Obzornik mat. fiz. 64 (2017) 2 75