Independent Component Analysis (ICA) menggunakan Python

Tulisan kali ini membahas hal yang sama persis dengan tulisan sebelumnya yakni tentang ICA. Letak perbedaannya adalah pada tulisan sebelumnya (bisa dilihat di link ini) saya menggunakan Octave. Sedangkan pada kesempatan kali ini saya mencoba menggunakan Python. Jadi, isi tulisan kali ini hanya berisi pengaplikasian Python code pada ICA saja tanpa disertai dengan penjelasan terkait ICA itu sendiri. Tidak ada alasan khusus kenapa saya menuliskan ICA menggunakan Python meski sangat banyak sekali code tentang ICA di internet yang lebih cepat dan efisien. Ini hanyalah sekedar tulisan seorang "beginner" yang sedang belajar tentang Python untuk Scientific Computing. Oke let's start!

Simpan code berikut dengan nama "ica.py". Code ini kita anggap sebagai modul ICA dimana modul ini nantinya akan dipanggil.

import numpy as np
import numpy.linalg as LA
from scipy.linalg import sqrtm

def fG(wn,xw):
  a1 = 1
  u = np.dot(xw.transpose(),wn)
  g = np.tanh(a1*u)
  dg = (1-g**2)*a1
  dg = np.dot(np.ones((wn.shape[0],1)),np.sum(dg,axis=0).reshape(1,dg.shape[1]))
  return g, dg

def fast_ica(x):
  # centering
  xc = x - np.tile(np.mean(x,axis=1).reshape(x.shape[0],1),(1,x.shape[1]))
 
  # whitening
  D,E = LA.eig(np.cov(xc,bias=1))
  D = np.diag(D)
  xw = np.dot(np.dot(E,LA.inv(np.sqrt(D))),np.dot(E.transpose(),xc))
 
  # Algoritma Fast ICA
  wn = np.random.rand(xc.shape[0],xc.shape[0])
  for i in range(5000):
    w_plus = np.dot(xw,fG(wn,xw)[0]) - fG(wn,xw)[1]*wn
    wn1 = np.dot(w_plus,np.real(LA.inv(sqrtm(np.dot(w_plus.transpose(),w_plus)))))
    if 1 - np.min(np.abs(np.diag(np.dot(wn,np.transpose(wn1))))) < np.spacing(1):
      print 'W konvergen pada iterasi ke - %i' % (i+1)
      break
    else:
      wn = wn1
   
  Aest = np.dot(np.dot(LA.inv(E.transpose()),np.sqrt(D)),np.dot(LA.inv(E),wn1))
  # Estimasi s
  Sest = np.dot(LA.inv(Aest),x)
  return Sest, Aest

Kemudian, simpan lagi code di bawah ini dengan nama "main.py"

import numpy as np
import matplotlib.pyplot as plt
import time
import ica

t0 = time.time()

fs = 10**4
t = np.arange(0,2,1./fs)

f1 = 5  # frekuensi sinyal 1
f2 = 20 # frekuensi sinyal 2
f3 = 2  # frekuensi sinyal 3

s1 = np.sin(2*np.pi*f1*t) # sinyal 1
s2 = np.sin(2*np.pi*f2*t) # sinyal 2
s3 = np.sin(2*np.pi*f3*t) # sinyal 3

# tiga sinyal dalam satu matrik
s = np.array([s1,s2,s3])

# plot sinyal asli
plt.figure()
for i in range(3):
  plt.subplot(3,1,i+1)
  plt.plot(t,s[i])
 
# Mixing matrix
A = np.random.rand(3,3)
# sinyal tercampur
x = np.dot(A,s)

# plot sinyal tercampur
plt.figure()
for i in range(3):
  plt.subplot(3,1,i+1)
  plt.plot(t,x[i])

# Fast ICA
Sest, Aest = ica.fast_ica(x)

# plot sinyal estimasi
plt.figure()
for i in range(3):
  plt.subplot(3,1,i+1)
  plt.plot(t,Sest[i])

t1 = time.time()
plt.show()

print 'Elapsed time %.5f seconds.' % (t1-t0)

Langkah terakhir, jalankan "main.py" di komputer anda.

Di Ubuntu, saya lebih suka menggunakan Ipython. Jadi biasanya cukup ketik

In [1]: %run main.py
W konvergen pada iterasi ke - 6
Elapsed time 0.98528 seconds.

Klo melalui terminal,

user@my-pc$ python main.py 

W konvergen pada iterasi ke - 7
Elapsed time 0.49387 seconds.
Contoh gambar yang muncul adalah seperti di bawah ini

Sinyal sinyal sebelum dicampur

Sinyal sinyal hasil pencampuran

Sinyal hasil permisahan menggunakan ICA

Untuk code-nya bisa didownload di link ini. 

Selamat mencoba