Radar Kompresif Dipamerkan pada Seminar Radar Nasional 2009

Hari Kamis tanggal 30 April 2009, LIPI bekerjasama dengan ITB menyelenggarakan Seminar Radar Seminar III di Hotel Savoy Homann, Bandung. Pada kesempatan ini prototip radar SFCW-GPR (700-1400 MHz) dengan kemampuan kompresif dipamerkan. Sebelumnya, telah dilakukan ujicoba scanning radar di LTRGM-STEI ITB, baik dengan loopback maupun antena-ke-antena. Menurut Prof. Leo P. Ligthart, yang saat ini sedang melakukan kunjungan ke Indonesia, radar kompresif berpotensi membuat terobosan baru di bidang radar, khususnya radar jenis SFCW.

Pak Koes dan Tommi menjaga stand pameran

Prototip Radar Kompresif 700-1400 MHz

Radar penembus permukaan (GPR/Ground Penetrating Radar) memiliki berbagai aplikasi sipil dan militer yang luas. Dalam bidang sipil, radar ini dipakai pada eksplorasi geofisika untuk mencari batubara atau sumber air;  dalam arkeologi untuk mencari situs purbakala yang terbenam dibawah tanah; dan juga untuk mendeteksi keretakan beton pada bangunan, kerusakan jalan raya atau jembatan, serta menentukan lokasi pipa,  kabel, atau utulitas di bawah tanah. Pada bidang militer, radar penembus permukaan berguna untuk mencari ranjau, tempat penimbunan senjata, ruang dibawah tanah, maupun untuk mendeteksi musuh dibalik dinding.

Sejak awal tahun 2000-an, IRCR-TU Delft berusaha meningkatkan kecepatan dari SFCW-GPR, agar dapat bersaing dengan radar impuls. SFCW memiliki banyak keunggulan dalam berbagai hal dibanding sistem impuls, kecuali satu kelemahan dalam kecepatan akuisisi data. Peneliti di IRCTR berusaha mengatasi hal ini, antara lain dengan membuat radar 8 kanal yang bekerja secara serempak. Inovasi baru melalui penginderaan kompresif (Compressive Sensing) memberi harapan baru berhasilnya program ini.

Pencuplikan kompresif jika rapat spektral daya sinyal diketahui

Salah satu jalur penelitian dalam CS (Compressive Sensing/Compressive Sensing) berhubungan dengan upaya peningkatan kinerja; baik  berupa hasil rekonstruksi yang lebih baik dan/ atau jumlah cuplikan yang lebih sedikit. Ini dilakukan, antara lain, dengan memanfaatkan pengetahuan tentang sifat dari sinyal atau prior.

Bagi seorang praktisi CS,  salah satu pertanyaan yang lebih penting untuk dijawab adalah:  “Jika rapat spektral daya dari sebuah sinyal diketahui dan jumlah cuplikan maksimum yang diperbolehkan dibatasi, bagaimana cara memilih lokasi titik cuplikan yang memberikan kinerja obyektif CS yang terbaik?”

Berdasarkan hasil penelitian, jika lokasi cuplikan mengikuti sebaran energi, ternyata kinerja rekonstruksi sinyal dapat ditingkatkan—atau, dengan persyaratan kinerja rekonstruksi yang diberikan, jumlah cuplikan yang diperlukan menjadi jauh lebih sedikit. Hasil ini dapat dipakai untuk mempercepat proses pengumpulan data pada radar SFCW.

Compressive SFCW-GPR/ Pencitraan Kompresif (3)

Radar pensintesa frekuensi SFCW-GPR memiliki kelemahan dalam hal kecepatan pengumpulan data. Inilah salah satu alasan mengapa sistem radar ini kurang berkembang. Berbagai upaya telah dilakukan untuk mengatasi masalah ini, misalnya dengan mengirimkan beberapa frekuensi sekaligus secara serempak. Munculnya teknik pencitraan kompresif memungkinkan masalah ini dapat diatasi dengan pendekatan lain. Tulisan ini menjelaskan prinsip kerja SFCW-GPR kompresif dan hasil-hasil penelitian yang telah dieproleh hingga saat ini.

Cerita lengkap dapat dibaca dari presentasi di sini atau di situ.

In English

Data acquisition speed is an inherent problem of the stepped-frequency continuous wave (SFCW) radars, which discouraging further usage and development of this technology. We propose an emerging paradigm called the compressed sensing (CS), to manage this problem. In the CS, a signal can be reconstructed exactly based on only a few samples below the Nyquist rate. Accordingly, the data acquisition speed can be increased significantly. A novel design of SFCW ground penetrating radar (GPR) with a high acquisition speed capability is proposed and evaluated. Simulation by a monocycle waveform and actual measurement by a Vector Network Analyzer in a GPR test-range confirm the implementability of the proposed system.

You can read the full story here or there.

Pencitraan Kompresif (1.b)

Pencitraan kompresif pada kamera piksel-tunggal dapat dilakukan dalam kawasan spasio-temporal maupun kawasan frekuensi. Pencitraan kawasan spasio-temporal merupakan tafsiran langsung dari pencuplikan kompresif; lebih sederhana, dan intuitif. Namun demikian, penggunaan DCT sebagai transformasi penjarang tidak memberikan hasil yang memuaskan. Tulisan ini memberikan alternatif transformasi penjarang berbentuk wavelet paling sederhana, yaitu transformasi Haar.

Transformasi Haar untuk Pencuplikan Kompresif

Ketidaksesuaian DCT sebagai transformasi penjarang terlihat dari peluruhan koefisien dari citra berbentuk ”H” yang relatif kurang cepat. Hal ini telah dapat diperbaiki dengan meminimumkan nilai TV (total variance). Cara ini cukup ampuh, namun kurang intuitif karena tidak secara langsung memberikan transformasi pada kawasan basis penjarang-nya. Tulisan ini menjelaskan penggunaan transformasi Haar, salah satu bentuk transformasi wavelet paling tua dan paling sederhana, sebagai transformasi penjarang-nya.

Konstruksi matriks penjarang PSI pada kode sebelumnya diubah dari transformasi DCT menjadi transformasi Haar dengan mengubah baris terkait sbb :

PSI = haar_T(N)’;%(eye(N,N)); % ** buat matriks transformasi Haar **

dimana ”haar_T” diatas adalah fungsi untuk membentuk matriks transformasi Haar yang kode Matlab-nya diperlihatkan di akhir tulisan.

Hasil dari transformasi Haar dengan citra dan parameter kompresi sebesar 1.4 (seperti pada tulisan pertama) diberikan pada Gambar 1 berikut ini.

(a)

(b)

Gambar 1. Hasil rekonstruksi pencitraan kompresif dengan transformasi Haar: (a) K-buah koefisien pertama, (b) 50 buah koefisien pertama.

Gambar 2. Plot koefisien terurut dari citra ”H” terhadap basis Haar

Hasil rekonstruksi terlihat masih kurang baik, meskipun bagian utama dari citra sudah terambil. Ini bisa dijelaskan dengan melihat kembali plot koefisien penyusun citra hasil optimasi yang ditampilkan pada Gambar 2. Terlihat bahwa meskipun koefisien meluruh jauh lebih cepat dari DCT, nilai koefisien nol tercapai pada indeks urutan ke sekitar 50. Dengan demikian, asumsi derajat sparsity sebesar 32 masih kurang mencukupi. Bahkan, melibatkan 50 buah koefisien hasil optimisasi-pun masih belum menghasilkan rekonstruksi yang sempurna, seperti diperlihatkan pada Gambar 1.b.

Gambar 3. Citra asli papan catur

Sebagai perbandingan, berikut ini disajikan hasil pencitraan kompresif dengan basis Haar untuk citra papan catur (checker board) yang diperlihatkan pada Gambar 3. Relatif terhadap basis Haar, citra ini sangat sparse, seperti diperlihatkan pada Gambar 4.

Gambar 4. Koefisien terurut citra papan-catur hasil pemilihan basis

Kali ini tingkat kompresi diperbesar menjadi 2.3 kalinya, dengan cara memperkecil K menjadi sekitar . Gambar 5 menunjukkan citra teramati, dimana piksel abu-abu menyatakan piksel yang tak terlihat oleh pengamatan DMD. Hasil rekonstruksi diperlihatkan pada Gambar 5. Terlihat bahwa rekonstruksi eksak sudah tercapai.

(a)

(b)

Gambar 5. (a) Citra teramati dan (b) citra rekonstruksi

Kode Matlab

% ——————————————————————

% ** fungsi membentuk matriks transformasi Haar **

% ——————————————————————

function [H]=haar_T(h_size);

n=ceil(log2(h_size));N=2^n;

H=zeros(N,N);

norm_factor=1/sqrt(N);

%difine z-axis

z=(0:1:N-1)/N;

H(1,:)=1;

for p=0:n-1;

for q=1:2^p;

k=2^p+q-1;

%disp(sprintf(‘k=%d , p=%d, q= %d \n’,k,p,q)); % debug p and q

hk=zeros(N,1);

hk_1=2^(p/2); hk_2=-(2^(p/2));

lb_1=ceil(N*(q-1)/2^p);ub_1=ceil(N*(q-0.5)/2^p);

lb_2=ceil(N*(q-0.5)/2^p);ub_2=ceil(N*q/2^p);

% *debug z-bound

% disp(sprintf(‘lb_1=%d , ub_1=%d, lb_2= %d, ub_2= %d\n’,lb_1,ub_1,lb_2,ub_2));

H(k+1,lb_1+1:ub_1)=hk_1;

H(k+1,lb_2+1:ub_2)=hk_2;

end

end;

H=H*norm_factor;

Pencitraan Kompresif (2)

Pemilihan transformasi penjarang yang tepat sangat berpengaruh pada hasil rekonstruksi pencitraan kompresif. Tulisan ini menjelaskan pemakaian teknik minimisasi variansi total (TV-total variance), salah satu estimasi kokoh (robust) yang sebenarnya sudah sering dipakai namun baru disadari hubungannya dengan pencuplikan kompresif. Karena parameter dari TV minimum saja tidak bisa dipakai untuk melakukan inversi, maka harus ada basis lain yang terhubung dengan seleksi TV minimum. Diagram ranah-jamak dapat dipakai untuk menjelaskan keterkaitan antar transformasi dan parameter yang digunakan dalam pencitraan kompresif kawasan frekuensi.

Pencuplikan Kompresif Kawasan Frekuensi

Operasi pengamatan pada kawasan spasial yang diberikan pada persamaan (1) didalam tulisan sebelumnya dapat ditafsirkan sebagai pencuplikan terhadap vektor basis penjarang jika dituliskan sebagai

(1)

Setiap baris pada hanya berisi satu buah angka bernilai satu, sedangkan sisanya nol. Perkalian dengan mengakibatkan seleksi vektor baris pada indeks yang bersangkutan dengan letak kolom elemen aktif tersebut. Sebagai akibatnya, vektor basis dalam akan berkurang dan letaknya diacak.

Kinerja rendah yang dihasilkan akibat penggunaan DCT sebagai transformasi penjarang dapat ditingkatkan dengan transformasi bentuk lain. Salah satu yang telah sering dipakai tetapi baru disadari kaitannya dengan pencuplikan kompresif adalah variansi total (TV-total variance). Pendekatan dengan cara ini cukup berbeda jika dibandingkan dengan metoda sebelumnya yang merupakan tafsiran langsung pencuplikan spasial. Kaitan berbagai transformasi dan proses recovery sinyal dilukiskan dalam diagram ranah jamak berikut ini.

Gambar 1. Diagram ranah-jamak [1] pencuplikan kompresif kawasan frekuensi

Mengacu pada diagram ranah-jamak, sinyal s terlebih dahulu ditransformasikan ke kawasan frekuensi menjadi . Pencuplikan kompresif kawasan ini dengan transformasi pengukur akan menghasilkan cuplikan kompresif . Karena basis penjarang ada dalam kawasan lain yang perhitungannya didasarkan pada sinyal kawasan spasio-temporal, maka hasil cuplikan kompresif terlebih dahulu dikembalikan ke kawasan semula menjadi sinyal . Selanjutnya yang akan diminimumkan adalah nilai variansi total:

(2)

Basis pursuit untuk kasus ini dapat dirumuskan sebagai

(P1): min |S| s.t. (3)

Optimisasi diatas mengatakan:

“pilih yang inverse transform Fouriernya, yaitu memiliki variansi total |S| terkecil dan bisa menjelaskan hasil pengamatan .

Pada diagram, pengaruh hasil pemilihan |S| terkecil pada koefisien Fourier digambarkan sebagai panah bergaris ganda pada bagian kiri diagram. Operasi pengambilan nilai variansi total digambarkan sebagai transformasi merugi karena pengetahuan nilai variansi total saja tidak dapat dipakai untuk merekonstruksi sinyal semula.

Eksperimen dan Analisis

Pada eksperimen dipakai citra yang sama seperti sebelumnya, tetapi mekanisme rekonstruksi-nya sangat berbeda; yang dapat disimpulkan dengan membandingkan diagram ranah jamak pada Gambar 1 dengan yang ada pada tulisan sebelumnya. Fungsi optimisasi yang digunakan ada dalam paket program l₁-magic dari Caltech. Program utama sudah sedikit dimodifikasi supaya sesuai dengan keperluan, seperti diperlihatkan pada lampiran di akhir tulisan ini.

Gambar 2. Citra asli

Gambar 3. Citra teramati

Gambar 2 menunjukkan citra asli, sedangkan Gambar 3 adalah citra teramati dilihat dari kawasan spasial. Piksel berharna abu-abu menyatakan ketiadaan data pada posisi bersangkutan, sedangkan piksel hitam dan putih menyatakan hasil pengamatan terkompresi sebanyak 4 kali. Tidak seperti pada tulisan sebelumnya, sangat sulit menebak citra sebenarnya dari citra teramati ini karena tingkat kompresi yang lebih tinggi.

Gambar 4. Citra hasil rekonstruksi berdasarkan TV

Hasil rekonstruksi yang diperlihatkan pada Gambar 4 jauh lebih bagus dari yang sebelumnya—bahkan tepat sama dengan citra asalnya. Disini terlihat betapa pentingnya pemilihan basis penjarang yang tepat didalam pencuplikan kompresif. Nilai TV menggambarkan smoothness karena berhubungan dengan variansi atau keragaman citra. Meminimumkan nilai ini berarti memilih citra dengan variansi rendah, atau smoothness yang tinggi. Ini sangat cocok untuk citra-citra tomografi dalam biomedika, karena organ tubuh manusia sudah semestinya smooth.

Daftar Pustaka

1. AB. Suksmono, ”A graphical representation of multi-domain signal processing,” Proc. of SICE-ICASE 2006, Busan, Korea.

2. EJ. Candes, J. Romberg, and T. Tao, ”Robust uncertainty principles: Exact signal recovery from highly incomplete frequency information,” IEEE Trans. Information Theory, Vol.52, no.2, Feb.2006, pp.489-509.

3. MF. Duarte, MA. Davenport, D. Takhar, JN. Laska, T. Sun, KF. Kelly, and RG. Baraniuk, “Single-pixel imaging via compressive sampling,” IEEE Signal Proc. Magazine (83), March 2008, pp.83-90.

4. AB. Suksmono, “A brief review on compressive imaging,” Doc. RGQ14-2/2/ -E/Jul 2008, ITU-D, Rapp. Group Meeting on Q.14-2/2, Tokyo, 3-4 July 2008.

Kode Matlab

% *******************************************

% kamera piksel tunggal: cs_camera.m

% ditulis kembali berdasarkan kode L1-magic

% fungsi-fungsi lain dapat di-download dari Caltech

% *******************************************

% use implicit, matrix-free algorithms ?

largescale = 1;

path(path, ‘./Optimization’); path(path, ‘./Measurements’); path(path, ‘./Data’);

% read original image

I=double(imread(‘IMG\H.bmp’,'bmp’));

%tampilkan citra asal;

figure(1); imagesc(I); colormap(gray); title(‘citra asli’);

[n_row,n_col] = size(I); n = n_row; % panjang baris/kolom

N = n_row*n_col; % total cuplikan asal/ jumlah piksel

I = I/norm(I(:)); %normalisasi energi citra

I = I – mean(I(:)); %mean dibuat nol

x = reshape(I,N,1); %ubah menjadi larik 1-D /vektor

% jumlah observasi yang diinginkan

f_comp=4;%faktor kompresi

K = N/f_comp;%jumlah pengamatan

% Simulasi DMD kawasan frekuensi

P = randperm(N)’;

q = randperm(N/2-1)+1;

OMEGA = q(1:K/2)’;

% measurement matrix

if (largescale)

A = @(z) A_f(z, OMEGA, P);

At = @(z) At_f(z, N, OMEGA, P);

% obsevations

b = A(x);

% initial point

x0 = At(b);

else

FT = 1/sqrt(N)*fft(eye(N));

A = sqrt(2)*[real(FT(OMEGA,:)); imag(FT(OMEGA,:))];

A = [1/sqrt(N)*ones(1,N); A];

At = [];

% observations

b = A*x;

% initial point

x0 = A’*b;

end

I0=reshape(x0, n, n);

figure(2);imagesc(I0);colormap(gray);title(‘citra teramati’);

tvI = sum(sum(sqrt([diff(I,1,2) zeros(n,1)].^2 + [diff(I,1,1); zeros(1,n)].^2 )));

disp(sprintf(‘Original TV = %.3f’, tvI)); %TV asal

time0 = clock;

xp = tveq_logbarrier(x0, A, At, b, 1e-3, 5, 1e-8, 200);

Ip = reshape(xp, n, n);

disp(sprintf(‘Total elapsed time = %f secs\n’, etime(clock,time0)));

figure(3);imagesc(Ip);colormap(gray);title(‘citra rekonstruksi’);

Pencitraan Kompresif (1)

Sebuah kamera dijital terdiri dari suatu susunan optik dan sebuah CCD (Charged-Coupled Device) yang berfungsi untuk menangkap bayangan benda, lalu mengubahnya menjadi citra dijital. Jika tidak ada CCD, satu buah fotodioda-pun dapat dipakai untuk membuat kamera dijital. Bagaimana caranya?

Kamera Piksel-Tunggal

Technology Review adalah sebuah majalah bulanan yang sering dijadikan acuan untuk melihat perkembangan teknologi terkini. Seperti pada tahun-tahun sebelumnya, sekitar bulan Maret majalah ini menampilkan sepuluh buah teknologi terkini yang akan menjadi trend di masa datang. Pada tahun 2006, satu dari 10 teknologi ini menarik perhatian khalayak luas karena dianggap sangat tidak biasa, yaitu Kamera Piksel-Tunggal. Kamera ini tidak memakai CCD sebagai pengubah bayangan benda menjadi gambar dijital, melainkan hanya satu buah fotodioda saja.

Sebuah kamera piksel-tunggal, seperti yang diperlihatkan pada Gambar 1, memanfaatkan susunan cermin renik (DMD-Digital Micromirror Device) yang dapat dikendalikan komputer untuk melakukan penyapuan bayangan benda secara acak. Bagian bayangan yang jatuh pada elemen cermin yang sedang ”aktif” diteruskan kesebuah fotodioda untuk diubah menjadi sinyal listrik. Perangkat DMD biasanya dipakai pada televisi proyeksi.

Gambar 1. Kamera piksel-tunggal

Informasi elemen DMD aktif dan intensitas cahaya yang sudah didijitasi selanjutnya disampaikan ke sebuah pengolah data. Karena tidak seluruh cermin diaktifkan untuk menangkap bayangan benda, jumlah piksel yang direkam lebih kecil jika dibandingkan dengan kamera dijital biasa. Selanjutnya rekonstruksi dilakukan menurut prinsip pencuplikan kompresif.

Dari sudut pandang pencuplikan kompresif, operasi pengambilan elemen bayangan secara acak oleh DMD adalah proses pengamatan yang dapat dinyatakan sebagai proyeksi s dengan matriks pengukur . Jika N jumlah total dari cermin DMD, sedangkan M adalah cermin aktif dalam satu kali pengambilan gambar, maka akan berupa sebuah matriks berukuran MxN. Setiap baris dari matriks ini akan berisi satu buah angka satu yang letaknya acak, sedangkan elemen sisanya pada baris ini bernilai nol. Operasi pengukuran ini dapat dituliskan sebagai:

(1)

Berdasarkan persamaan diatas, proses pengamatan dapat dipandang sebagai pencuplikan kompresif dengan dua cara yang berbeda, yaitu:

1. Pencuplikan kawasan spasial terhadap sinyal s secara acak
2. Pencuplikan kawasan transform sebagian dari basis penjarang secara acak

Dua cara pandang ini dapat memberikan dua pendekatan berbeda dalam mencari solusi atau merekonstruksi citra hasil pengamatan. Tulisan ini akan menjelaskan terlebih dahulu cara pertama, sedangkan cara kedua akan dijelaskan pada tulisan berikutnya.

Pencuplikan kompresif kawasan spasial

Sebelumnya telah diperkenalkan formulasi teknik pencuplikan kompresif dalam kawasan ruang-waktu. Secara prinsipil, teknik pencitraan yang pertama ini tidak berbeda dari yang telah dijelaskan dalam pencuplikan kompresif. Dengan demikian, diagram ranah-jamaknya masih tetap sama seperti sebelumnya.

Operasi proyeksi dilakukan dengan mengalikan sinyal asli s dengan operator proyeksi yang dibuat dengan ”mengaktifkan” satu elemen saja dari setiap vektor baris matriks ini. Setelah mengubah ”bayangan” benda dari larik dua dimensi menjadi satu dimensi, instruksi untuk mensimulasikan proses pengukuran dengan Matlab adalah:

% *** Simulasikan Digital Micromirror Device ***

rp=randperm(N); act_cell(1:M_sub)=rp(1:M_sub); %randomize mirror in DMD

PHI=zeros(M_sub,N);

for k=1:M_sub; PHI(k,act_cell(k))=1;end;

PHI=orth(randn(M_sub,N)’)'>0;

%*** Lakukan pengukuran x_sub=PHI*x0 ***

x_sub=PHI*x0;

Pada contoh ini akan diamati sebuah objek sederhana, yaitu citra yang menyerupai huruf “H” seperti yang diperlihatkan pada Gambar 2. Hasil pengukuran akan merupakan citra tak-sempurna karena piksel yang terambil hanya sedikit dan letaknya acak. Hasil ini diperlihatkan pada Gambar 3.

Gambar 2. Citra sederhana dalam pencitraan kompresif

Gambar 3. Citra hasil pengukuran

Berdasarkan hasil pengukuran, gambar dapat direkonstruksi dengan melalui algoritma basis pursuit. Sebagai contoh bisa diambil DCT sebagai transformasi penjarang. Jika menyatakan hasil pengamatan, vektor basis penjarang ke-k, maka pencarian basis dapat dirumuskan sebagai berikut:

P₁: min ||S||₁ s.t. (2)

Setelah seluruh vektor basis ditemukan, citra bisa direkonstruksi dengan menjumlahkan vektor-vektor basis ini yang diboboti dengan S_k yang sesuai. Hasil rekonstruksi ditunjukkan pada Gambar 4.

Gambar 4. Citra rekonstruksi

Meskipun citra huruf ”H” masih terlihat pada hasil rekonstruksi, Gambar 4 memperlihatkan bahwa secara visual hasilnya masih kurang baik. Hal ini bisa dijelaskan dengan melihat seluruh komponen S_k yang terambil pada pencarian basis.

Gambar 5. Magnitudo koefisien terpilih yang diurutkan

Gambar 5 memperlihatkan kurva magnitudo S_k yang telah diurutkan. Meskipun koefisien ini meluruh, terlihat bahwa peluruhan ini kurang cepat. Pengambilan koefisien sebanyak derajat sparsity yang diperkirakan masih akan menyisakan informasi yang cukup penting. Dengan demikian, masalah utamanya adalah pemilihan transformasi penjarang yang kurang tepat. Bagaimana cara mengatasinya?

Daftar Pustaka

1. Michael Wakin, Jason Laska, Marco Duarte, Dror Baron, Shriram Sarvotham, Dharmpal Takhar, Kevin Kelly, and Richard Baraniuk, “An Architecture for Compressive Imaging,” Proc. International Conference on Image Processing — ICIP 2006, Atlanta, GA, Oct. 2006.
2. Dharmpal Takhar, Jason Laska, Michael Wakin, Marco Duarte, Dror Baron, Shriram Sarvotham, Kevin Kelly and Richard Baraniuk, “A New Compressive Imaging Camera Architecture using Optical-Domain Compression,” Proc. of Computational Imaging IV at SPIE Electronic Imaging, San Jose, CA, Jan. 2006.
3. Marco F. Duarte, Mark A. Davenport, Dharmpal Takhar, Jason N. Laska, Ting Sun, Kevin F. Kelly and Richard G. Baraniuk, “Single Pixel Imaging via Compressive Sampling”, IEEE Signal Processing Magazine, March 2008.

Kode Matlab

%———————————————————-

% Pencitraan Kompresif: kamera dijital dng 1 foto-dioda

%———————————————————-

clear;%close all;clc;

% read original image

%cam01=double(imread(‘cameraman.tif’,'tif’));

%I=cam01(140:155,35:50); %load cam01; I=cam01;

%baca citra

I=double(imread(‘IMG\H.bmp’,'bmp’));

%tampilkan citra asal;

figure(1); imagesc(I); colormap(gray); title(‘citra asli’);

[N_col,N_row]=size(I);

N=N_col*N_row; %panjang sinyal

%Pembentukan basis sparsity /DCT

PSI=dct(eye(N,N))’; %buat matriks DCT

K_sparse=2*round(sqrt(N));%tingkat sparsity

x = reshape(I,N,1); %ubah menjadi vektor

x0=(x-mean(x))/(max(x)-min(x));%normalisasi

%proses pengamatan/pengukuran

%estimasikan minimal jumlah cuplikan M>CKlog(N)

mu=1*(1/sqrt(N)); %tentukan koherensi->perkiraan

C=sqrt(N)*mu; %koherensi tak ternormalisasi=sqrt(N).koh

M_sub=ceil(C*K_sparse*log(N));

% *** Simulasikan Digital Micromirror Device ***

rp=randperm(N); act_cell(1:M_sub)=rp(1:M_sub); %randomize mirror in DMD

PHI=zeros(M_sub,N);

for k=1:M_sub; PHI(k,act_cell(k))=1;end;

PHI=orth(randn(M_sub,N)’)'>0;

% *** Lakukan pengukuran x_sub=PHI*x0 ***

x_sub=PHI*x0;

x_obs=zeros(N,1);x_obs(act_cell)=x0(act_cell);

%Lakukan recovery

%Bentuk kamus basis D=PHI*PSI

D=PHI*PSI;

%x_sub=D*x0;

figure(4);imagesc(reshape(x_obs,N_row,N_col));colormap(gray);title(‘DMD aktif’);

[D_row,D_col]=size(D);

D_pos=-D;

%ubah L1 menjadi LP

A=zeros(D_row,2*D_col);

A(:,1:D_col)=D;

A(:,D_col+1:2*D_col)=D_pos;

f=ones(2*D_col,1);

lb=zeros(2*D_col,1);

f1=ones(2*D_col,1);

lb1=zeros(2*D_col,1);

A1=-A;

x_sub1=x_sub;

%*******************************

%* Jalankan algoritma simpleks *

%*******************************

alpha0=linprog(f,[],[],A,x_sub,lb); %get positive part

alpha1=linprog(f1,[],[],A1,x_sub1,lb1);%get negative part

alpha=alpha0-alpha1;

% Urutkan

[z,imax]=sort(abs(alpha(1:D_col)),’descend’);

x_abs=1:D_col;

figure(2); plot(x_abs, z,’b-’);

%Reconstrusi Citra

%===

x_hat=zeros(N,1);

for k=1:K_sparse;

x_hat=x_hat+alpha(imax(k))*PSI(:,imax(k));

end;

Ip = reshape(x_hat, N_row, N_col);

figure(3);imagesc(Ip);colormap(gray);title(‘citra rekonstruksi’);

Pencuplikan Kompresif (2)

Ada dua transformasi penting dalam pencuplikan kompresif, yaitu sparsity transform dan projection transform. Transformasi pertama dipakai untuk mencari komponen sparse dari sinyal, sedangkan yang kedua dipakai dalam operasi pengukuran atau pengamatan. Tulisan ini memberikan contoh sederhana cara melakukan pencuplikan kompresif dan rekonstruksi sinyal secara eksak dari cuplikan terkompresi.

Aspek Penting dalam Pencuplikan Kompresif

Pada pembahasan mengenai dekomposisi sinyal dan pencarian basis ideal, telah dijelaskan cara menemukan vektor basis penyusun suatu sinyal hasil sintesis K buah vektor kamus basis , atau sebuah sinyal K-sparse. Berdasarkan prinsip ketidakpastian yang diperumum, jika jumlah komponen K kurang dari dengan menyatakan koherensi dari kamus basis, maka semua komponen penyusun sinyal akan dapat ditemukan oleh optimisasi terkendala P₁ (basis pursuit). Pada akhir pembahasan ditunjukkan, bahkan ketika sinyal sudah tercampur derau, seluruh komponen masih dapat ditemukan sehingga sinyal dapat direkonstruksi secara eksak.

Pencuplikan kompresif bekerja berdasarkan prinsip yang mirip seperti diatas, yakni: diberikan informasi M buah hasil pengamatan dari sinyal K-sparse, sinyal asal sepanjang N>>M>K akan dapat direkonstruksi jika dengan C bergantung pada sistem basis yang dipilih. Hal yang perlu dicermati disini, jumlah cuplikan yang diperlukan menjadi jauh lebih kecil daripada N, bahkan hanya dalam orde logaritma-nya. Disamping itu, nilai M juga dipengaruhi oleh tingkat sparsity K dari sinyal yang dapat dianggap menggambarkan kandungan informasi dari sinyal yang sebenarnya. Pemilihan kamus-basis akan menentukan kinerja kompresi karena dapat menentukan besar-kecilnya jumlah data cuplikan yang diperlukan. Berdasarkan penjelasan ini, ada beberapa hal penting yang perlu diperhatikan dalam melakukan pencuplikan secara kompresif.

Pertama, sinyal yang akan dicuplik haruslah bersifat K-sparse. Kebanyakan sinyal hasil pengamatan di dunia nyata pada dasarnya adalah sinyal sparse terhadap suatu basis tertentu, meskipun belum tentu demikian di kawasan sinyal tersebut diukur. Sistem basis atau transformasi yang membuat suatu sinyal menjadi bersifat sparse disebut sebagai transformasi penjarang (sparsity transform). Termasuk dalam kategori ini antara lain, berbagai sistem transformasi ortogonal yang telah disinggung pada tulisan-tulisan sebelumnya seperti, DCT, Hadamard, DFT, wavelet, … dll.

Aspek penting kedua adalah observasi atau pengukuran (observation/ measurement). Proses pencuplikan klasik dapat dianggap sebagai proses observasi dari suatu sinyal kontinyu. Secara umum, observasi dapat dinyatakan sebagai proyeksi sinyal pada basis tertentu. Pencuplikan klasik adalah proyeksi dari sinyal kontinyu ke basis impuls, sistem pencitraan SFCW-GPR (stepped-frequency continuous wave ground penetrating radar), MRI (magnetic resonance imaging), dan teleskop VLBI (very large base line interferometry) melakukan observasi data spasio-temporal didalam basis Fourier, sedangkan CT (computed tomography) scanner melakukan observasi dalam kawasan Radon. Sistem basis atau transformasi dimana suatu sinyal diamati disebut sebagai transformasi proyeksi (projection transform).

Yang terakhir, besarnya konstanta C yang akan menentukan jumlah minimum data observasi M berkaitan erat dengan derajat koherensi sistem basis. Kedua macam transformasi yang disebutkan akan menentukan besarnya nilai ini karena yang diukur adalah koherensi antara transformasi penjarang dengan transformasi proyeksi.

Gambar 1. Diagram ranah-jamak [1] dalam sistem pencuplikan kompresif sederhana

Hubungan berbagai transformasi yang terlibat dalam pencuplikan kompresif dilukiskan sebagai diagram ranah-jamak [1] pada Gambar 1. Sinyal s sepanjang N-cuplikan bersifat K-sparse dalam sistem basis . Pencuplikan atau pengamatan akan menghasilkan sinyal yang hanya mengandung sebagian data dari sinyal asalnya karena ada informasi yang hilang, yaitu . Jika jumlah pengamatan mencukupi, meskipun jauh dibawah N, K-buah basis sparse dapat ditemukan melalui P₁. Dengan demikian, karena seluruh basis S dapat ditemukan, maka sinyal asal dapat direkonstruksi secara eksak. Pada diagram ini, transformasi penjarang adalah , sedangkan transformasi proyeksi adalah .

Pencuplikan Kompresif Sederhana

Sebagai ilustrasi, akan ditinjau suatu sinyal s yang bersifat sparse terhadap basis DCT dengan derajat sparsity sebesar 4. Bagaimana cara melakukan pencuplikan kompresif dari sinyal ini dan cara melakukan rekonstruksinya?

Pada kasus ini bisa dipilih suatu basis acak Gaussian yang dinormalisasi sebagai transformasi proyeksi-nya. Karena pencuplikan langsung dari sinyal analog melibatkan beberapa hal teknis yang baru akan dibahas pada tulisan mendatang, anggap bahwa pengamatan ini dapat direpresentasikan sebagai subsampling terhadap sinyal waktu diskrit s sepanjang N-cuplikan. Dengan demikian, tahap pertama yang harus dilakukan adalah membentuk matriks pengamatan atau transformasi proyeksi . Dimensi matriks ini adalah MxN, dimana M<<N adalah jumlah pengamatan yang diharapkan. Berdasarkan teorema pencuplikan kompresif, agar sinyal s bisa direkonstruksi secara eksak, maka jumlah ini minimal-nya adalah sekitar CK log(N).

Berbagai peneliti bidang ini memberi nilai C yang berlainan. Pada umumnya nilai ini sebesar koherensi kamus-basis (yang tak-ternormalisasi), dimana kamus basis adalah perkalian antara matriks proyeksi dengan matriks sparsity, sehingga . Pada percobaan ini akan diambil C=1, dengan demikian yang kemudian menentukan M adalah derajat sparsity (K=4) dan logaritma dari cuplikan Nyquist N. Sebagai contoh diambil N=64, jadi M~17 yang berarti adalah faktor kompresi sekitar 3.8 kalinya.

Gambar 2. Sinyal asli dan sinyal cuplikan kompresif

Gambar 2 menunjukkan sinyal s dan cuplikan kompresif-nya . Terlihat bahwa hasil pencuplikan kompresif memberikan data yang lebih ringkas atau lebih pendek dari sinyal asalnya. Dalam skema pencuplikan yang sedang dibahas, sinyal ini adalah .

Gambar 3. Vektor basis terpilih hasil optimisasi P₁

Rekonstruksi dilakukan dengan terlebih dahulu mencari vektor sparse S dalam kawasan DCT melalui optimisasi P₁:

(P₁): min||S||₁ s.t.

Optimisasi ini memiliki arti: pilihlah vektor S paling sparse dalam basis yang dapat menjelaskan hasil pengamatan . Hasilnya ditampilkan dalam Gambar 3, dimana semua basis terpilih tepat sama dengan basis penyusun sinyal.

Gambar 4. Perbandingan sinyal asli dan hasil rekonstruksi

Rekonstruksi dilakukan dengan menyusun kembali sinyal dari vektor basis terpilih. Hasilnya ditampilkan pada Gambar 4, dimana sinyal rekonstruksi identik dengan sinyal asalnya dan sinyal kesalahan bernilai nol.

Pemilihan basis acak sebagai transformasi proyeksi sangat menguntungkan jika dilihat dari beberapa sisi. Yang pertama, setiap data cuplikan didalam s dapat dianggap sama pentingnya. Dengan demikian tidak perlu skema pencuplikan adaptif dengan memilih komponen dominan saja karena semuanya sama penting. Yang kedua, basis acak bersifat universal karena hampir ortogonal ke semua basis lain. Basis ini bisa diterapkan ke berbagai kasus pencuplikan kompresif. Dan yang terakhir, modulasi dengan basis acak untuk mendapatkan membuat hasil pengamatan sedikit terenkripsi (weakly encrypted). Hal ini sangat menguntungkan jika pencuplikan kompresif akan dipakai didalam sistem sensor tersebar kompresif (distributed compressed sensing) dan data yang diamati bersifat rahasia.

Daftar Pustaka

1. AB Suksmono, “A graphical representation of the multi-domain signal processing,” Proc. of SICE-ICASE 2006, Busan, Korea.

Program Matlab

%———————————————–
%Contoh Pencuplikan Kompresif
%———————————————–
clear;%close all;clc;
N=64; %panjang sinyal
%Pembentukan sinyal 4-sparse
PSI=dct(eye(64,64)); %buat matriks DCT
K_sparse=4;%tingkat sparsity
%pN=randperm(nAtoms);
pN=[3 13 42 58];
x=zeros(N,1);
X=zeros(N,1);%’transformed-domain signal
for k=1:K_sparse;
x=x+(-1)^k*PSI(:,pN(k)); %masukkan faktor polaritas basis
X(pN(k))=1;
end;
x0=x;

%lakukan pengukuran atau subsampling
%estimasikan minimal jumlah cuplikan M>CKlog(N)
mu=1*(1/sqrt(N)); %tentukan koherensi->perkiraan
C=sqrt(N)*mu; %koherensi tak ternormalisasi=sqrt(N).koh
M_sub=ceil(C*K_sparse*log(N))
%bentuk matriks pengukuran
PHI=orth(randn(M_sub,N)’)';
%lakukan pengukuran x_sub=PHI*x
x_sub=PHI*x;
x_sub_xtd=zeros(N,1); x_sub_xtd(1:M_sub)=x_sub(1:M_sub);
figure(1);plot(1:N, x0,’k-’,1:N,x_sub_xtd,’b:’);title(’subsampling’);
legend(‘original’,’sub-samples’); xlabel(‘time’);ylabel(‘amplitude’);
%Lakukan recovery
%Bentuk dictionary D=PHI*PSI
D=PHI*PSI;

[D_row,D_col]=size(D);
D_pos=-D;
%define dictionary A
A=zeros(D_row,2*D_col);
A(:,1:D_col)=D;
A(:,D_col+1:2*D_col)=D_pos;

%Pecahkan P1
f=ones(2*D_col,1);
lb=zeros(2*D_col,1);

%*******************
%alpha1=alpha;
f1=ones(2*D_col,1);
lb1=zeros(2*D_col,1);
A1=-A;
x_sub1=x_sub;
%*******************

alpha0=linprog(f,[],[],A,x_sub,lb); %ambil solusi positif
alpha1=linprog(f1,[],[],A1,x_sub1,lb1); %ambil solusi negatif

alpha=alpha0-alpha1; %gabungkan keduanya
%Tampilkan hasil
figure(2);plot(1:D_col,X,’r.’,1:D_col,abs(alpha(1:D_col)),’b-’);
title(‘extracted dct vectors by P1′);
legend(‘original components’,'extracted components’);
xlabel(‘DCT vector index’);ylabel(‘magnitude’);
%obtain a few largest component, by sorting alpha
[z,imax]=sort(abs(alpha(1:D_col)),’descend’);
%Reconstruct the estimated values
x_hat=zeros(N,1);
for k=1:K_sparse;
x_hat=x_hat+alpha(imax(k))*PSI(:,imax(k));
end;
figure(3);plot(1:N,x0,’k.’,1:N,x_hat,’r:’,1:N,x0-x_hat,’b-’);
xlabel(‘time’);ylabel(‘amplitude’);
legend(‘original’,'reconstructed’,'error’);

Pencuplikan Kompresif (1)

Pencuplikan adalah salah satu proses terpenting dalam tahapan pengolahan sinyal dijital, dimana sinyal analog diubah menjadi sinyal diskrit (dijital). Teorema pencuplikan Shannon menetapkan aturan batas minimal laju pencuplikan agar sinyal dapat direkonstruksi ke asalnya, yakni sebesar dua kali lebar-pita dari sinyal. Berdasarkan paradigma pencuplikan kompresif, aturan ini dapat dilanggar karena yang berperanan dalam menentukan rekonstruktivibilitas bukanlah kandungan frekuensi sinyal, melainkan kandungan informasinya yang tercermin dalam derajat sparsity.

Batas Nyquist untuk Pencuplikan Sinyal

Didalam Pengolahan sinyal dijital, pencuplikan (sampling) adalah proses mengambil cuplikan dari sinyal analog secara periodik. Proses ini dapat dinyatakan sebagai perkalian antara sinyal analog dengan deretan impuls , dimana n adalah bilangan bulat dan T adalah perioda pencuplikan. Besaran yang lebih populer dibandingkan dengan T, yaitu kecepatan pencuplikan atau fekuensi pencuplikan .

Seperti yang tergambar didalam blok diagram pengolahan sinyal dijital, pita frekuensi dari sinyal akan terlebih dahulu dibatasi (bandlimited) dengan sebuah tapis anti-aliasing sebelum pencuplikan dilakukan. Secara matematis, proses pencuplikan dapat dirumuskan sebagai perkalian antara (modulasi) sinyal analog dengan sebuah deretan impuls

(1)

Gambar 1 menjelaskan proses pencuplikan yang dilihat dari kawasan waktu. Tinjauan kawasan waktu berguna untuk melihat terbentuknya sinyal waktu diskrit (dijital) dan aspek pencuplikan selang-seragam (uniform sampling) terhadap keluarannya. Namun demikian, untuk mengetahui persyaratan laju pencuplikan minimum atau selang pencuplikan maksimum agar sinyal s_a(t) dapat diperoleh kembali, maka proses ini perlu dilihat dari kawasan frekuensi.

Gambar 1. Pandangan proses pencuplikan dari kawasan waktu.

Ekspresi suatu sinyal dalam kawasan frekuensi dapat diperoleh melalui transformasi Fourier dari persamaan dalam kawasan waktu-nya. Karena perkalian dalam kawasan waktu ekivalen dengan konvolusi dalam kawasan frekuensi, maka hasil transformasi Fourier dari persamaan (1) adalah

(2)

Gambar 2 memperlihatkan proses pencuplikan jika dilihat dari kawasan frekuensi. Karena transform Fourier dari deretan impuls adalah juga suatu deretan impuls, maka konvolusi antara spektrum sinyal dengan impuls menghasilkan pergeseran spektrum sejauh . Sebagai akibatnya akan terjadi pengulangan (tiling) spektrum di seluruh rentang frekuensi pada posisi kelipatan dari frekuensi pencuplikan. Gambar 2 bagian kiri bawah menunjukkan spektrum dari sinyal yang lebar pitanya yang kemudian mengalami proses tiling akibat proses pencuplikan.

Gambar 2. pencuplikan dilihat dari kawasan frekuensi

Jika jarak antar pengulangan atau grid tiling cukup lebar, seperti diperlihatkan pada Gambar 3 bagian atas, yang juga berarti bahwa frekuensi pencuplikannya cukup besar, maka tidak akan terjadi tumpang tindih antar spektrum yang bertetangga. Kondisi ini disebut sebagai non-aliasing. Selanjutnya sifat keunikan dari transformasi Fourier akan menjamin bahwa sinyal asal dapat diperoleh secara sempurna. Sebaliknya, jika kurang besar, maka akan terjadi tumpang tindih antar spektrum yang mengakibatkan hilangnya sebagian dari informasi. Peristiwa ini disebut aliasing, seperti diperlihatkan pada Gambar 3 bagian bawah. Pada kondisi ini, sinyal tidak dapat lagi direkonstruksi secara eksak.

Gambar 3. Kondisi non-aliasing dan aliasing pada proses pencuplikan

Dengan memahami peristiwa aliasing dalam kawasan frekuensi, maka batas minimum laju pencuplikan atau batas Nyquist dapat diperoleh, yaitu sebesar . Hasil ini durumuskan sebagai teorema Shannon untuk pencuplikan sebagai berikut:

Teorema Pencuplikan Shannon. Suatu sinyal pita-terbatas dengan lebar dapat direkonstruksi secara eksak dari cuplikannya jika laju pencuplikan minimum dua kali dari lebar pita tersebut, atau .

Sebagai contoh, manusia dapat mendengar suara dari frekuensi 20 Hz sampai dengan sekitar 20kHz, artinya lebar pita dari suara yang mampu didengar manusia adalah sekitar 20 KHz. Dengan demikian, pengubahan suara menjadi data dijital memerlukan laju pencuplikan sedikitnya 2×20kHz = 40 kHz atau 40.000 cuplikan/detik supaya sinyal suara dapat direkonstruksi secara sempurna, yang berarti juga kualitas dari suara hasil perekaman dijital dapat dimainkan tanpa distorsi.

Pada teori pencuplikan klasik, rekonstruksi sinyal s_r(t) akan ditempuh dengan cara melakukan konvolusi antara sinyal disktrit s(n) hasil cuplikan, dengan tapis rekonstruksi

(3)

Dengan demikian,

(4)

adalah sinyal terekonstruksi yang diinginkan. Tanggapan impuls dari tapis rekonstruksi pada persamaan (3) diperlihatkan pada Gambar 4 berikut ini.

Gambar 4. Tanggapan impuls dari tapis rekonstruksi

Pencuplikan Sub-Nyquist dengan Compressive Sensing

Teorema pencuplikan Shannon memberikan batas minimal laju pencuplikan atau jumlah cuplikan persatuan waktu supaya sinyal kontinyu dapat direkonstruksi secara eksak. Pencuplikan dibawah laju ini akan mengakibatkan terjadinya distorsi atau cacat aliasing. Apakah ada cara untuk memperkecil batas ini tanpa mengakibatkan distorsi?

Teorema pencuplikan diatas didasarkan pada analisis Fourier, dimana fungsi basis dalam transformasi sinyal berupa fungsi sinusoid dari berbagai frekuensi. Bahkan, makna sebenarnya dari transformasi Fourier adalah penguraian suatu sinyal kedalam basis sinusoid. Batas Nyquist menjaga supaya semua komponen frekuensi tetap utuh didalam hasil cuplikan sehingga sinyal dapat dikembalikan ke bentuk asalnya. Jika basis selain sinusiod dipakai dalam mendekomposisi sinyal, tentu aturan ini akan berubah karena yang sebenarnya menggambarkan kandungan informasi dari sinyal bukanlah lebar pita fekuensi-nya.

Gambar 5. Pencuplikan dan rekonstruksi sebuah impuls

Kaitan antara pencuplikan dan kandungan informasi diperlihatkan pada Gambar 5. Sebuah impuls akan memiliki transform Fourier S(k) yang tersebar ke seluruh kawasan frekuensi. Berdasarkan informasi kandungan frekuensi sinyal, maka seluruh komponen frekuensi didalam S(k) diperlukan untuk merekonstruksi sebuah impuls. Penghilangan sebagian dari cuplikan, dalam gambar diperlihatkan hanya 25% nya saja, akan menghasilkan rekonstruksi yang tak sempurna. Sinyal impuls yang telah terdistorsi akibat hilangnya sebagian dari komponen frekuensi diperlihatkan pada Gambar 5 bagian bawah.

Seandainya diketahui ada satu buah impuls di kawasan waktu diskrit dalam selang 32 cuplikan, berapa banyakkah informasi yang dibutuhkan untuk mendeskripsikan sinyal ini?

Gambar diatas menyatakan bahwa seluruh 32 buah komponen frekuensi perlu. Tetapi yang sebenarnya diperlukan hanya satu buah saja, yaitu posisi dari impuls. Informasi inilah yang sesungguhnya lebih esensial dibandingkan dengan lebar pita frekuensi. Satu buah informasi, dalam hal ini posisi impuls, disebut juga sebagai derajat kebebasan dari sinyal s(n). Dalam terminologi CS, besaran ini dinamakan tingkat sparsity dari sinyal. Dengan demikian, sinyal s(n) diatas adalah sebuah sinyal 1-sparse. Hasil penting dari pencuplikan kompresif dirumuskan dalam teorema berikut ini:

Teorema Pencuplikan Kompresif. Suatu sinyal K-sparse sepanjang N-cuplikan dapat direkonstruksi secara eksak berdasarkan buah cuplikannya, dimana C sebuah konstanta kecil yang nilainya bergantung pada sistem-basis yang digunakan.

Jumlah M ini jauh lebih kecil daripada N didalam teknik pencuplikan klasik. Inilah alasan mengapa pencuplikan kompresif dapat dianggap sebagai pencuplikan Sub-Nyquist. Bagaimana melakukan rekonstruksi sinyal jika cuplikannya jauh lebih kecil dari batas Nyquist? Yang jelas, proses rekonstruksi klasik melalui konvolusi sinyal pengamatan dengan tapis rekonstruksi tidak lagi bisa dipakai.

Referensi

1. D. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory, 52, 2006, pp.1289-306.
2. E. J. Candes and J. Romberg, “Sparsity and incoherence in compressive sampling,” Inverse Problems, 23 (2007), pp.696-985.
3. E.J. Candes, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory, Vol. 52, No. 2, Feb. 2006, pp. 489-509.

Kode Matlab

Kode Matlab untuk membuat kurva tanggapan impuls dari tapis rekonstruksi

%MATLAB codes to generate hr(t)

%define sampling period

T = 0.1 ;%0.1 seconds, sampling at 10 Hz

%define time axis

t=-1.0:0.001:1.0; %2 second

hr=sin(pi*t/T)./(pi*t/T);

figure(1);plot(10*t,hr);

xlabel(‘time in Ts’); ylabel(‘hrec’)

title(‘ideal reconstruction filter’);

Kode Matlab untuk melihat pengaruh penghilangan sebagian cuplikan terhadap hasil rekonstruksi sinyal impuls.

%generate an impulse in N=32 sample intervals;

clear;close all; clc;

N=32; s=zeros(N,1);

s(16)=1; %put impulse in the middle

figure(1);stem(s);title(‘Time Domain’);

xlabel(‘time’); ylabel(‘magnitude’);

S=fft(s); %transform s into freq-domain S

figure(2); stem(abs(S));title(‘Frequency Domain’);

xlabel(‘time’); ylabel(‘magnitude’);

%discard some 25% of the components

S1=S; for k=1:round(N/4); S1(4*k)=0;end;

s1=ifft(S1); %signal inversion by IDFT

figure(3); stem(abs(S1));title(‘Frequency Domain’);

xlabel(‘time’); ylabel(‘magnitude’);

figure(4);stem(abs(s1));title(‘Time Domain’);

xlabel(‘time’); ylabel(‘magnitude’);

Dongeng dari Negeri Tanahjarang

oooOOOooo

Tanahjarang adalah nama sebuah negeri yang diperintah oleh seorang Raja yang pandai, adil, bijaksana, dan selalu memikirkan kesejahteraan dan kebahagiaan rakyatnya.

Musim kemarau tahun ini membuat hasil pertanian menurun tajam hingga mengancam kesejahteraan negeri. Raja yang bijak sudah lama memperkirakan kejadian seperti ini dan bersiap menghadapinya dengan membuat lumbung padi, gudang penyimpan ikan dan daging yang diawetkan, dan sebuah fasilitas istimewa untuk menyimpan berbagai sayuran dan buah-buahan.  Selama beberapa musim panen belakangan, rakyat telah diminta menyisihkan hasil sawah dan kebun mereka ke gudang makanan untuk menghadapi masa sulit yang akan segera tiba.

Ketika akhirnya paceklik datang, Kerajaan Tanahjarang dan seluruh rakyat di penjuru negeri siap menghadapinya. Raja yang pandai tahu bahwa masa sulit akan berlangsung lama. Dengan demikian, makanan didalam lumbung dan gudang harus dibagikan dengan sehemat mungkin. Tapi Raja yang bijak juga tahu pentingnya nutrisi yang lengkap untuk rakyatnya, supaya anak negeri ini kelak tumbuh menjadi pemuda yang sehat dan cerdas. Menurut penasihat kerajaan, rakyat tidak hanya cukup diberi nasi saja, tetapi juga lauk-pauk, buah-buahan, dan sayur-sayuran. Bukan hanya itu, setiap desa di pelosok negeri memiliki masakan khas masing-masing, meskipun semuanya terbuat dari bahan pokok yang sama. Ada desa yang suka masakan pedas, masakan asam, masakan manis, masakan asin, dan lain-lain rasa yang jenisnya sebanyak jumlah desa di negeri ini.

Raja ingin selalu membuat rakyat bahagia, bahkan dalam hal selera makan mereka. Karena ketersediaan bahan pangan yang terbatas, tidak mungkin memberikan langsung bahan-bahan pokok ini ke masing-masing penduduk. Sang Raja juga telah belajar dari sahabatnya, Maharaja Hindia, bahwa masalah pembagian makanan bukanlah hal yang sederhana. Lumbung padi di pelosok Negeri Hindia pernah habis karena sang Maharaja ceroboh bertaruh dengan seorang pecatur ulung yang berhasil mengalahkannya; meskipun hadiah yang dimintanya nampak sepele: berikan satu butir beras ke kotak catur pertama, dua butir ke kotak kedua, empat butir di kotak ketiga, …, dan seterusnya sampai ke kotak ke-enampuluh empat. Permintaan sederhana ini telah menyedot persediaan beras di seluruh negeri. Jadi, jika tidak hati-hati,  maka seluruh Negeri Tanahjarang bisa tertimpa bencana.

Raja yang bijak akhirnya memutuskan bahwa rakyat tidak akan diberikan bahan makanan, melainkan masakan yang sudah jadi. Dengan demikian, seluruh juru masak yang mewakili desa dipanggil ke Kerajaan Tanahjarang untuk secara bersama-sama mempersiapkan masakan bagi penduduk desa masing-masing. Hal ini dilakukan supaya tidak ada bahan makanan yang terbuang percuma karena takaran adonan untuk setiap jenis masakan berlainan.

Beberapa bulan kemudian musim kemarau berlalu. Hujan kembali turun, tanah menjadi subur. Petani kembali kesawah menanam padi, peternak menggiring sapi ke ladang mencari rumput. Penduduk Negeri Tanahjarang bekerja dengan penuh sukacita. Masa sulit itu telah berhasil mereka lalui berkat Raja yang bijak dan dicintai rakyatnya.

oooOOOooo

Apa hubungan dongeng diatas dengan compressive sensing … ? Legenda menyebutkan bahwa Sang Raja adalah keturunan langsung dari Maharaja Heisenberg, seorang Raja yang berhasil menjadikan negerinya makmur sejahtera dengan membuat aturan ketat tentang pembagian jatah makanan di masa paceklik.

Dekomposisi Sinyal dan Pencarian Basis Ideal (4)

Masalah pencarian vektor basis ideal dari sebuah sinyal telah dapat diformulasikan sebagai proses optimisasi berkendala (constrained optimization) P₁. Tulisan ini akan terlebih dahulu membahas secara sekilas konsep optimisasi linier, khususnya LP (Linear Programming), dan metode pemecahannya. Meskipun kendala dalam P₁ 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 dan yang meminimumkan nilai , dengan syarat bahwa , , dan .

Masalah diatas dapat dituliskan kedalam bentuk baku sebagai berikut:

minimize

subject to:

,

,

.

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

Gambar 1. Ilustrasi pencarian titik optimum dalam LP

Pada Gambar 1 diatas, daerah yang diarsir adalah pasangan 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 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

subject to:

bahkan dengan dengan menggunakan slack variables

, i = 1, 2, …, m

LP dapat dituliskan dalam bentuk matriks yang lebih ringkas, yaitu

minimize subject to

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 .

Metode Basis Pursuit dengan Pemrograman Linier

Metoda basis pursuit merumuskan pencarian basis ideal didalam kamus-basis berdimensi NxM sebagai masalah pemrograman konveks:

(P₁):

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

jika dilakukan substitusi sebagai berikut

M↔2M; x↔ ; c↔[1 1 …. 1; 1 1… 1]; A↔ ; 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 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.

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).

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.