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