import numpy as np
import matplotlib.pyplot as plt
plt.style.use('ggplot')
from scipy import signal
from joblib import Parallel, delayed
from spkit import ICA
from spkit.data import load_data
X,ch_names = load_data.eegSample()
band =[1,63.5]
order,fs = 5,128
b,a = signal.butter(order,np.array(band)/(0.5*fs),btype='bandpass')
Xf = np.array(Parallel(n_jobs=-1)(delayed(signal.lfilter)(b,a,X[:,i]) for i in range(X.shape[1]))).T
t = np.arange(X.shape[0])/128.0
plt.figure(figsize=(12,5))
plt.plot(t,Xf+np.arange(-7,7)*200)
plt.xlim([t[0],t[-1]])
plt.xlabel('time(sec)')
plt.yticks(np.arange(-7,7)*200,ch_names)
plt.show()
methods = ('fastica', 'infomax', 'extended-infomax', 'picard')
icap = ['ICA'+str(i) for i in range(1,15)]
x = Xf[128*10:128*12,:]
t = np.arange(x.shape[0])/128.0
myICA = ICA(n_components=14,method='fastica')
myICA.fit(x.T)
s1 = myICA.transform(x.T)
myICA = ICA(n_components=14,method='infomax')
myICA.fit(x.T)
s2 = myICA.transform(x.T)
myICA = ICA(n_components=14,method='picard')
myICA.fit(x.T)
s3 = myICA.transform(x.T)
myICA = ICA(n_components=14,method='extended-infomax')
myICA.fit(x.T)
s4 = myICA.transform(x.T)
plt.figure(figsize=(15,15))
plt.subplot(321)
plt.plot(t,x+np.arange(-7,7)*200)
plt.xlim([t[0],t[-1]])
plt.yticks(np.arange(-7,7)*200,ch_names)
plt.title('X : EEG Data')
plt.subplot(322)
plt.plot(t,s1.T+np.arange(-7,7)*700)
plt.xlim([t[0],t[-1]])
plt.yticks(np.arange(-7,7)*700,icap)
plt.title('FastICA')
plt.subplot(323)
plt.plot(t,s2.T+np.arange(-7,7)*700)
plt.xlim([t[0],t[-1]])
plt.yticks(np.arange(-7,7)*700,icap)
plt.title('Infomax')
plt.subplot(324)
plt.plot(t,s3.T+np.arange(-7,7)*700)
plt.xlim([t[0],t[-1]])
plt.yticks(np.arange(-7,7)*700,icap)
plt.title('Picard')
plt.subplot(325)
plt.plot(t,s4.T+np.arange(-7,7)*700)
plt.xlim([t[0],t[-1]])
plt.yticks(np.arange(-7,7)*700,icap)
plt.title('Extended-Infomax')
plt.show()
Computing Extended Infomax ICA
where $\mu$ is mean, computed before applying PCA
mu = myICA.pca_mean_[:, None]
W = myICA.get_sMatrix()
A = myICA.get_tMatrix()
s1 = np.dot(A,(x.T-mu))
x1 = np.dot(W,s4)+mu
plt.figure(figsize=(15,5))
plt.subplot(121)
plt.plot(x1.T+np.arange(-7,7)*400)
plt.yticks(np.arange(-7,7)*400,ch_names)
plt.title('Reconstructed X')
plt.subplot(122)
plt.plot(s1.T+np.arange(-7,7)*400)
plt.title('Computed ICA Ex-InfoMax')
plt.yticks(np.arange(-7,7)*400,icap)
plt.show()
print('Error X',np.sum(x1-x.T))
print('Error S',np.sum(s1-s4))
Error X 8.51877856450356e-12 Error S -2.443648061678516e-11