1/28/2015

Independent Component Analysis (ICA) menggunakan Octave

Di signal processing, ICA adalah salah satu metode komputasi numerik untuk memisahkan sumber-sumber independen atau sinyal-sinyal independen yang tercampur secara linier dan sinyal tersebut direkam oleh beberapa sensor. Ada banyak algoritma yang digunakan untuk ICA dan salah satu algoritma yang efisien and cukup populer adalah FastICA yang ditemukan oleh Aapo Hyvärinen. Algoritma ini dibangun  berdasarkan fixed-point iteration scheme dengan memaksimalkan non-Gaussianity sebagai parameter untuk mengukur independensi suatu sinyal dan dapat juga diturunkan dari pendekatan iterasi newton.

Cukup komplek jika kita harus menurunkan persamaan ICA, namun cukup mudah untuk menuliskannya dalam Octave. Bahasan Kali ini, kita akan mencoba menuliskan algoritma ICA menggunakan Octave dan mengujinya dengan memisahkan sinyal-sinyal yang tercampur secara linier.

Model Pencampuran
 
Persamaan di bawah ini menunjukkan model pencampuran (mixing model) di mana $s$ adalah sumber-sumber sinyal (sebelum dicampur), $A$ adalah matrik pencampuran, dan $x$ adalah sinyal-sinyal yang tercampur.

$ x = A s $

Untuk mendapat $s$ cukup mudah seperti persamaan di bawah ini

$ s = A^{-1} x $

Namun akan menjadi sangat sulit untuk mendapatkan $s$ jika kita hanya mengetahui $x$ saja tanpa mengetahui $A$. Nah fungsi ICA adalah untuk mengestimasi nilai $A$ untuk mendapatkan $s$ dengan hanya menggunakan informasi dari $s$.

Ada 5 step untuk menyusun algoritma ini di Octave:

Pertama, Centering
Secara sederhana bisa dikatakan sebagai proses untuk membuat rataan suatu data atau sinyal menjadi nol yakni dengan mengurangi sinyal tersebut dengan rataann dirinya sendiri. Hasil pengurangan tersebut yakni sinyal baru yang memiliki rataan nol (zero-mean).

$ x_c = x - \bar{x} $

di mana $x_c$ adalah $x$ setelah proses centering, dan $\bar{x}$ adalah rataan dari $x$.

Kedua, Whitening
Yakni proses transformasi data atau sinyal ke bentuk data baru di mana kovarian matrixnya adalah matrik identitas. Metode yang biasa digunakan untuk proses whitening adalah eigen-value decomposition (EVD).

$ x_w = E \ D^{-\frac{1}{2}} \ E^T \ x_c $

dimana $E$ adalah orthogonal matrix dari eigenvector dan $D$ adalah diagonal matrix dari eigenvalue.

Ketiga, Algoritma FastICA
   a) tentukan nilai $w_n$ secara random
   b) update atau perbaharui nilai $w_n$ berdasarkan rumus berikut

 $w^{+} = E\{x_w g(w_n^T x_w)\} - E\{g^{'}(w_n^T  x_w)\}w$


$ w_{n+1} = (w^{+} w^{+T})^{-\frac{1}{2}} w^+$
      
       di mana $g(u)$ yang digunakan adalah sebagai berikut

$g(u) = \tanh(a_1 u)$


$g^{'}(u) = a_1 (1-\tanh(a_1 u)^2)$

   c) cek konvergensi

$1 - \min{(\textrm{diag}(|w_n w_{n+1}^T|))} < \textrm{eps}$

       jika belum konvergen, kembali ke step b, di mana $w_n = w_{n+1}$

Keempat, estimasi nilai $A$
Untuk mendapatkan nilai $A$, kita harus mentransformasi balik nilai $w_{n+1}$ yang sebelumnya telah melalui tahap whitening. Proses ini bisa kita sebut juga sebagai dewhitening.

$A = {E^T}^{-1} \sqrt{D} E^{-1} w_{n+1}$


Kelima, estimasi $s$
Tahap terakhir ini sangat mudah

$ s = A^{-1} x $

Octave Code

Simpan code berikut dengan nama "fast_ica.m"

----------------------------------------------------------------------------------------------------------------------------------

function [s A] = fast_ica(x)

% Centering
xc = x - repmat(mean(x')',1,size(x,2));

% whitening
[E D] = eig(cov(xc',1));
xw = E * inv(sqrt (D)) * E' *xc;

% Algoritma Fast Ica
wn = rand(size(xc,1));
for n = 1:5000
    w_plus = xw * func_g(wn,xw) - diff_func_g(wn,xw) .* wn;
    wn1 = w_plus * real(inv(w_plus' * w_plus)^(1/2));
    if 1 - min(abs(diag(wn*wn1'))) < eps
        fprintf (stderr, ['W konvergen pd iterasi ke - ' num2str(n)  "\n"]);
        break
    else
        wn = wn1;
    end
end

% Estimasi nilai A
A = inv(E')*sqrt(D)*inv(E)*wn1;

% Estimasi s
s = inv(A)*x;

% ======================================
% subfunction

function g = func_g(wn,xw)
a1 = 1;
u = xw'*wn;
g = tanh(a1*u);
dg = (1 - (tanh(a1*u).^2))*a1;


function dg = diff_func_g(wn,xw)
a1 = 1;
u = xw'*wn;
dg = (1 - (tanh(a1*u).^2))*a1;
dg = ones(size(wn,1),1)*sum(dg);

----------------------------------------------------------------------------------------------------------------------------------


Contoh Program untuk menguji Algoritma ICA

Simpan code berikut di folder yang sama dengan "fast_ica.m", kemudian F5 (run)

close all
clear
clc

fs = 10000;
t = 0:1/fs:2;

f1 = 5;   % frekuensi sinyal 1
f2 = 20;  % frekuensi sinyal 2
f3 = 2;   % frekuensi sinyal 3

s1 = sin(2*pi*f1*t);  % sinyal 1
s2 = sin(2*pi*f2*t);  % sinyal 2
s3 = sin(2*pi*f3*t);  % sinyal 3
% tiga sinyal dalam satu matrik
s = [s1;s2;s3];

% plot sinyal asli
figure
for i = 1:3
  subplot(3,1,i)
  plot(t,s(i,:),'k','LineWidth',2)
end

% mixing matrik
A = rand(3);  % 3 menunjukkan tiga buah sinyal dicampur
% sinyal tercampur
x = A*s;

% plot sinyal tercampur
figure
for i = 1:3
  subplot(3,1,i)
  plot(t,x(i,:),'k','LineWidth',2)
end

% fast ICA
[s_estimasi A_estimasi] = fast_ica(x);
 
% plot sinyal estimasi
figure
for i = 1:3
  subplot(3,1,i)
  plot(t,s_estimasi(i,:),'k','LineWidth',2)
end

Sinyal-sinyal sebelum tercampur

Sinyal-sinyal setelah proses pencampuran

Sinyal-sinyal estimasi

Jika kalian malas menyalin code di atas, kalian bisa download langsung di link ini.
Selamat mencoba dan semoga bermanfaat.

Referensi
Aapo Hyvärinen and Erkki Oja, Independent Component Analysis: Algorithms and Applications, Neural Networks , 13(4-5): 411-430, 2000

5 comments:

  1. Knpa u = xw'*wn
    Bukannya u = wn'*xw

    ReplyDelete
    Replies
    1. Terima kasih atas pertanyaannya. Mudahah program ini sudah anda download dan dijalankan.

      Saya tidak seberapa ahli dalam penurunan rumus ICA, sekedar menterjemahkan dari persamaan ICA ke bahasa pemrograman octave. Semoga jawaban berikut ini bisa membantu.

      Referensi utama diambil dari pdf di link ini http://mlsp.cs.cmu.edu/courses/fall2012/lectures/ICA_Hyvarinen.pdf

      Saya pribadi agak bingung saat pertama kali menterjemahkannya ke Octave, dan bagian yang saya bingungkan adalah yang sedang anda tanyakan tersebut.

      dalam program yang saya buat di atas, data yang dicampur terdiri dari tiga sinyal sehingga "x" saya definisikan berukuran 3 x 20001 dan otomatis data setelah whitening "xw" juga berukuran 3 x 20001.

      pada persamaan "w_plus = xw * func_g(wn,xw) - ........"

      dimana "w" adalah matrik pembebanan (keterangan lengkap ada di pdf di link yang saya berikan) berukuran 3x3. "w" ini adalah parameter penting yang ingin kita estimasi agar nantinya kita bisa mendapatkan mixing matrix.

      karena

      w_plus (3 x 3) = xw (3 x 20001) * func_g ( .... x ....)

      berdasar persamaan di atas, saya simpulkan func_g harus berukuran 20001 x 3.
      sehingga "u" yang seharusnya "u = wn'*xw" (3x3 * 3x20001) dibalik menjadi "u = xw'*wn" (20001x3 * 3x3) agar func_g berukuran 20001x3. Apakah hal ini BENAR? saya tidak berani mengatakan jawaban ini 100% benar. Namun jika persamaan di dalam "u" tersebut saya balik, saya bisa mendapatkan hasil simulasi yang benar (beberapa kali saya mencoba dengan cara yang lain, namun masih gagal).

      Mudahan jawaban ini bisa membantu dan monggo jika anda mempunyai jawaban yang tepat beserta kode-nya agar dishare di sini,. Nanti postingan di atas akan saya edit. Terima kasih

      Delete
  2. Minta referensi buku tentang ICA ini dong kk

    ReplyDelete
  3. kalo saya punya suara orang lagi ngobrol,,,
    jadi dengan ICA ini bisa dipisahkan setiap gelombang suara nya..
    gitu gan ?

    ReplyDelete
  4. kalau rumus ica dengan metode maksimum adan diferensial gimana?

    ReplyDelete