Chaotic Pearls

Lamunan dari seberang GKU lama …

Archive for the ‘Theory and Concept’ Category

Dekomposisi Sinyal dan Pencarian Basis Ideal (4)

Posted by suksmono on May 22, 2008

Masalah pencarian vektor basis ideal dari sebuah sinyal telah dapat diformulasikan sebagai proses optimisasi berkendala (constrained optimization) P1. Tulisan ini akan terlebih dahulu membahas secara sekilas konsep optimisasi linier, khususnya LP (Linear Programming), dan metode pemecahannya. Meskipun kendala dalam P1 sebenarnya bersifat tak-linier, pengaturan parameter yang tepat dan batas-batas yang sesuai dapat mengubahnya menjadi permasalahan pemrograman linier. Suatu contoh pencarian basis ideal dari sinyal hasil bentukan sepasang ortobasis dipakai untuk menunjukkan kegunaan LP. Jika batasan sparsity terpenuhi, LP dapat menemukan seluruh basis penyusun sinyal sehingga sinyal dapat direkonstruksi secara eksak.

Pemrograman Matematika dan Optimisasi

Didalam optimisasi ada suatu fungsi objektif yang akan dioptimumkan; dibuat minimum atau maksimum. Kendala memberi batas wilayah dimana fungsi ini berlaku. Dari berbagai macam metoda optimisasi yang ada, salah satu yang memiliki pemakaian luas adalah pemrograman linier (LP/Linier Programming), suatu bentuk khusus dari pemrograman matematika (Mathematical Programming). Bahkan metoda Simpleks yang ditemukan oleh G. Dantzig untuk mencari solusi masalah LP dinobatkan sebagai satu dari sepuluh algoritma terbaik abad keduapuluh.

LP pada dasarnnya merupakan masalah pencarian nilai optimum dari suatu fungsi linier dengan kendala yang juga berbentuk linier. Berikut ini contoh sederhana dari permasalahan LP:

Tentukan harga {x}_{1}  dan {x}_{2}  yang meminimumkan nilai {x}_{1}-{x}_{2}  , dengan syarat bahwa {x}_{1}+{x}_{2} \leq 10  , {x}_{1}-{x}_{2} \leq 5  , dan {x}_{1} \geq1 \: , {x}_{2} \geq 1  .

Masalah diatas dapat dituliskan kedalam bentuk baku sebagai berikut:

minimize {x}_{1}-{x}_{2}

subject to:

{x}_{1}+{x}_{2} \leq 10  ,

{x}_{1}-{x}_{2} \leq 5  ,

{x}_{1} \geq1 \: , {x}_{2} \geq 1  .

Nilai yang akan dioptimalkan disebut sebagai fungsi obyektif, dalam soal diatas adalah {x}_{1}-{x}_{2}  . Pemecahan masalah ini dapat dilakukan secara grafis dengan cara menggambarkan semua kendala, lalu mencari titik dimana nilai obyektif mencapai optimum dalam daerah solusi:

Linear Programming

Gambar 1. Ilustrasi pencarian titik optimum dalam LP

Pada Gambar 1 diatas, daerah yang diarsir adalah pasangan \left( {x}_{1},{x}_{2} \right)  yang memenuhi semua kendala yang diberikan atau disebut juga feasible region. Kita harus mencari satu titik dimana nilai fungsi obyektif mencapai optimal. Pada LP, solusi selalu terletak di pojok wilayah daerah atau titik simpul. Program Matlab berikut dapat dipakai untuk mencari solusinya:

%pemrograman linier untuk masalah sederhana
% min c’x s.t. Ax<=b; x1>=1, x2>=1
%——————————————
c=[1;-1]; % koefisien fungsi obyektif
A=[1 1; 1 -1]; %matriks A
b=[10; 5]; % kendala
lb=[1 1] ; %batas bawah x1 dan x2
[x, fval,exitflag,output,lambda]=linprog(c,A,b,[],[],lb);
disp(sprintf(‘x1= %f x2=%f\n’,x(1),x(2))); %tampilkan

Dari eksekusi program akan diperoleh nilai {x}_{1}=1, {x}_{2}=9  yang akan menghasilkan nilai fungsi objektif paling kecil, yaitu -8. Posisi titik optimum ini terletak dipojok paling kiri atas dari politop yang menyatakan daerah solusi.

Metoda simpleks dari Dantzig melakukan perunutan setiap tepian wilayah hingga ditemukan titik optimal yang selalu terletak di sudut atau verteks dari politop. Dengan demikian, dapat ditunjukkan bahwa pada kondisi terburuk, metoda ini akan memiliki kompleksitas yang tumbuh secara eksponensial. Namun kebanyakan masalah praktis dapat dipecahkan dengan cukup cepat, karena kompleksitas rata-ratanya cukup kecil hingga dapat ditangani.

Masalah pertumbuhan kompleksitas perhitungan dalam metoda Simpleks akhirnya dapat diatasi oleh Narendra Karmakar dengan metoda interior point (IPM), yang memiliki kompleksitas polinomial. IPM tidak hanya dipakai pada optimisasi linier, kasus lain yang taklinier-pun dapat dipecahkan dengan metoda ini asalkan daerah solusi bersifat konveks. IPM mencapai kompleksitas polynomial karena tidak lagi merunut tepian daerah solusi, tapi masuk menerobos kedalam wilayah ini sampai solusi yang dicari dapat ditemukan.

Masalah optimisasi seperti diatas dapat dituliskan dalam bentuk umum sebagai berikut:

minimize {c}_{1} {x}_{1} + {c}_{2} {x}_{2} + ... + {c}_{n} {x}_{n}

subject to:

{a}_{11} {x}_{1}+ {a}_{12} {x}_{2} + ... + {a}_{1n} {x}_{n}

{a}_{21} {x}_{1}+ {a}_{12} {x}_{2} + ... + {a}_{2n} {x}_{n}

.........................................

{a}_{m1} {x}_{1}+ {a}_{m2} {x}_{2} + ... + {a}_{mn} {x}_{n}

{x}_{1}, {x}_{2}, ..., {x}_{n}   \geq 0

bahkan dengan dengan menggunakan slack variables

{x}_{n+1} = {b}_{j} \sum_{j=1 .. n} {a}_{ij} {x}_{j}  , i = 1, 2, …, m

LP dapat dituliskan dalam bentuk matriks yang lebih ringkas, yaitu

minimize {c}^{T} x  subject to Ax=b, \: x \geq 0

Masalah mencari nilai minimum suatu fungsi obyektif dapat diubah menjadi pencarian maksimum fungsi obyektif yang lain dan demikian pula sebaliknya. Dengan demikian, pencarian nilai minimum diatas memiliki solusi yang sama dengan masalah dual-nya, yaitu pencarian nilai maksimum. Jika fungsi objektif dan kendala tidak lagi linier, fungsi dari peubah x ini dapat ditulis sebagai c(x) dan kendala dapat ditulis sebagai {a}_{i}(x) \geq {b}_{i}  .

Metode Basis Pursuit dengan Pemrograman Linier

Metoda basis pursuit merumuskan pencarian basis ideal {\phi}_{m}  didalam kamus-basis \Phi  berdimensi NxM sebagai masalah pemrograman konveks:

(P1): min \; {||\alpha ||}_{0} \; s.t. \; \Phi \alpha = s

Telah diketahui sejak lama bahwa optimisasi non-linier l1 seperti ini ekivalen dengan LP

min \; {c}^{T}x \; s.t. \; Ax=b, x \geq 0

jika dilakukan substitusi sebagai berikut

M↔2M; x↔( {\alpha}_{1},{\alpha}_{2})  ; c↔[1 1 …. 1; 1 1… 1]; A↔ \left[ \Phi, -\Phi \right]  ; b↔s

Artinya, kamus-basis digabung dengan nilai negatif dari kamus tersebut untuk menggantikan peran matriks A, batas kendala b digantikan dengan sinyal s, panjang vektor solusi dibuat menjadi dua-kalinya, sedangkan fungsi objektif adalah penjumlahan masing-masing koefisien vektor solusi sparse \alpha  dan dimensi kolom diperbesar menjadi duakali-nya.

Cara ini selanjutnya akan dipakai untuk memecahkan masalah dekomposisi sinyal s terdahulu yang terdiri dari dua impuls dan dua vektor basis DCT. Sebagai tambahan, sinyal ini tercampur dengan derau acak Gaussian. Berikut ini kode Matlab-nya.

%———————————————–
%Basis pursuit on overcomplete basis
%Using Matlab Optimization Toolbox
%To search for solution … modified
%———————————————–
%length of signal
clear;close all;clc;
N=64; %signal dimension
over_factor=2; %over complete factor
nAtoms=N*over_factor;
%create pair of orto-basis dictionary
D1=eye(N,N); D2=dct(eye(N,N)); D=[D1;D2]‘;
%create a signal with N_support’
N_support=4;
pN=[8 32 73 125]; %dpt diganti dng indeks acak jika perlu
x=zeros(N,1); X=zeros(nAtoms,1);%’transformed-domain signal
for k=1:N_support;
x=x+D(:,pN(k)); X(pN(k))=1;
end;
x0=x;
figure(1);plot(x0);title(‘Original Signal’);
%perform measurement, >> add noise
derau=randn(N,1);
y = x0 + 5*derau/(derau’*derau);
min_y=min(y);
D_pos=-D;
%define dictionary A
A=zeros(N,2*over_factor*N);
A(:,1:over_factor*N)=D;
A(:,over_factor*N+1:2*over_factor*N)=D_pos;
%Basis Pursuit on random basis
f=ones(2*nAtoms,1); lb=zeros(2*nAtoms,1);
alpha=linprog(f,[],[],A,y,lb);
figure(2);plot(1:nAtoms,X,’r:’,1:nAtoms,alpha(1:nAtoms),’b–’);
legend(‘original’,'alpha’);
%obtain a few largest component, by sorting alpha
[z,imax]=sort(alpha(1:nAtoms));
%Reconstruct estimated signal
x_hat=zeros(N,1);
for k=1:N_support;
x_hat=x_hat+ D(:,imax(nAtoms-k+1));
end;
x_hat1=D*alpha(1:nAtoms);
figure(3);plot(1:N,x0,’k.’,1:N,y,’b-’,1:N,x_hat,’r:’,1:N,y-x0,’g-’,1:N,x_hat1,’k–’);
legend(‘original’,'noisy’,'estimated’,'noise-level’,'D*alpha’);
figure(4);plot(1:N,y-x0,’r-’,1:N,x_hat-x0,’b-’,1:N,x_hat1-x0,’k-’);
legend(‘initial error’,'estimation error’,'D*alpha error’);

Gabungan pasangan sistem-basis impulse-DCT menghasilkan sebuah kamus-basis berukuran 64×128. Batas sparsity Elad-Bruckstein untuk pasangan basis ini adalah sekitar (0.9).8 atau sekitar 7. Dengan hanya mengambil 4 buah vektor basis penyusun, yaitu vektor pada urutan ke 8, 32, 73, dan 125 haruslah P1 dapat menemukan solusi eksak-nya.

Selected Basis

Gambar 2. Vektor basis terpilih didalam kamus basis

Gambar 2 menunjukkan vektor basis terpilih oleh P1 yang tepat berhimpit dengan basis penyusun sinyal. Dengan demikian seluruh komponen sinyal telah berhasil ditemukan. Disamping itu, koefisian basis lainnya jauh lebih kecil daripada komponen basis penyusun. Komponen ini biasanya dihilangkan dengan peng-ambangan (tresholding).

Sinyal rekonstruksi

Gambar 3. Sinyal asal dan hasil rekonstruksi

Gambar 4. Sinyal kesalahan

Berbagai sinyal dibandingkan dalam Gambar 3. Terlihat bahwa sebenarnya tingkatan sinyal derau (kurva berwarna hijau) cukup tinggi jika dibandingkan dengan sinyal asalnya (titik-titik hitam), menghasilkan distorsi yang signifikan (sinyal warna biru). Namun demikian, karena semua koefisien berhasil ditemukan, sinyal dapat direkonstruksi secara eksak. Hal ini ditunjukkan pada Gambar 4, dimana sinyal kesalahan (biru) berharga nol.

Catatan: Beberapa browser menerjemahkan karakter ‘prime’ secara salah sehingga program yang ada di blog ini tidak bisa langsung di-copy paste dari halaman web ke Command Window Matlab untuk dijalankan. Perbaiki terlebih dahulu kesalahan ini sebelum menjalankan program.

Posted in Compressive Sensing, Signal Decomposition and Search of Basis, Theory and Concept, Uncategorized | Leave a Comment »

Dekomposisi Sinyal dan Pencarian Basis Ideal (3)

Posted by suksmono on May 19, 2008

Pencarian basis untuk pasangan transformasi waktu-frekuensi dapat diperluas untuk sebarang pasangan basis ortonormal. Berbagai macam basis baru yang telah ditemukan seiring maraknya penelitian wavelet dan analisis multiresolusi, baik yang ortonormal (mis. wavelet, Fourier, DCT, Hadamard, ridgelets), yang tak-ortonormal (mis. Gabor, pseudo-polar FFT), maupun transformasi non-kuadratis (Laplacian pyramid, wavelet invarian-geser) perlu dielaborasi supaya bisa dipakai untuk menyatakan berbagai macam bentuk sinyal yang ditemui di dunia nyata. Pencarian vektor basis ”ideal” didalam gabungan berbagai macam basis memerlukan generalisasi sifat keunikan dari representasi sinyal sparse. Tulisan ini memberikan penjelasan secara garis besar hasil-hasil yang telah dicapai oleh Elad-Bruckstein dan perbandingannya dengan formulasi dari Donoho.

Batas Keunikan Elad-Bruckstein dari Pasangan Orto-basis

Pada tulisan terdahulu telah diberikan batas Donoho-Huo (DHB: Donoho-Huo Bound) dari jumlah komponen penyusun sinyal {|| \alpha ||}_{0}  untuk setiap pasangan sistem basis ortonormal (ortobasis) \Theta  dan \Psi  supaya sinyal s=[ \Psi , \Theta ] \alpha  dapat dinyatakan secara unik dan P1 sanggup menemukannya, yaitu

Batas P1 Donoho-Huo: {|| \alpha ||}_{0} < 0.5 (1+{\mu}^{-1})

dimana \mu = \mu (\Psi , \Theta )  adalah koherensi antar basis. Pada kasus pasangan ortobasis ini, Elad dan Bruckstein memberikan batas yang lebih baik, sehingga cukup penting untuk dijelaskan terlebih dahulu.

Sebuah sinyal s sepanjang N-cuplikan dapat dinyatakan dalam sebarang sistem ortobasis, misalnya sebagai s= \Theta {\alpha}_{1}  atau s= \Psi {\alpha}_{2}  . Sebagaimana yang telah disinggung pada tulisan sebelumnya, koherensi antar sistem basis dapat dinyatakan sebagai

Koherensi pasangan basis: 1/\sqrt{N} \leq \mu( \Psi , \Theta) = \max_{\theta \in \Theta , \psi \in \Psi} | < \psi , \theta> | \leq 1  (1)

Jumlah support sinyal s dari kedua kawasan adalah {|| {\alpha}_{1}||}_{0} + {|| {\alpha}_{2}||}_{0}  . Untuk setiap dua bilangan bulat positif x dan y, maka akan berlaku sifat {x}^{2}+{y}^{2} \geq 2xy  . Karena {(x+y)}^{2} = ( {x}^{2} + {y}^{2} ) + 2xy  maka akan berlaku pula { \left( x+y \right)}^{2} \geq 4xy  , sehingga x+y \geq 2 \sqrt{xy}  . Dengan demikian relasi berikut akan berlaku:

Teorema 1. {|| {\alpha}_{1}||}_{0} + {|| {\alpha}_{2} ||}_{0} \geq 2 \sqrt{ {|| {\alpha}_{1} ||}_{0}{|| {\alpha}_{2}||}_{0} } \geq \frac{2}{\mu}  (2)

Hubungan antara dua suku terakhir pada pertidaksamaan diatas, yaitu 2 \sqrt{ {|| {\alpha}_{1} ||}_{0}{|| {\alpha}_{2} ||}_{0} } \geq \frac{2}{\mu}  , dapat dijelaskan dengan analisis sederhana sebagai berikut:

  1. Nilai \sqrt{ {|| {\alpha}_{1} ||}_{0}{|| {\alpha}_{2} ||}_{0} }  minimum terjadi jika keduanya basis impuls yang hanya memiliki satu komponen vektor bernilai tak-nol {\alpha}_{1}  dan {\alpha}_{2}  , sehingga {|| {\alpha}_{1}||}_{0}= {|| {\alpha}_{2} ||}_{0} =1  dan \sqrt{{||{\alpha}_{1} ||}_{0}{|| {\alpha}_{2} ||}_{0}} = 1  . Untuk pasangan basis ini, nilai koherensi mencapai maksimum atau \mu=1  sehingga \sqrt{ {|| {\alpha}_{1} ||}_{0}{|| {\alpha}_{2} ||}_{0} } = 2/ \mu  berlaku.
  2. Disisi lain, nilai {|| {\alpha}_{1}||}_{0} + {|| {\alpha}_{2} ||}_{0}  maksimum terjadi jika semua komponen dari vektor {\alpha}_{1}  dan {\alpha}_{2}  tak nol, misalnya jika keduanya basis Hadamard. Pada saat tersebut {|| {\alpha}_{1} ||}_{0}= {|| {\alpha}_{2} ||}_{0} =N  dan \sqrt{{||{\alpha}_{1} ||}_{0}{||{\alpha}_{2} ||}_{0}} = N  . Nilai koherensi akan mencapai nilai maksimum juga, yaitu \mu=1  . Dengan demikian \sqrt{ {|| {\alpha}_{1} ||}_{0}{|| {\alpha}_{2} ||}_{0} }=N \geq 2/ \mu  juga berlaku.

Akhirnya diperoleh hubungan antara koherensi \mu  dengan batas dari jumlah koefisien dikedua kawasan, yaitu:

Batas P1 Elad-Bruckstein: {|| {\alpha}_{1} ||}_{0} + {|| {\alpha}_{2} ||}_{0} \geq \frac{2}{\mu}  (3)

Batas atas ini tentu lebih baik daripada yang telah diperoleh Donoho-Huo sebelumnya, yang untuk kasus sinyal s ini dapat dituliskan:

Batas P1Donoho-Huo: {|| {\alpha}_{1} ||}_{0} + {|| {\alpha}_{2} ||}_{0} \geq 0.5 \left( 1+ {\mu}^{-1} \right)  (4)

Andaikan bahwa sinyal ternormalisasi (berenergi satu satuan) s dapat dinyatakan dengan dua cara berbeda didalam kamus-basis \Phi = \Psi \cup \Theta  , yaitu s = \Phi {\alpha}_{1}  = \Phi {\alpha}_{2}  . Dengan demikian

0 = \Phi \left( {\alpha}_{1} - {\alpha}_{2} \right) = \begin{pmatrix} \Psi & \Theta \end{pmatrix} \begin{pmatrix} {x}_{1} \\ {x}_{2} \end{pmatrix}

Akibatnya \Psi {x}_{1} + \Theta {x}_{2}=0  atau \Psi {x}_{1} = -\Theta {x}_{2}  . Berdasarkan prinsip ketidakpastian, maka

\frac{2}{\mu} \leq {||{x}_{1}||}_{0} + {||{x}_{2}||}_{0} = {||{\alpha}_{1}-{\alpha}_{2}||}_{0} \leq {||{\alpha}_{1}||}_{0} + {||{\alpha}_{2}||}_{0}  , atau {||{\alpha}_{1}||}_{0} + {||{\alpha}_{2}||}_{0} \geq \frac{\mu}{2}

Ini berarti bahwa dua representasi dari sinyal yang sama tidak akan sekaligus sparse. Akhirnya diperoleh hasil berikut:

Teorema 2. Jika ada suatu representasi yang memenuhi syarat

batas keunikan P0 Elad-Bruckstein: {|| \alpha ||}_{0} < \frac {1}{\mu}  (5)

maka haruslah representasi itu unik atau yang paling sparse.

Sampai dengan tahap ini, pencarian basis ideal \alpha  masih mengandalkan P0. Tahap berikutnya adalah mencari batas sparsity dimana P1 sanggup menemukan solusi yang identik dengan solusi P0.

Untuk suatu sinyal s=[ \Psi ; \Theta] \alpha  , dengan asumsi sparsity pada \alpha  dan jumlah koefisien k1<k2 sehingga vektor

\alpha= \begin{pmatrix} {\alpha}_{1} & {\alpha}_{2} & ... &{\alpha}_{N} & ; & {\alpha}_{N+1} & {\alpha}_{N+2} & ... {\alpha}_{2N} \end{pmatrix}

memiliki k1 buah koefisien tak nol dalam N koefisien pertama (yaitu di dalam \Psi  ) dan k2 buah tak nol untuk N koefisien kedua (didalam \Theta  ), Elad-Bruckstein mendapatkan hasil berikut:

Jika k1 dan k2 memenuhi hubungan 2{\mu}^{2}{k}_{1}{k}_{2} + \mu {k}_{2} - 1 <0  , maka P1 akan dapat menemukannnya.

Hasil yang lebih lemah tapi bentuknya sederhana diberikan oleh {k}_{1} + {k}_{2} < \frac{\left (\sqrt{2} - 0.5 \right)}{\mu}  , sehingga

Batas Keunikan P1 Elad-Bruckstein: {|| \alpha ||}_{0}<\frac{\sqrt{2}-0.5}{\mu} \sim \frac{0.9142}{\mu}  (6)

Sebagai ilustrasi, untuk sinyal s sepanjang N=1024 cuplikan yang tersusun dari komponen pasangan basis impulse-Fourier, maka P0 akan sanggup menemukan solusi untuk 32 vektor basis. Bandingkan dengan P1 yang memberikan harga lebih kecil lagi, yaitu 16 menurut DHB. Untuk kasus pasangan basis ini, EBB memberikan angka untuk P1 yang lebih besar, yaitu 29 buah.

Batas Keunikan P0 dan P1 untuk Kamus Basis

Keunikan untuk sepasang ortobasis telah diturunkan secara langsung dari prinsip ketidakpastian. Sebuah kamus basis bersifat overcomplete, yakni jumlah vektor basisnya jauh lebih banyak dari yang seharusnya cukup untuk membentuk sebuah vektor berdimensi N. Andaikan suatu sinyal s sepanjang N-cuplikan dapat dinyatakan dalam dua representasi didalam kamus-basis \Phi  yang dinyatakan sebagai matriks berdimensi NxK, dimana K>>N s= \Phi {\alpha}_{1} = \Phi {\alpha}_{2} , dengan demikian \Phi ({\alpha}_{1}-{\alpha}_{2}) = 0  atau \Phi v = 0

Hasil terakhir menyatakan kombinasi linier vektor-vektor kolom dari \Phi  yang saling bergantung dan akan dicari sekelompok vektor yang demikian. Untuk keperluan ini akan diperkenalkan konsep spark dari suatu matriks yang dilambangkan sebagai \sigma  sebagai berikut

Definisi. Spark dari suatu matriks \Phi  yang dituliskan sebagai \sigma = spark( \Phi )  adalah bilangan bulat terkecil yang menyatakan banyaknya anggota dari sekelompok vektor kolom dari \Phi  yang saling bebas.

Berikut ini ilustrasi dari konsep nilai spark untuk dua buah matriks yang sederhana.

spark\left[\begin{pmatrix}1 & 0 & ...&0 & 1\\ 0 & 1 & ... &0 & 1\\ ... & ... & ...&... & ... \\ 0 & 0 &0 & 1 & 1 \end{pmatrix} \right]=N+1

sedangkan

spark\left[\begin{pmatrix}1 & 0 & ...&0 & 1\\ 0 & 1 & ... &0 & 0\\ ... & ... & ...&... & ... \\ 0 & 0 &0 & 1 & 0 \end{pmatrix} \right]=2

Apa hubungannya dengan rank? Rank sebuah matriks adalah jumlah maksimum vektor kolom yang saling-bebas, sedangkan spark adalah jumlah minimum vektor kolom yang tidak-saling-bebas. Secara umum hubungan berikut ini akan berlaku:

2 \leq \sigma=spark (\Phi) \leq rank {{}\Phi{}} + 1

Berdasarkan definisi, untuk sinyal s yang diatas akan berlaku {||v||}_{0} \geq \sigma  , dengan demikian ||{\alpha}_{1}-{\alpha}_{2} {||}_{0} \geq \sigma  atau \sigma \leq || {\alpha}_{1}-{\alpha}_{2} {||}_{0} \leq ||{\alpha}_{1} {||}_{0} + ||{\alpha}_{2} {||}_{0}  .Akhirnya untuk kamus basis \Phi  akan diperoleh hubungan yang mirip dengan prinsip ketidakpastian yaitu

||{\alpha}_{1} {||}_{0} + ||{\alpha}_{2} {||}_{0} \geq \sigma  (7)

yang artinya bahwa dua buah representasi berbeda dari suatu sinyal s didalam kamus-basis \Phi  tidak boleh keduanya sekaligus sparse. Selanjutnya akan diperoleh persyaratan keunikan sinyal s sebagai berikut

Teorema 3. Representasi sinyal s=\Phi \alpha  dimana || \alpha {||}_{0} < \sigma/2  akan bersifat unik.

Koherensi suatu kamus-basis \Phi  dapat didefinisikan dengan cara yang mirip seperti pasangan orto-basis sebagai berikut

Koherensi kamus-basis: 2 \leq \mu (\Phi) = {max}_{m \neq n} \left| <{\phi}_{m}, {\phi}_{n}> \right| \leq 1  [8]

Berdasarkan teorema piringan Gershgorin, dapat ditunjukkan bahwa batas bawah (pesimistik) spark adalah \sigma \geq (1+{\mu}^{-1})  . Dengan demikian, dari [8] akan diperoleh

||{\gamma}_{1} {||}_{0} +||{\gamma}_{1} {||}_{0} \geq \sigma \geq (1+{\mu}^{-1})

dan dari Teorema 3 akan diperoleh

Teorema 4. Representasi sinyal s=\Phi \alpha dimana ||\alpha {||}_{0} \leq 0.5( 1+ {\mu}^{-1} ) \leq \sigma/2  akan bersifat unik.

Dengan cara yang sama seperti sebelumnya, keunikan sinyal s didalam kamus basis \Phi  terhadap optimasi P1 dapat dinyatakan sbb;

Teorema 5. Diberikan suatu sinyal s dengan representasi s=\Phi \alpha  , jika ||\alpha {||}_{0} < 0.5( 1+{\mu}^{-1})  maka P1 atau algoritma Basis Pursuit akan dapat menemukannya.

Apa bedanya dengan yang diperoleh Donoho-Huo pada tulisan sebelumnya? Perbedaan pertama, jumlah support dari sinyal yang dibentuk pasangan orto-basis dalam penurunan Elad-Bruckstein lebih tinggi dari yang diberikan Donoho-Huo. Yang kedua, Elad-Bruckstein memberikan generalisasi dari pasangan basis ortogonal ke kamus-basis dengan batas sparsity yang diberikan oleh (9).

Contoh Permasalahan

Soal-1: Tentukan koherensi antara ortobasis DCT dan ortobasis acak tersebar normal, masing-masing berdimensi 512×512. Tentukan batas sparsity untuk P1 berdasarkan DHB (Donoho-Huo Bound) dan EBB (Elad-Bruckstein Bound).

Jawab: Perhatikan bahwa sistem yang ditinjau adalah dua orto-basis, bukan kamus-basis. Dengan demikian, perhitungan koherensi akan memakai rumus (2) sedangkan batas sparsity adalah rumus (1) untuk DHB dan rumus (6) untuk EBB. Kode Matlab berikut dapat dipakai.

% progam menghitung koherensi dan batas sparsity
% antara basis DCT dengan basis acak Gaussian
clear;close all;clc;
N=512; %panjang sinyal yang diasumsikan
%buat basis DCT berdimensi NxN
THETA=dct(eye(N,N));
%buat basis acak Gaussian NxN
PSI=orth(randn(N,N)’)';
%hitung koherensi antara DCT dengan Gauss
mu_DCT_Gauss=mu_THETA_PSI(THETA,PSI);
disp(sprintf(‘koherensi DCT-Gauss = %f\n’, mu_DCT_Gauss));
%Tentukan batas sparsity
%kasus-1: DHB (Donoho-Huo bound)
N_DHB=floor(0.5*(1+1/mu_DCT_Gauss));
disp(sprintf(‘batas sparsity Donoho-Huo= %d\n’, N_DHB));
%kasus-2: EBB (Elad-Bruckstein bound)
N_EBB=floor((sqrt(2)-0.5)/mu_DCT_Gauss);
disp(sprintf(‘batas sparsity Elad-Bruckstein= %d’, N_EBB));
%==FUNGSI ===
%fungsi menghitung koherensi dua orto basis
function [mu]=mu_THETA_PSI(THETA,PSI);
[M,N]=size(PSI); %THETA diasumsikan berdimensi sama
CORR=transpose(conj(THETA))*PSI;
mu=max(max(CORR));


Hasil perhitungan memberikan nilai koherensi \mu =0.219030  , DHB =2, dan EBB = 4. Terlihat perbedaan cukup signifikan antara EBB dengan DHB.

Soal-2: Tentukan tingkat koherensi dari sebuah kamus-basis acak tersebar normal yang berdimensi 64×1024 dan hitunglah batas sparsity dari sinyal s sepanjang 64 cuplikan yang terbentuk dari kamus-basis ini!

Jawab: Didalam compressive sensing, basis acak dapat dianggap sebagai basis universal karena bersifat hampir ortogonal ke semua sistem basis lain. Koherensi dihitung menurut rumus (7) sedangkan batas sparsity memakasi rumus (9). Program Matlab berikut ini dapat dipakai untuk membantu pemecahan masalah ini.

% progam menghitung koherensi dan
% batas sparsity basis acak Gaussian
clear;close all;clc;
N=64;M=1024; %dimensi kamus basis
%buat basis acak Gaussian MxN
D=orth(randn(N,M)’)';
%hitung koherensi kamus basis D
mu=mu_PHI(D);
disp(sprintf(‘koherensi kamus-basis = %f\n’, mu));
%Tentukan batas sparsity
N_sparse=floor(0.5*(1+1/mu));
disp(sprintf(‘batas sparsity = %d’, N_sparse));
% === fungsi yang dipakai ===
%fungsi menghitung koherensi kamus basis
function [mu]=mu_PHI(PHI);
[M,N]=size(PHI);
CORR=transpose(conj(PHI))*PHI;
mu=max(max(CORR));

Hasil running program Matlab diatas menghasilkan nilai koherensi \mu \sim 0.094  dan batas sparsity sebesar 5. Dengan demikian, sebarang sinyal sepanjang 64 cuplikan yang terbentuk dari kombinasi linier kamus basis berukuran 64×1024 akan dapat ditemukan dengan optimasi P1 atau algoritma Basis Pursuit.

Referensi

  1. M.Elad and A.M. Bruckstein, “A generalized uncertainty principle and sparse representation in pairs of bases”, IEEE Trans. on Information Theory, Vol. 48, pp. 2558-2567, September 2002.
  2. D. L. Donoho and M. Elad, “Maximal sparsity representation via l1 minimization”, Proc. of the National Academy of Science, Vol. 100, pp. 2197-2202, March 2003.
  3. S.S. Chen, D. Donoho, and MA Saunders, ”Atomic decomposition by basis pursuit,” SIAM J. on Scientific Computing, 20(1):33-61, 1999.

Posted in Compressive Sensing, Signal Decomposition and Search of Basis, Theory and Concept, Uncategorized | Leave a Comment »

Dekomposisi Sinyal dan Pencarian Basis Ideal (2)

Posted by suksmono on May 16, 2008

Penggabungan berbagai basis ortogonal menjadi sebuah kamus basis-lewat-lengkap (overcomplete) diharapkan dapat menaikkan efisiensi dekomposisi sinyal. Peningkatan kinerja sistem kompresi terjadi karena dua hal: (1) tingkat kompresi tinggi karena koefisien yang diambil hanya sedikit, dan (2) sejumlah kecil koefisien tadi sudah cukup untuk merekonstruksi sinyal secara baik. Namun demikian pencarian vektor didalam basis overcomplete dengan kriteria cacah koefisien terkecil memerlukan sumberdaya komputasi yang sangat besar. Melalui prinsip ketidakpastian, masalah ini dapat diubah menjadi optimisasi konveks yang lebih mudah dipecahkan. Syaratnya, sinyal harus bersifat sparse.

Tinjau gabungan M buah sistem basis ortogonal {\Phi}_{1}, {\Phi}_{2}, ...,{\Phi}_{M}  yang masing-masing berukuran N×N menjadi satu buah basis overcomplete \Phi={\Phi}_{1} \cup {\Phi}_{2} ...\cup {\Phi}_{M}  . Basis baru \Phi  ; yang juga disebut sebagai dictionary (kamus basis), akan memiliki dimensi N×MN. Suatu sinyal s=[s1 s2 ... sN]T sepanjang N cuplikan dapat dinyatakan dengan berbagai macam kombinasi linier vektor-vektor basis yang ada didalam \Phi  , bahkan dari masing-masing sistem basis ortogonal {\Phi}_{m}  karena bersifat lengkap. Dengan pertimbangan efisiensi yang merupakan tuntutan kompresi sinyal, haruslah dipilih kombinasi yang berbentuk:

s=\sum_{\gamma} {\alpha}_{\gamma}{\phi}_{\gamma}  (1)

dimana {\phi}_{\gamma}  adalah vektor basis didalam \Phi  . Pada penjumlahan tersebut \gamma = (m,i)  menyatakan bahwa komponen tersebut adalah vektor basis ke–i didalam sistem basis {\Phi}_{m}  . Pertimbangan efisiensi mengharuskan cacah {\alpha}_{\gamma}  , yaitu

{||\alpha||}_{0} = \sum_{\forall \gamma} {|{\alpha}_{\gamma}|}_{0}  (2)

sekecil mungkin. Sebagai ilustrasi arti norm diatas , tinjau contoh perhitungan berikut ini.

Suatu dictionary \Phi  yang merupakan hasil penggabungan 2 buah basis ortogonal memiliki dimensi 4×2, dimana
{\Phi}_{1}=\begin{pmatrix}1 & 0 \\ 0 & 1\end{pmatrix}  , {\Phi}_{2}= \frac{1} {\sqrt{2}} \begin{pmatrix}1 & 1 \\ 1 & -1\end{pmatrix}  , dan \Phi=\begin{pmatrix}1 & 0 & \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \\ 0 & 1 & \frac{1}{\sqrt{2}} & -\frac{1}{\sqrt{2}} \end{pmatrix}

Tiga dari empat basis didalam \Phi  ini dapat dipakai untuk menyatakan sebuah sinyal yang panjangnya dua cuplikan, misalnya {s=\begin{pmatrix} -0.75/\sqrt{2} &  1-0.75/\sqrt{2} \end{pmatrix}}^{T} . Sinyal ini dapat terbentuk dari \gamma \in \Gamma= \begin{Bmatrix}(1,2), (2,1), (2,2) \end{Bmatrix}  , masing-masing dengan koefisien {\alpha}_{(1,2)}=0.25, {\alpha}_{(2,1)}=-0.25, dan \: {\alpha}_{(2,2)}= -0.5  . Dengan demikian {\alpha}_{(1,2)}=0.25, {\alpha}_{(2,1)}=-0.25, {\alpha}_{(2,2)}= -0.5   , sehingga
{||\alpha||}_{0}={0.25}^{0}+{0.25}^{0}+{0.5}^{0}=1+1+1=3
Terlihat bahwa norm diatas hanya menghitung cacah dari koefisien terpilih.

Pencarian basis representasi sinyal s didalam \Phi  dengan cacah terkecil disebut sebagai masalah (P0) yang dapat dinyatakan sebagai optimisasi berikut:

({P}_{0}): min {||\alpha||}_{0}, \: \; s.t.\; \; s=\sum_{\gamma} {\alpha}_{\gamma}{\varphi}_{\gamma} (3)

Mengingat ada sejumlah tak-berhingga solusi dari sistem persamaan linier yang underdetermined, apakah nilai minimum ini dapat ditemukan? Untuk pasangan basis waktu-frekuensi, Donoho memberikan batasan supaya pencarian menemukan solusi yang unik dengan menggunakan prinsip ketidakpastian sebagai berikut.

Teorema 1. Jika s(n) adalah sinyal waktu diskrit sepanjang N cuplikan yang memiliki Nt buah support (koefisien tak-nol), sedangkan S(k) adalah transform Fourier sinyal tersebut yang memiliki Nw buah support, maka {N}_{t} {N}_{\omega} \geq \sqrt{N}  atau
{N}_{t}+{N}_{\omega} \geq 2 \sqrt{N}  (4)

Jika cuplikan yang bernilai tak-nol (support) dari sinyal di kedua wilayah kurang dari setengah batas prinsip ketidakpastian diatas, maka P0 akan memiliki solusi yang unik.

Teorema 2. Tinjau sinyal s(n) yang tersusun dari himpunan T vektor basis kawasan waktu dan himpunan \Omega  basis dari kawasan frekuensi. Jika
|T|+|\Omega| \leq \sqrt{N}  (5)
maka P0 akan memiliki solusi yang unik. Sementara itu akan ada sinyal s(n) lain sedemikian hingga |T|+|\Omega| = \sqrt{N}  dan P0 akan memiliki solusi yang tidak unik.

Pada teorema ini, notasi |.| menyatakan kardinalitas atau banyaknya anggota himpunan. Donoho menunjukkan hasil diatas dengan meninjau kasus fungsi ekstremal sisir Dirac (Dirac comb) yang telah dijelaskan sebelumnya. Untuk N bilangan kuadrat sempurna (perfect square), misalnya 4, 9, 16, … dst, sisir Dirac memiliki support sebesar \sqrt{N}  di kawasan waktu maupun kawasan frekuensi. Ini berarti bahwa sebuah sisir Dirac dapat dinyatakan sebagai \sqrt{N}  buah vektor basis impuls atau \sqrt{N}  buah basis Fourier, sehingga representasinya didalam basis gabungan impuls-Fourier tidaklah unik. Akibatnya solusi dari P0 juga tidak unik. Sebaliknya, jika jumlah support-nya kurang dari K= \sqrt{N}  , maka akan hanya ada satu buah solusi saja, sehingga solusi P0 akan bersifat unik.

Meskipun solusinya unik, memecahkan masalah P0 bukanlah hal yang mudah. Untuk mengatasi hal ini, Chen dan Donoho mengajukan teknik Basis Pursuit, yaitu dengan menggantikan P0 dengan optimisasi P1 sebagai berikut:

({P}_{1}): min {||\alpha||}_{1}, \: \; s.t.\; \; s=\sum_{\gamma} {\alpha}_{\gamma}{\varphi}_{\gamma} (6)

Keberadaan solusi P1 memerlukan batasan support sinyal didalam \Phi  yang lebih kecil lagi, yaitu kurang dari setengah batas P0. Untuk basis waktu-frekuensi, ketentuan ini dinyatakan dalam teorema berikut.

Teorema 3. Untuk sinyal s seperti pada Teorema 2, jika
|T| + |\Omega| \leq \frac{1}{2}\sqrt{N}
maka solusi dari P1 adalah unik, yang juga merupakan solusi unik untuk P0. Sementara itu, ada sinyal yang jika |T| + |\Omega| = \frac{1}{2}\sqrt{N}  dimana P1 memiliki solusi yang tidak unik.

Sinyal dengan support dibawah batas ketentuan P1 ini disebut sebagai sinyal sparse karena dapat dibentuk hanya dari sejumlah kecil vektor basis didalam \Phi  . Secara singkat dapat dirangkum sebuah hasil penting berikut ini:

Jika suatu sinyal s bersifat sparse dalam kawasan waktu-frekuensi, maka dekomposisinya bersifat unik dan algoritma Basis Pursuit atau P1 akan dapat menemukannya.

Bagaimana dengan pasangan basis selain waktu-frekuensi? Pada tulisan sebelumnya dikatakan bahwa prinsip ketidakpastian di sebarang kawasan telah diperumum oleh Elad dan Bruckstein dengan memasukkan faktor koherensi antara sistem basis \mu ({\Phi}_{1}, {\Phi}_{2})  yang secara konsisten akan kembali menjadi WHUP jika kedua basis adalah pasangan waktu-frekuensi. Dasar pemikiran ini dipakai oleh Elad-Bruckstein maupun Donoho-Huo untuk menghasilkan teorema keunikan dari solusi P1 berikut.

Teorema 4. Andaikan \Phi= {\Phi}_{1} \cup {\Phi}_{2}  adalah gabungan kedua basis. Jika suatu sinyal waktu diskrit s dapat dinyatakan sebagai s=\Phi \alpha  dengan \alpha  memenuhi syarat
{||\alpha||}_{0} \leq \frac {1}{2} \left( 1+ \frac {1} {\mu ({\Phi}{1}, {\Phi}{2})} \right)  (7)
maka \alpha  adalah solusi unik dari P1 dan juga sekaligus merupakan solusi unik dari P0.

Dalam bahasa yang kurang-teknis, teorema ini mengatakan bahwa suatu sinyal s yang terbentuk dari kombinasi linier vektor-vektor basis didalam kamus-basis \Phi  dapat ditemukan vektor basis penyusunnya dengan memecahkan masalah optimisasi (P1), asalkan jumlah vektor basis penyusun sinyal ini kurang dari nilai tertentu.

Seberapa besar signifikansi perubahan optimisasi P0 menjadi P1? Meskipun sekilas optimisasi ini hanya mengubah norm dari l0 menjadi l1, perubahan yang terjadi sangat besar. Saat ini, kebanyakan algoritma optimisasi didalam pengolahan sinyal justru menggunakan kriteria atau fungsi biaya kuadratis yang tak lain adalah l2. Penjelasan diatas bahkan menginginkan pencarian solusi ideal dengan kriteria l0. Dari sudut pandang pencarian ruang solusi, optimisasi dengan kriteria l2 dan l1 bersifat konveks, sedangkan l0 bukan konveks. Algoritma untuk mencari solusi dari optimisasi konveks yang efisien telah banyak tersedia. Transformasi dari masalah l0 menjadi l1 ini disebut juga sebagai proses konveksifikasi.

Salah satu contoh sederhana didalam statistik mengenai berbagai norm ini adalah perhitungan mencari nilai tengah dari populasi. Meminimumkan nilai kesalahan kuadrat terkecil, yang tak lain adalah l2, akan menghasilkan nilai mean cuplikan, sedangkan jika yang diminumkan adalah nilai absolut kesalahan atau l1, solusinya akan berupa median. Pada pengolahan sinyal atau citra dijital, para peneliti menemukan bahwa tapis median lebih kokoh (robust) terhadap derau jika dibandingkan dengan tapis mean. Analogi yang sama dapat dipakai untuk menjelaskan mengapa P1 ini memiliki peranan yang sangat penting dalam pencarian basis sinyal yang ideal.

Referensi

[1]. M. Elad and A.M. Bruckstein, “A generalized uncertainty principle and sparse representation in pairs of basis,” IEEE Trans. Info. Theory, 49. 2558-3567.

[2]. S.S. Chen, D. Donoho, and MA Saunders, ”Atomic decomposition by basis pursuit,” SIAM J. on Scientific Computing, 20(1):33-61, 1999.

[3]. DL Donoho & PB Stark, ”Uncertainty principles and signal recovery”, SIAM J. on Applied Math., 49(3):906-31, 1989.

[4]. DL Donoho & X. Huo, Uncertainty Principles and Ideal Atomic Decomposition, Tech. Report, Stanford University.

Posted in Compressive Sensing, Signal Decomposition and Search of Basis, Theory and Concept, Uncategorized | 3 Comments »

Dekomposisi Sinyal dan Pencarian Basis Ideal (1)

Posted by suksmono on May 14, 2008

Pemanfaatan transformasi ortogonal untuk kompresi sinyal telah ditunjukkan pada tulisan sebelumnya. Secara intuitif dapat difahami bahwa sistem kompresi menghendaki terkumpulnya energi di kawasan transformasi hanya ke sejumlah kecil koefisien supaya tingkat kompresi nya tinggi dengan kualitas hasil rekonstruksi yang terjaga. Namun demikian hal ini tidak selalu bisa dilakukan karena bagaimanapun setiap sistem basis hanya cocok untuk segolongan sinyal-sinyal tertentu saja. Basis sinusoid dalam transformasi Fourier atau DCT melakukan pengkompakan energi dengan baik jika sinyal asalnya tersebar di kawasan waktu. Sebaliknya, sinyal yang terlokalisir dalam kawasan ini sebaiknya memakai basis impuls atau wavelet.

Impulse Base

Gambar 1. Beberapa vektor basis impuls

DCT Base

Gambar 2. Beberapa vektor basis DCT

Gambar diatas memperlihatkan beberapa vektor basis dari sistem basis impuls \Phi  dan sistem basis DCT \Psi  . Sifat lokalisasi sinyal terlihat dengan jelas pada kedua gambar; basis impuls terlokalisir sedangkan DCT tersebar. Program Matlab berikut dapat dipakai untuk menampilkan vektor-vektor basis seperti diatas.

%———————————————–
% Sistem Basis Impuls dan DCT
%———————————————–
N=64; %dimensi vektor basis
x=1:N; %buat sumbu x
PHI=eye(N,N); %Basis Impulse
PSI=dct(eye(N,N)); %Basis DCT
figure(1);plot(x,PHI(2,:),x,1+PHI(4,:),x,2+PHI(8,:),…
x,3+PHI(16,:),x,4+PHI(32,:),x,5+PHI(60,:));
title(‘Vektor basis impuls’);
figure(2);plot(x,PSI(2,:),x,0.5+PSI(4,:),x,1+PSI(8,:), …
x,1.5+PSI(16,:),x,2+PSI(32,:),x,2.5+PSI(60,:));
title(‘Vektor basis DCT’);

Suatu sinyal biasanya merupakan campuran dari komponen yang terlokalisir dan yang tersebar. Basis mana yang harus dipilih dalam kasus yang demikian? Tinjau suatu sinyal s(n) yang terbentuk dari sistem basis \Phi  berupa vektor basis ke-4 dan ke-32 dan dua buah lainnya dari sistem basis \Psi , yaitu vektor basis ke-8 dan ke-60. Prinsip kompresi sederhana yang sudah dijelaskan dapat dicobakan untuk setiap sistem basis, yaitu dengan mengalikan matriks transformasi-nya dengan vektor s(n),

Dekomposisi sinyal ke basis impuls: {S}_{1}= \Phi .s

Dekomposisi sinyal ke basis DCT: {S}_{2}=\Psi.s

Langkah berikutnya, ambil beberapa koefisien dengan energi (magnitudo) terbesar. Tentu cara lain dapat dilakukan, misalnya tresholding dengan mencari koefisien-koefisien yang energi/magnitudo-nya berada diatas ambang nilai tertentu. Untuk indeks dari koefisien <em>k </em></span>\<span>in K</span><em><span> </span></em><span> & </span><span>vektor basis terpilih, rekonstruksi dilakukan dengan cara menjumlahkan vektor dari basis yang telah diboboti, yaitu:</span></p> <p class=\"MsoNormal\" style=\"text-align:justify;\">Rekonstruksi dengan basis impuls: latex {\hat{s}}_{1}(n)=\sum_{k \in K}{S}_{1}(k) \Phi(k) &s=-2 $

Rekonstruksi dengan basis DCT: {\hat{s}}_{2}(n)=\sum_{k \in K}{S}_{s}(k) \Psi(k)

Kinerja dari rekonstruksi sinyal dapat dilihat dari energi sinyal kesalahan (error) rata-rata atau RMSE (Root Mean Square Error) masing-masing, yang dinyatakan sebagai:

RMSE=\frac{1}{N}\sqrt{\sum_{n=0}^{N-1}{e(n)}^{2}}   , dimana e(n)=s(n)-\hat{s}(n)

Semakin kecil nilai ini, maka secara obyektif hasil rekonstruksi semakin baik. Nilai RMSE yang besar menandakan terjadinya distorsi signifikan pada sinyal.

Original Signal

Gambar 3. Sinyal asli s(n)

Ordered-magnitude

Gambar 4. Magnitudo koefisien terurut

Gambar 3 menunjukkan sinyal asal s(n) berupa gabungan dua sinusoid dan dua sinyal impuls. Sinyal ini akan didekomposisi dengan basis impuls dan kemudian dengan basis DCT. Dekomposisii sinyal ini kedalam sistem basis tersebut akan menghasilkan koefisien yang dinyatakan sebagai vektor S1 dan S2. Untuk mengamati sebaran energi, magnitudo atau nilai mutlak S1 dan S2 diurutkan dan kemudian dibuatkan kurvanya. Kedua kurva diberikan pada Gambar 4. Terlihat bahwa, meskipun beberapa koefisien awal dominan, penurunan nilai energi koefisien lainnya tidaklah terlalu tajam. Dengan demikian, pengambilan sejumlah kecil komponen dominan akan tetap menghasilkan sinyal kesalahan yang cukup tinggi energinya. Ini artinya rekonstruksi sinyal kurang berhasil.

reconstructed signals

Gambar 5. Perbandingan sinyal asli dengan hasil rekonstruksi

Error Signal

Gambar 6. Sinyal kesalahan

Gambar 5 memperlihatkan hasil rekonstruksi dengan menggunakan delapan buah vektor basis. Kurva biru merupakan hasil dari basis impuls, sedangkan kurva merah dari basis DCT. Perbandingan dengan sinyal asli yang dinyatakan sebagai titik-titik berwarna hitam memperlihatkan beda yang cukup signifikan. Perbedaan ini semakin jelas setelah sinyal kesalahan e1(n) dan e2(n) diperlihatkan pada Gambar 6.

RMSE Curves

Gambar 7 Kinerja RMSE

Kinerja rekonstruksi dapat dievaluasi dengan mengukur RMSE berbagai jumlah vektor basis seperti yang diperlihatkan pada Gambar 7. Penurunan RMSE terjadi secara perlahan dan mencapai nol ketika seluruh N vektor basis diikutsertakan; yang artinya tidak ada kompresi. Program Matlab berikut dapat dipergunakan untuk menunjukkan proses proyeksi dan rekonstruksi diatas.

% Pencarian basis
clear;%close all;clc;
N=64; %dimensi vektor basis
x=1:N; %buat sumbu x
PHI=eye(N,N); %Basis Impulse
PSI=dct(eye(N,N)); %Basis DCT
%bentuk sinyal dari basis DCT saja
s1=zeros(N,1);s1=zeros(N,1);
s1=(PHI(4,:)+PHI(32,:))’;
s2=(PSI(8,:)+PSI(60,:))’;
s=zeros(N,1); s=s1+s2;
figure(1);plot(x,s); title(’sinyal kombinasi impuls dan DCT’);
xlabel(‘n’);ylabel(’s(n)’);
S1=PHI*s; S2=PSI*s; %cari koefisien
%urutkan koefisien dan ambil dua terbesar
[S1_oval,S1_oidx]= sort(abs(S1),’descend’);
[S2_oval,S2_oidx]= sort(abs(S2),’descend’);
% rekonstruksi sinyal
s_hat1=zeros(N,1); s_hat2=zeros(N,1);
for k=1:8; s_hat1=s_hat1+(S1(S1_oidx(k))*PHI(S1_oidx(k),:))’; end;
for k=1:8; s_hat2=s_hat2+(S2(S2_oidx(k))*PSI(S2_oidx(k),:))’; end;
%hitung kesalahan
e1=s-s_hat1;e_rms1=sqrt(e1′*e1)/N;
disp(sprintf(‘rms error-Impulse %f’,e_rms1));
e2=s-s_hat2;e_rms2=sqrt(e2′*e2)/N;d
isp(sprintf(‘rms error-DCT %f’,e_rms2));
figure(2); plot(x,S1_oval,’r-’, x,S2_oval,’b-’);
title(‘magnitudo koefisien terurut’);
xlabel(‘index’);ylabel(‘|S|-terurut’);
legend(‘proyeksi Impuls’,'proyeksi DCT’);
figure(3);plot(x,s,’k.’,x,s_hat1,’b:’,x,s_hat2,’r-’);
title(‘perbandingan hasil rekonstruksi’);
xlabel(‘n’);ylabel(’s(n)’); legend(‘asli’,'rek. impuls’,'rek. DCT’);
figure(4);plot(x,e1,’r-’,x,e2,’b–’);title(’sinyal kesalahan’);
xlabel(‘n’);ylabel(‘e(n)’);
legend(’sinyal error impulse’,’sinyal error DCT’);
%data error beberapa nilai k
x_e=[2 4 8 16 32 64]‘;
e_PHI=[2.18e-2 2.05e-2 1.78e-2 1.33e-2 5.94e-3 0.0]‘;
e_PSI=[2.16e-2 2.02e-2 1.77e-2 1.32e-2 6.28e-3 0.0]‘;
figure(5);plot(x_e,e_PHI,’r-’,x_e,e_PSI,’b-’); title(‘kinerja RMSE pengkodean’);
xlabel(‘jumlah koefisien’);ylabel(‘RMS Error’); legend(‘Impulse’,'DCT’);

Melihat hasil yang kurang baik saat menggunakan satu buah basis ortogonal saja, apakah mungkin kinerja-nya meningkat seandainya dua atau lebih sistem basis digabungkan? Memang cara inilah yang paling masuk akal, akan tetapi cara ini menimbulkan satu masalah besar lainnya. Gabungan kedua basis akan menghasilkan suatu basis baru yang dimensinya menjadi Nx2N. Dilihat dari sistem persamaan linier, permasalahan menjadi lebih sulit dipecahkan karena kita harus mencari solusi dari N buah persamaan dengan 2N peubah yang tak diketahui . Sistem persamaan linier yang demikian disebut sebagai underdetermined.

Sebagai ilustrasi, tinjau suatu matriks transformasi A berukuran 2×2 yang bisa dianggap sebagai transformasi vektor pada bidang atau transformasi sinyal 2 cuplikan x=[x1 x2] menjadi y=[y1 y2]. Pencarian koefisien transformasi x(n) yang tepat dari sinyal teramati y(n) dapat dinyatakan sebagai persamaan matriks: Ax=y dengan x tak diketahui atau

\begin{pmatrix}{a}_{11} & {a}_{12} \\ {a}_{21} & {a}_{22} \end{pmatrix}\begin{pmatrix}{x}_{1}\\ {x}_{2} \end{pmatrix}= \begin{pmatrix} {y}_{1} \\ {y}_{2}\end{pmatrix}

Bentuk sistem persamaan liniernya adalah

{a}_{11}{x}_{1}+{a}_{12}{x}_{2}={y}_{1}

{a}_{21}{x}_{1}+{a}_{22}{x}_{2}={y}_{2}

Karena yang tak-diketahui ada dua {x1, x2} dan ada dua persamaan, jika matriks A tak-singular tau memiliki inverse, maka sistem persamaan linier diatas dapat dipecahkan dengan eliminasi Gauss atau mengalikan inverse A dengan vektor [y1 y2 ]T dari kiri.

Sebaliknya, penggabungan sistem basis A dengan basis baru B menjadi D=[A;B] akan membuat persamaan matriks diatas menjadi

\begin{pmatrix}{a}_{11} & {a}_{12} & {b}_{11} & {b}_{12} \\ {a}_{21} & {a}_{22} & {a}_{21} & {a}_{22} \end{pmatrix} \begin{pmatrix}{x}_{1}\\ {x}_{2}\\ {x}_{3}\\{x}_{4} \end{pmatrix}=\begin{pmatrix} {y}_{1}\\ {y}_{2} \end{pmatrix}

Atau dinyatakan sebagai sistem persamaan linier berikut:

{a}_{11}{x}_{1}+{a}_{12}{x}_{2}+{b}_{11}{x}_{3}+{b}_{12}{x}_{4}={y}_{1}

{a}_{21}{x}_{1}+{a}_{22}{x}_{2}+{b}_{21}{x}_{3}+{b}_{22}{x}_{4}={y}_{2}

Sistem persamaan linier ini terdiri dari empat buah variabel, yaitu {x1, x2, x3, x4} sedangkan jumlah persamaan hanya ada dua. Karena masing-masing sistem basis bersifat lengkap, akan ada banyak sekali kombinasi yang dapat muncul untuk menyatakan suatu vektor y. Pencarian basis ke seluruh ruang solusi memerlukan sumberdaya komputasi yang sangat besar karena termasuk permasalahan kombinatorial. Bagaimana cara mengatasinya?

Posted in Compressive Sensing, Signal Decomposition and Search of Basis, Theory and Concept, Uncategorized | Leave a Comment »

Prinsip Ketidakpastian dan Kompresi Data (3)

Posted by suksmono on May 12, 2008

Pada kawasan diskrit, WHUP dinyatakan sebagai penjumlahan RoS (Region of Support) dari waktu- dan frekuensi- diskrit sinyal. Prinsip ini dapat diperumum untuk sebarang pasangan kawasan melalui hubungan matriks transformasi atau himpunan fungsi basis-nya. Batas daerah dukung ditentukan oleh tingkat koherensi antar basis ini.

Prinsip Ketidakpastian untuk Kawasan Waktu-Frekuensi Diskrit

Prinsip ketidakpastian Weyl-Heisenberg (WHUP) melarang sebuah sinyal kontinyu untuk terlokalisir pada kawasan waktu dan frekuensi sekaligus. Derajat lokalisasi suatu sinyal dinyatakan sebagai lebar/luas wilayah dimana sinyal ini memiliki nilai yang tidak nol, atau daerah-dukung (RoS-region of support) dari fungsi atau sinyal tersebut. Sebagai contoh, fungsi delta Dirac \delta\left(t \right)  hanya memiliki nilai tak nol pada waktu t=0; dengan demikian RoS-nya adalah satu. Sebaliknya, sinyal sinusoidal Asin\left(\omega t+\theta \right)  terbentang ke seluruh kawasan waktu sehingga memiliki RoS di kawasan waktu yang tak berhingga.

Pengertian RoS untuk sinyal (waktu-) diskrit mirip dengan penjelasan diatas, yaitu daerah atau sekumpulan indeks waktu n dimana sinyal diskrit s(n) bernilai tak nol. Sinyal diskrit s(n) sepanjang N cuplikan akan memiliki transform Fourier

S\left( k \right)={1}_{\sqrt N}\sum_{n=0}^{N-1} s \left( n \right) {e}^{-j\frac{2 \pi kn}{N}}  ; dimana k = 0, 1, 2, …, N-1.

Untuk sinyal ini, WHUP menyatakan bahwa penjumlahan RoS (yang dilambangkan sebagai supp) sinyal di kawasan waktu dan RoS di kawasan frekuensi dibatasi oleh suatu nilai tertentu, yakni:

\left|supp\left(s \right) \right|+\left|supp\left(S \right) \right|\geq 2\sqrt{N}

Nilai terkecil dari jumlah kedua RoS ini, yaitu 2 \sqrt{N}  , dicapai oleh fungsi Dirac Comb yang support-nya bersifat invarian terhadap transformasi Fourier. Fungsi ini dinyatakan sbb:

s\left(n \right)=\left\{\begin{matrix} 1 & ; n=m \sqrt{N} ;m=0, 1, 2, ..., \sqrt{N}-1 \\ 0 & ; lainnya \end{matrix}\right

Gambar berikut ini memperlihatkan sinyal s(n) dan transform Fouriernya untuk kasus N=16.

Dirac Comb

Sifat ini menunjukkan bahwa batas bawah dari prinsip ketidapkastian adalah ketat (tight), dan telah dipenuhi oleh sebuah fungsi ekstremal berupa Dirac comb. Kedua gambar sinyal diatas dapat dibuat dengan memanfaatkan script Matlab berikut ini:

%————————————————–
%Sifat Invarian Sisir Dirac (Dirac Comb)
%terhadap Transformasi Fourier
%—————————————————
clear; clc;
N=16; s=zeros(N,1);
sqrt_N = round(sqrt(N));
maxIdx = floor(N/sqrt_N);
for m=1:maxIdx; s(m*sqrt_N) =1; end
figure(1);stem(s);title(‘Time-domain Dirac Comb’);
xlabel(‘Time’);ylabel(‘Amplitude’);
%Dirac Comb in Frequency Domain
S=fft(s)/sqrt(N);
figure(2);stem(abs(S)); title(‘Frequency domain abs(Dirac Comb)’)
xlabel(‘Frequency’);ylabel(‘Magnitude’);

Prinsip Ketidakpastian untuk Sebarang Kawasan Diskrit

Apakah prinsip ketidakpastian berlaku juga untuk pasangan kawasan selain waktu-frekuensi? Seperti telah diuraikan pada tulisan sebelumnya, suatu sinyal dapat dilihat dari berbagai kawasan melalui transformasi. Transformasi pada dasarnya adalah perubahan basis dan operasi untuk sinyal diskrit dapat dilihat sebagai perkalian antara vektor sinyal dengan matriks transformasi \Phi  .

GUP Diagram

Suatu matriks transformasi \Phi  akan berisi vektor-vektor basis lengkap yang membentang ruang vektor. Ini berarti bahwa setiap sinyal dalam kawasan ini dapat dilihat sebagai kombinasi linier dari vektor-vektor basis tersebut. Untuk sepasang basis ortonormal (saling tegak-lurus dan panjangnya satu) \Phi  dan \Psi  , Elad dan Bruckstein memperkenalkan sebuah prinsip ketidakpastian yang diperumum atau GUP (Generalized Uncertainty Principle):

\left| {\Gamma }_{1} \right| + \left| {\Gamma }_{2} \right| \geq \frac{2}{\mu \left(\Phi ,\Psi  \right)}

dimana:

\mu \left(\Phi ,\Psi  \right)=\max_{\phi \in \Phi , \psi \in \Psi }\left| \left< \Phi,\Psi \right> \right|

Pada ekspresi diatas, {\Gamma}_{1}  adalah support dari sinyal dalam basis (kawasan) \Phi  , sedangkan, {\Gamma}_{2}  adalah support dari sinyal dalam basis (kawasan) \Psi  . Besaran \mu  disebut koherensi (mutual coherence) antara kedua basis dan haruslah memenuhi syarat

1/\sqrt N \leq \mu \left( \Phi , \Psi \right) \leq 1

Koherensi menggambarkan tingkat keserupaan kedua basis: nilainya akan kecil jika keduanya berbeda dan akan bernilai satu jika mereka identik.

Sebagai contoh, tinjau kasus WHUP kawasan diskrit yang merupakan kasus khusus dari GUP. Sinyal s(n) berada pada kawasan waktu memiliki pasangan S(k) pada kawasan frekuensi. Bisa dianggap bahwa s(n) memiliki basis berupa matriks identitas I, sedangkan S(k) memiliki basis Fourier diskrit. Dengan demikian bisa dianggap bahwa pasangan basis ini adalah

\Phi =\begin{pmatrix} 1 & 0 & ... & 0\\ 0 & 1 & ... & 0\\ ... & ... & ... & ...\\ 0 & 0 & ... & 1<br /> \end{pmatrix}

\Psi =\frac{1}{\sqrt{N}}\begin{pmatrix}<br /> 1 & 1 & ... & 1\\<br /> 1 & W & ... ^{N-1}\\<br /> ... & ... ^{kn} &... \\<br /> 1 ^{N-1} & ...  ^{{N-1}^{2}}<br /> \end{pmatrix}

dimana: W={e}^{-i\frac{2\pi }{N}}

Nilai maksimum dari perkalian skalar vektor basis \left| \left< \phi,\psi \right> \right|  akan tercapai untuk sebarang pasangan vektor dan semuanya akan menghasilkan koherensi yang bernilai sama dengan 1/ \sqrt{N}  akibat faktor normalisasi dari basis Fourier. Dengan demikian GUP akan menjadi \left| {\Gamma }_{1} \right| + \left| {\Gamma }_{2} \right| \geq 2 \sqrt{N}  atau \left| supp(s) \right| + \left| supp(S) \right| \geq 2 \sqrt{N}  seperti pada WHUP untuk kasus diskrit yang telah diuraikan sebelumnya. Tingkat koherensi beberapa pasangan basis lain dapat dihitung dengan menggunakan kode Matlab berikut ini.

%Hitung koherensi beberapa pasangan basis
N=4; %dimensi basis
I=eye(N,N); %basis identitas / waktu
DFT=fft(eye(N,N))/(sqrt(N)); %basis DFT
H2=[1 1 1 1; 1 -1 1 -1; 1 1 -1 -1; 1 -1 -1 1]/2;
mu1=mu_GUP(I,DFT); %koherensi waktu-frekuensi
disp(sprintf(‘Koherensi basis waktu-frekuensi %f\n’,mu1));
mu2=mu_GUP(H2,DFT); %koherensi Hadamard-Fourier
disp(sprintf(‘Koherensi basis Hadamard-Fourier %f\n’,mu2));
% === FUNGSI %====
function [mu_value]= mu_GUP(PHI, PSI);
% fungsi untuk menghitung koherensi dua buah sistem basis
% untuk memperlihatkan Generalized Uncertainty Principle
% masukan: PHI, PSI dua buah matriks berdimensi sama
[M,N]=size(PHI);
G=PHI’*conj(PSI);
mu_value = sqrt(M)*max(max(abs(G)));

Seperti pada kompresi dengan matriks ortogonal, metoda CS (compressed sensing) juga mengambil komponen dominan dari suatu sistem basis. Namun demikian, basis yang dipakai bersifat overcomplete, dimana ada banyak sekali kemungkinan memilih basis yang sesuai untuk menyatakan suatu sinyal. Oleh karena itu perlu kendala tambahan agar solusi bisa ditemukan, yakni sinyal yang direkonstruksi dianggap bersifat sparse. Metoda CS memerlukan dua buah basis pada saat melakukan pencuplikan atau sensing, yaitu sparsity basis \Psi  dan projection basis \Phi  . Tingkat koherensi kedua basis menentukan batas minimum banyaknya cuplikan untuk merekonstruksi suatu sinyal secara eksak.

Posted in Compressive Sensing, Theory and Concept, Uncertainty Principles and Signal Compression | Leave a Comment »

Prinsip Ketidakpastian dan Kompresi Data (2)

Posted by suksmono on May 9, 2008

Prinsip ketidakpastian Weyl-Heisenberg memiliki konsekuensi yang penting untuk pengolahan sinyal, yaitu sinyal kontinyu yang tersebar dalam kawasan waktu akan terlokalisir dalam kawasan frekuensi dan demikian pula sebaliknya. Perubahan sinyal antar kawasan pada prinsip ini terhubung melalui transformasi Fourier. Tulisan ini memperlihatkan bahwa prinsip ketidakpastian dapat berlaku untuk sinyal waktu diskrit. Disini juga diperkenalkan adanya konsep kawasan yang lebih umum. Sifat penting transform Fourier berupa pengumpulan energi ke sejumlah kecil koefisien ternyata juga dimiliki oleh beberapa transformasi ortogonal atau uniter. Sifat ini dapat dimanfaatkan untuk melakukan kompresi sinyal secara sederhana, yaitu dengan memilih koefisien transform yang dominan dan membuang sisanya.

Prinsip Pengolahan Sinyal secara Dijital
Supaya dapat diolah oleh komputer atau pengolah DSP (Digital Signal Processor), sinyal-sinyal alami yang analog seperti suara atau gambar harus terlebih dahulu diubah menjadi sinyal dijital. Prinsip pengolahan sinyal secara dijital digambarkan dalam diagram blok berikut ini.

Diagram Blok Pengolahan Sinyal Dijital

Pengolah ini menerima masukan sinyal listrik analog dari pengindera atau transducer, yaitu sinyal s(t) yang merupakan berisikan campuran dari sinyal dikehendaki sa(t) dan sinyal lain s’a(t). Tapis (filter) antialiasing berfungsi untuk menghapus sa(t) sehingga tinggal sa(t) yang mengandung semua informasi didalam pita frekuensi yang akan diolah. Pada gambar diperlihatkan sinyal akan menjadi lebih ”halus” setelah melalui tapis ini karena komponen frekuensi tinggi diluar pita kerja sudah hilang.

Blok S/H (sampling and hold) mengambil cuplikan sinyal kontinyu dalam kawasan waktu secara periodik dan akibatnya sinyal berubah menjadi sinyal waktu diskrit sa(n), dimana n adalah bilangan bulat. Sinyal ini masih memiliki amplitudo kontinyu. Karena pengolah bekerja dengan register terbatas, sinyal ini terlebih dahulu diubah menjadi sinyal dijital oleh perangkat A/D (analog to digital converter). Sebuah kartu ADC biasanya telah dilengkapi dengan tapis antialiasing dan perangkat S/H. Keluaran dari A/D ini akan berupa sinyal dijital s(n).

Pengolah dijital (DSP processor) selanjutnya memproses sinyal sesuai dengan tujuan saat perancangan, misalnya penapisan lolos rendah, lolos tinggi, dst. Bagian ini mengubah sinyal asal s(n) menjadi sinyal y(n).

Supaya bisa disajikan kembali ke pengguna, sinyal dijital y(n) harus dikembalikan menjadi sinyal analog dengan memakai perangkat D/A (digital to analog converter). Selanjutnya sebuah tapis rekonstruksi untuk menghilangkan cacat akibat konversi dijital ke analog, dan akhirnya diperoleh keluaran ya(t) yang siap disajikan.

Aljabar Pengolahan Sinyal

Sinyal dijital dapat dianggap sebagai suatu deretan bilangan dan dapat dituliskan sebagai vektor kolom. Sebagai contoh, vektor s=[2 3 7]T menyatakan sebuah sinyal dijital sepanjang 3 cuplikan. Karena berbentuk vektor, manipulasi atau pengolahan sinyal ini dapat dilakukan dengan operasi-operasi matriks atau vektor, misalnya penjumlahan atau pengurangan, perkalian dengan skalar, perkalian dengan matriks, dll.

Salah satu proses yang sangat penting didalam DSP adalah transformasi sinyal. Pada dasarnya transformasi sinyal adalah proses mengubah sinyal dari suatu kawasan ke kawasan yang lain. Secara geometrik transformasi sinyal pada hakekatnya merupakan perubahan sistem koordinat.

Signal Transform

Pandangan Geometrik dari Transformasi Sinyal

Gambar diatas memperlihatkan sinyal asal s yang berada dalam sistem koordinat (x, y, z) diubah menjadi sinyal S dalam sistem koordinat baru (x’, y’, z’) melalui transformasi T. Kita bisa menganggap vektor basis dari koordinat asal sebagai {i, j, k}={(1,0,0), (0,1,0), (0,0,1)}, sedangkan vektor basis dari koordinat yang baru adalah {i’, j’, k’}.

Menurut Aljabar Linier, transformasi dari vektor s menjadi vektor S dapat dinyatakan sebagai perkalian antara vektor s dengan matriks transformasi T, yakni

S= Ts

Ada sekelompok transformasi penting didalam pengolahan sinyal yang disebut sebagai transformasi ortogonal, dimana semua vektor basis-nya memiliki panjang satu satuan dan bersifat saling tegak lurus. Transformasi ini bersifat mempertahankan panjang atau magnitudo vektor s. Secara geometrik transformasi yang demikian akan merupakan suatu perputaran sistem koordinat. Rotasi terhadap sumbu z sebesar sudut q dapat dinyatakan sebagai transformasi dengan matriks:

T=\begin{pmatrix} cos\left(\theta \right) & -sin\left(\theta \right) & 0\\ sin\left(\theta \right) & cos\left(\theta \right) &0 \\ 0 & 0 & 1 \end{pmatrix}

Akibat rotasi terhadap sumbu-z sebesar 45 derajat, maka sinyal s=[2 3 7]T pada contoh sebelumnya akan berubah menjadi sinyal S = [-0.71 3.54 7.00]T. Perhitungan ini dapat dilakukan dengan bantuan Matlab sebagai berikut:

%** Transformasi sinyal **
s=[2 3 7]‘, % sinyal asli
% ** Definisikan matriks transformasi **
T=[cos(pi/4) -sin(pi/4) 0; sin(pi/4) cos(pi/4) 0; 0 0 1];
S=T*s, % sinyal hasil transformasi

Matriks Transformasi Uniter dan Ortogonal

Transformasi uniter merupakan bentuk umum dari transformasi ortogonal dimana vektor basis-nya memiliki magnitudo satu satu satuan, saling tegak lurus, tetapi nilai komponennya berupa bilangan kompleks. Transformasi Fourier adalah contoh dari transformasi uniter. Matriks DFT dapat diperoleh secara mudah dengan Matlab dengan melakukan operasi fft terhadap suatu matriks identitas yang berukuran NxN

% cara-1: kalikan sinyal dengan matriks DFT
s=[2 3 7]‘, % sinyal asli
N=length(s);
U=fft(eye(N,N))/sqrt(N); %bentuk matriks DFT 3×3
S1= U*s, % hasil transformasi Fourier
%cara-2: langsung dengan perintah fft
S2=fft(s)/sqrt(N), %sama dengan S1

Tranformasi sinyal s pada contoh terdahulu memerlukan pembentukan matriks DFT berukuran 3´U. Perkalian dari matriks U dengan s menghasilkan sinyal kawasan frekuensi yaitu S, yang dalam kode Matlab diatas akan sama dengan S1 maupun S2.

U=\frac{1}{\sqrt{3}}\begin{pmatrix} 1 & 1 & 1\\ 1 & -0.5-0.9i & -0.5+0.9i \\ 1 & -0.5-0.9i & -0.5-0.9i \end{pmatrix}

S=U{s}^{T}=\begin{pmatrix}6.9\\ -1.7+2.0i\\ -1.7-2.0i \end{pmatrix}

Kita akan memakai hasil diatas untuk menunjukkan konsep-konsep penting yang dibahas sebelumnya. Yang pertama, suatu matriks U disebut uniter jika perkalian matriks ini dengan transpose dari konjugasi kompleksnya menghasilkan matriks satuan, atau U*T×U = I. Kedua, kita akan memeriksa berlakunya Teorema Parseval yang menyatakan bahwa pengukuran energi dikawasan waktu akan sama dengan hasil pengukuran energi dikawasan frekuensi, atau s*T×s = S*T×S. Dan yang terakhir, kita akan melihat sifat pengkompakan energi dengan memeriksa sebaran energi dikawasan frekuensi E = S×S*. Script Matlab berikut dapat ditambahkan pada script sebelumnya:

S=S1; % dan ini sama juga dng S=S2
% Cek sifat uniter, T. Parseval, dan sebaran
conj(transpose(U))*U, %jika uniter akan menghasilkan I berukuran NxN
E_waktu=sum(conj(transpose(s))*s), %Energi kawasan waktu
E_frekuensi=sum(conj(transpose(S))*S) , %Energi kawasan frekuensi
E_dist_frek = S.*conj(S), %sebaran energi

Analisis terhadap hasil perhitungan pada kode diatas memberikan kesimpulan bahwa U memang memiliki sifat uniter karena hasil perkalian dengan konjugasi transpose-nya menghasilkan matriks identitas. Perhitungan juga memberikan nilai energi kawasan waktu sebesar E_waktu=62 dan energi kawasan frekuensi sebesar E_frekuensi=62, jadi keduanya identik dan Teorema Parseval telah dipenuhi. Hal yang juga menarik terdapat pada sebaran energi dikawasan frekuensi, yakni E_dist_frek=[48.0 7.0 7.0]T. Disini terlihat bahwa energi hasil transformasi telah terkumpul pada sejumlah kecil koefisien, dalam hal ini hanya pada koefisien yang pertama saja.

Disamping matriks DFT yang dijelaskan diatas, beberapa matriks ortogonal atau uniter lain juga memiliki sifat pengkompakan energi. Disini akan ditinjau dua buah matriks ortogonal yang sering dipakai untuk melakukan transformasi sinyal, yaitu matriks DCT dan matriks Hadamard. Terlebih dahulu sinyal asal diperpanjang menjadi empat cuplikan, misalnya s=[2 3 7 11]T. Dengan demikian haruslah dipilih matriks DCT dan Hadamard yang ukurannya 4´4, masing-masing diberi notasi D dan H2.

{H}_{2}=\frac{1}{2} \begin{pmatrix} 1 & 1 & 1 & 1\\ 1 & -1 & 1 & -1 \\ 1 & 1 & -1 & -1\\ 1 & -1 & -1 & 1 \end{pmatrix}

Matriks DCT dapat diperoleh dengan perintah Matlab dct(eye(4,4)), sedangkan matriks Hadamard bisa dibentuk secara rekursif dengan konstruksi Sylvester, yaitu melalui perkalian Kronecker dengan matriks dasar H1.

{H}_{1}=\frac{1}{\sqrt{2}}\begin{pmatrix}1 & 1\\ 1 & -1 \end{pmatrix}
{H}_{n}={H}_{n-1} \otimes {H}_{1}

Script Matlab berikut ini dapat dipakai untuk mempelajari sifat-sifat penting kedua transformasi ortogonal tersebut.

%transformasi ortogonal: DCT dan Hadamard
s=[2 3 7 11]‘, % sinyal asli
N=length(s);
%bentuk matriks transformasi
D=dct(eye(N,N)); %bentuk matriks DCT 4×4
H2=[1 1 1 1; 1 -1 1 -1; 1 1 -1 -1; 1 -1 -1 1]/2;
S_dct= D*s, % hasil transformasi DCT
S_had= H2*s, % hasil transformasi Hadamard
E_asli= sum(transpose(s)*s), %energi sinyal asal
E_dct =sum(transpose(S_dct)*S_dct) , %Energi kawasan DCT
E_hadamard =sum(transpose(S_had)*S_had) ,
%Energi kawasan Hadamard
E_dist_dct = S_dct.*S_dct, %sebaran energi DCT
E_dist_had = S_had.*S_had, %sebaran energi Hadamard
%gambarkan sebaran energi
figure(1); bar(E_dist_dct); title(‘Sebaran energi di kawasan DCT’);
figure(2); bar(E_dist_had); title(‘Sebaran energi di kawasan Hadamard’);

Perhitungan dengan script diatas menegaskan kembali berlakuknya hukum konservasi energi yang tersirat dari Teorema Parseval, karena diperoleh hasil E_asli=E_dct=E_had = 183. Sifat pengkompakan energi juga ditemukan dikedua transformasi ini. Sebaran masing-masing energi digambarkan dalam histogram berikut ini

Energi DCT

Sebaran Energi Sinyal di Kawasan DCT

Energi Hadamard

Sebaran Energi Sinyal di Kawasan Hadamard

Kompresi Data dengan Transformasi Ortogonal

Kompresi data sederhana dapat dilakukan dengan cara memilih beberapa koefisien yang memiliki energi dominan dan membuang yang lain. Sebagai contoh, hasil transformasi DCT untuk sinyal s adalah S=[11.5 -6.9 1.5 0.1]T. Jika dipilih dua komponen dominan saja (untuk kompresi sebesar dua kali) dan mengisi komponen sisanya dengan nol, maka kita akan memperoleh pendekatan dari S , yaitu S_hat = [11.5 -6.9 0 0]. Hasil transformasi balik sinyal ini (dengan perintah Matlab: D’*[11.5 -6.9 0 0]‘ ) akan memberikan nilai pendekatan dari sinyal asal sebesar s_hat=[1.2 3.9 7.6 10.3]. Sebaliknya, jika dua komponen tak-dominan yang dipilih, maka dengan tingkat kompresi yang sama akan diperoleh s_hat1= D’*[0 0 1.5 0.1]‘=[0.8 -0.8 -0.7 0.7]T yang sama sekali berbeda dengan sinyal aslinya. Ukuran dari baik-buruk-nya pendekatan sinyal yang satu terhadap yang lain dapat ditentukan dengan menghitung jarak Euclidian kedua sinyal. Caranya: kurangkan kedua sinyal komponen demi komponen, lalu kuadratkan setiap beda komponen dan jumlahkan. Akar kuadrat dari hasil terakhir menunjukkan jarak yang kita cari. Analisis untuk transformasi Hadamard sejalan dengan penjelasan ini, tinggal menggantikan D dengan H2.

Pada tulisan sebelumnya telah diperoleh hasil berupa implikasi tak-langsung dari WHUP, yaitu sinyal yang tersebar dikawasan waktu akan terkumpul di kawasan frekuensi dan begitu pula sebaliknya. Tulisan ini menunjukkan hasil WHUP dengan Matlab dengan cara menghitung energi dari koefisien transformasi Fourier diskrit. Sifat ini juga dimiliki oleh transformasi ortogonal atau uniter yang lainnya, misalnya transformasi DCT dan transformasi Hadamard. Sebuah teknik kompresi data sederhana, yaitu dengan mengambil dua dari komponen dominan koefisien transformasi DCT, telah dijelaskan. Hasil inversi dari data kawasan transformasi yang telah dikompresi dengan cara ini lebih menyerupai sinyal aslinya jika dibandingkan dengan pilihan koefisien lainnya.

Kompresi JPEG (versi lama) memakai DCT dan bukan transformasi Fourier karena pertimbangan praktis. Menyimpan koefisien kompleks yang terdiri dari bagian riil dan imajiner terlalu merepotkan dan kurang efisien. Lagipula, DCT adalah hampiran yang cukup baik dari transformasi ideal KLT (Karhunen-Loeve Tranform) untuk sinyal acak Gauss-Markov orde satu, sebuah model yang mendekati sinyal alami. Cara kerja kompresi JPEG mirip dengan kompresi sederhana yang telah dijelaskan dengan tambahan adanya kompresi tak-merugi berupa entropy coding.

Penutup dan Tulisan Mendatang

Pada “kompresi” DCT diatas telah dipilih komponen yang dominan dan letaknya bisa berubah bergantung pada sinyal masukannya, meskipun pada JPEG letak komponen dan alokasi bit-nya telah ditentukan sebelumnya. Pemilihan yang mengikuti bentuk sinyal seperti ini disebut sebagai pemilihan koefisien secara adaptif; ini adalah masalah yang sebaiknya dihindari. Terlebih lagi kompresi dilakukan pada data dijital yang pada saat dikumpulkan (pemotretan gambar atau perekaman suara) menghasilkan sejumlah besar data yang kemudian dibuang pada saat kompresi. Teknik terbaru yang disebut sebagai Compressed Sensing/ Compressive Sampling (CS) menghindari cara yang tidak efisien seperti ini. Disamping itu, CS juga tidak memerlukan pencuplikan koefisien dominan secara adaptif karena semua komponen sama pentingya, dapat diambil yang manapun asalkan batas minimum jumlah cuplikan terpenuhi. Jumlah cuplikan yang diperlukan untuk rekonstruksi sinyal didalam CS ternyata jauh dibawah batas Nyquist yang dinyatakan oleh Teorema Shannon. Teori CS memerlukan generalisasi dari hubungan waktu-frekuensi supaya prinsip ketidakpastian tetap berlaku. Tulisan mendatang akan mengulas teknik CS ini dengan lebih rinci.

Posted in Compressive Sensing, Theory and Concept, Uncertainty Principles and Signal Compression | Leave a Comment »

Prinsip Ketidakpastian dan Kompresi Data (1)

Posted by suksmono on May 2, 2008

Selama ini teknik kompresi data (sinyal) telah banyak dipakai untuk mengefisienkan penyimpanan dan pengiriman suara atau gambar pada peralatan komputer maupun ponsel. Biasanya orang mengenal kompresi sebagai metoda yang berkembang dalam bidang pengolahan sinyal dijital. Tanpa disadari, sebenarnya ada kaitan erat antara paradigma kompresi sinyal ini dengan prinsip fundamental yang dikenal dalam bidang Fisika Kuantum yaitu prinsip ketidakpastian.

Siapapun yang pernah belajar Fisika Modern tentu mengenal prinsip ketidakpastian Hiesenberg atau HUP (Heisenberg Uncertainty Principle). Prinsip ini menyatakan bahwa pengukuran posisi dan momentum suatu partikel tidak mungkin dibuat teliti kedua-duanya sekaligus, yakni

\Delta p\Delta x\geq \frac{\hbar}{2 }

HUP adalah konsekuensi langsung dari sifat dualitas gelombang-partikel. Sebelum Teori Kuantum muncul, orang mengira bahwa ketelitian pengukuran posisi suatu partikel dan momentum (pada saat bersamaan) hanya akan dibatasi oleh ketelitian alat ukur. Prinsip diatas menyatakan larangan alam untuk mengetahui keduanya sekaligus secara teliti, sepresisi apapun alat ukur yang digunakan.

Selain relasi ketidakpastian dari momentum dengan posisi, HUP juga dapat dinyatakan sebagai hubungan ketidakpastian antara energi dengan waktu, yaitu

\Delta E\Delta t\geq \frac{\hbar}{2 }

Menurut Max Planck, suatu foton yang memiliki frekuensi \omega, akan memiliki energi sebesar \hbar\omega dimana \hbar adalah konstanta Planck. Dengan demikian HUP untuk energi-waktu akan setara dengan

\Delta \omega \Delta t\geq \frac{1}{2}

Pertidaksamaan diatas dinamakan juga prinsip ketidakpastian Weyl-Heisenberg (WHUP) yang berbunyi:

a continuous-time signal cannot be simultaneously well-localized in both time and frequency

Apakah bedanya dengan prinsip ketidakpastian yang sebelumnya? Jika HUP berbicara mengenai partikel dan gelombang sebagai objek fisik, WHUP menyatakan sifat umum dari prinsip ketidakpastian akan berlaku untuk sebarang sinyal atau fungsi kontinyu. Dalam bahasa pengolahan sinyal bisa ditafsirkan bahwa jika suatu sinyal s(t) terlokalisir dalam kawasan waktu maka transform Fourier dari sinyal ini,

F\left(s\left(t \right) \right)=S\left(\omega \right)

akan tersebar dikawasan frekuensi dan demikian pula sebaliknya, suatu sinyal yang terlokalisir dikawasan frekuensi akan memiliki transform Fourier yang tersebar dikawasan waktu.

Dirac Delta Function
Gambar 1. Pasangan transform Fourier dari sinyal delta Dirac

Contoh dari sinyal kontinyu terlokalisir waktu adalah sinyal delta Dirac. Transform Fourier dari sinyal ini akan tersebar keseluruh rentang frekuensi, seperti yang diperlihatkan oleh Gambar 1. Karena delta Dirac adalah fungsi riil, maka haruslah transform Fouriernya bernilai kompleks. Gambar 1 sebelah kanan hanya menunjukkan nilai magnitude dari koefisien Fourier-nya untuk memperlihatkan sebaran komponen dominan dan terlihat disini bahwa semua komponen sama dominan-nya.

Fungsi delta dapat dilihat sebagai fungsi (distribusi) Gaussian dengan limit variansi mendekati nol. Kita akan memakai fungsi Gaussian ini untuk memperlihatkan (bukan membuktikan) bahwa implikasi tidak langsung dari WHUP adalah,

sinyal kontinyu yang tersebar dalam kawasan waktu, akan terlokalisir pada kawasan frekuensi dan begitu pula sebaliknya.

Pasangan Fourier dari fungsi Gaussian dengan mean nol dapat dituliskan sebagai berikut

g\left(t \right)=\frac{1}{\sigma \sqrt{2\pi }}{e}^{-\frac{t^2}{\left( \sigma\right)^2}}\leftrightarrow G\left(\omega \right)={e}^{-\frac{\omega ^2}{\left(1/\sigma\right)^2}}

Dari relasi diatas terlihat bahwa peningkatan variansi di kawasan waktu atau pelebaran sinyal Gaussian akan mengakibatkan penyempitan spektrum dikawasan frekuensi dan demikian pula sebaliknya, penyempitan sinyal di kawasan waktu akan melebarkan spektrum-nya. Dengan demikian, implikasi tak langsung dari WHUP tersebut menjadi jelas, khususnya untuk fungsi Gaussian. Fungsi-fungsi atau sinyal-sinyal waktu-kontinyu yang lain juga akan memiliki sifat yang demikian. Matlab Script berikut ini dapat dipakai untuk memperlihatkan fenomena ini.

%—————————————————–———-
% Script Matlab untuk memperlihatkan WHUP
%—————————————————————–
clear;clc;
t=-5:0.1:5; w=t; %koordinat waktu dan frekuensi.
sigma1=1; sigma2=0.5; sigma3=0.25;

%Hitung Fungsi Gaussian
g1=normal_baku(sigma1,t); g2=normal_baku(sigma2,t); g3=normal_baku(sigma3,t);

%Hitung transform Fourier-nya
G1=F_normal_baku(sigma1,w);
G2=F_normal_baku(sigma2,w);
G3=F_normal_baku(sigma3,w);

figure(1);plot(t,g1,t,g2,t,g3);
legend(’sigma’,’0.5sigma’,’0.25sigma’)
figure(2);plot(w,G1,w,G2,w,G3);
legend(’sigma’,’0.5sigma’,’0.25sigma’)

%—————————————————–
% Fungsi Gaussian
%—————————————————–
function [g]= normal_baku(sigma,t)
g=(1/sqrt(2*pi*sigma*sigma))*exp(-t.*t/(sigma*sigma));

%————————————————————————–
% Transform Fourier dari fungsi Gaussian
%————————————————————————–
function [G]= F_normal_baku(sigma,w)

G=exp(-w.*w*(sigma*sigma));

Keluaran dari program Matlab diatas diperlihatkan pada Gambar 2 dan Gambar 3. Disini jelas terlihat penyempitan fungsi Gauss pada kawasan waktu akan mengakibatkan pelebaran pada kawasan frekuensi.


Gambar 2. Fungsi Gaussian dengan berbagai nilai variansi
Fourier Transform of Gaussian
Gambar 3. Transform Fourier dari fungsi Gaussian pada Gambar 2

Sebelum diskusi dilanjutkan, akan terlebih dahulu disinggung sekilas mengenai suatu teorema yang sangat penting dalam transformasi sinyal, yaitu Teorema Parseval. Teorema ini pada prinsipnya adalah Hukum Kekekalan Energi dari sinyal, yang menyatakan bahwa pengukuran energi dalam kawasan waktu (misalnya dengan bantuan Osciloscope) dan pengukuran energi pada kawasan frekuensi (misalnya dengan Spectrum Analyzer) akan memberikan hasil yang sama. Teorema ini juga berimplikasi bahwa penghilangan komponen frekuensi tertentu dari sinyal akan mendistorsi sinyal kawasan waktunya, sepadan dengan magnitudo koefisien Fourier tersebut. Dengan demikian, hilangnya komponen dengan magnitudo rendah tidak akan mengubah bentuk sinyal semula terlalu banyak.

Untuk aplikasi pengolahan sinyal atau citra dijital, perhitungan dilakukan terhadap data dijital yang tak lain adalah fungsi-fungsi diskrit. Fungsi ini diperoleh dari fungsi kontinyu dengan melalui proses pencuplikan dan kuantisasi yang dapat dilakukan dengan sebuah ADC (Analog to Digital Converter).

  • Apakah WHUP juga berlaku untuk fungsi diskrit?
  • Apakah hubungan antara prinsip ketidakpastian, Teorema Parseval, dan kompresi data?
  • Selain transformasi Fourier, adakah transformasi lain yang punya sifat mirip?
  • Mengapa justru DCT dan bukannya DFT yang dipilih untuk kompresi JPEG?

Jawaban dan penjelasan lebih lanjut dari pertanyaan-pertanyaan diatas akan disajikan pada tulisan mendatang di blog ini.

Posted in Compressive Sensing, Image Processing, Theory and Concept, Uncertainty Principles and Signal Compression | 5 Comments »