Lecture Notes on class

Ex01_testPCA.m

@@ -1,5 +1,4 @@
-% Principal component analysis
-%
+% Boyang Bao (bbao5@emory.edu)
 % BMI500 Course
 % Lecture:  An Introduction to Blind Source Separation and Independent Component Analysis
 %           By: R. Sameni
@@ -36,3 +35,55 @@
 N = size(x, 1); % The number of channels
 T = size(x, 2); % The number of samples per channel
 
+PlotECG(x, 4, 'b', fs, 'Raw data channels');
+
+x_demeaned = x - mean(x, 2) * ones(1, size(x, 2));
+
+Cx = cov(x_demeaned');
+
+[V, D] = eig(Cx, "vector"); % in vector form, eigenvalues in reversed form
+
+figure
+subplot(121)
+plot(D(end:-1:1));
+grid
+xlabel('Index');
+ylabel('Eigenvalue');
+title('Eigenvalues in linear scale');
+subplot(122)
+plot(10*log10(D(end:-1:1)/D(end)));  % log-normalize, see the decimal
+grid
+xlabel('Index');
+ylabel('Eigenvalue ratios in dB');
+title('Normalized eigenvalues in log scale');
+
+
+x_var = var(x_demeaned, [], 2);
+x_var2 = diag(Cx);  % the variance
+
+y = V' * x_demeaned;
+Cy = cov(y');
+y_var = diag(Cy);
+
+x_total_energy = sum(x_var);
+Cx_trace = trace(Cx);
+eigenvalue_sum = sum(D);
+Cy_trace = trace(Cy);
+
+x_partial_energy = 100.0 * cumsum(D(end:-1:1))./x_total_energy;
+
+th = 99.9;
+N_eigs_to_keep = find(x_partial_energy <= th, 1, 'last');
+
+x_compressed = V(:, N-N_eigs_to_keep+1:N) * y(N-N_eigs_to_keep+1:N, :);
+
+t = (0: T-1)/fs;
+for ch = 1:N
+    figure
+    hold on
+    plot(t, x(ch, :));
+    plot(t, x_compressed(ch, :));
+    legend(['channel' num2str(ch)], 'compressed');
+    grid
+end
+

Ex02_testEigenAnalysisPowerMethod.m

@@ -1,5 +1,4 @@
-% The power method for eigenvalue decomposition
-%
+% Boyang Bao (bbao5@emory.edu)
 % BMI500 Course
 % Lecture:  An Introduction to Blind Source Separation and Independent Component Analysis
 %           By: R. Sameni
@@ -24,3 +23,19 @@
 % Cx = x * x';
 Cx = cov(x');
 
+[V,D] = eig(Cx)
+
+Itr = 100;
+
+v0 = rand(N,1);
+v1 = EigenAnalysisPowerMethod(Cx, v0, Itr);
+scale1 = (Cx*v1) ./v1;
+lambda1 = mean(scale1);
+
+C = Cx - lambda1 * (v1 * v1');
+v2 = EigenAnalysisPowerMethod(C, v0, Itr);
+scale2 = (Cx*v2)./v2;
+lambda2 = mean(scale2);
+
+C = C - mean(lambda2) * (v2 * v2');
+

Ex02_testEigenAnalysisPowerMethod.m~

@@ -0,0 +1,32 @@
+% Boyang Bao (bbao5@emory.edu)
+% BMI500 Course
+% Lecture:  An Introduction to Blind Source Separation and Independent Component Analysis
+%           By: R. Sameni
+%           Department of Biomedical Informatics, Emory University, Atlanta, GA, USA
+%           Fall 2020
+%
+% Dependency: The open-source electrophysiological toolbox (OSET):
+%       https://github.com/alphanumericslab/OSET.git
+%   OR
+%       https://gitlab.com/rsameni/OSET.git
+%
+
+close all;
+clear
+clc;
+
+% Build a random signal
+N = 3;
+T = 1000;
+a = randn(1, N);
+x = diag(a) * randn(N, T);
+% Cx = x * x';
+Cx = cov(x');
+
+[V,D] = eig(Cx)
+
+Itr = 100;
+
+v0 = rand(N,1);
+v1 = EigenAnalysisPowerMethod(Cx, v0, Itr)
+

Ex03_testICAmethods.m

@@ -1,5 +1,4 @@
-% Independent component analysis using classical methods
-%
+% Boyang Bao (bbao5@emory.edu)
 % BMI500 Course
 % Lecture:  An Introduction to Blind Source Separation and Independent Component Analysis
 %           By: R. Sameni
@@ -16,7 +15,7 @@
 clear
 close all
 
-example = 3;
+example = 1; %change here with 1 or 3
 switch example
     case 1 % A sample EEG from the OSET package
         load EEGdata textdata data % Load a sample EEG signal
@@ -46,3 +45,33 @@
 N = size(x, 1); % The number of channels
 T = size(x, 2); % The number of samples per channel
 
+PlotECG(x, 4, 'b', fs, 'Raw data channels');
+
+approach = "symm";  % defl: extract component one by one
+g = "tanh";
+lastEigfastica = N; % PCA stage
+numOfIC = N; % ICA stage
+interactivePCA = 'off';
+[s_fastica, A_fastica, W_fastica] = fastica(x, 'approach', approach, ...
+        'g', g, 'lastEig', lastEigfastica, 'numOfIC', numOfIC, ...
+        'interactivePCA', interactivePCA);
+
+Cs = cov(s_fastica');
+
+% JADE
+lastEigJADE = N;
+W_JADE = jadeR(x, lastEigJADE);
+s_jade = W_JADE * x;
+
+% SOBI
+lastEigSOBI = N;
+num_cov_matrices = 100;
+[W_SOBI, s_sobi] = sobi(x, lastEigSOBI, num_cov_matrices);
+
+
+PlotECG(s_fastica, 4, 'r', fs, 'Fastica');
+
+PlotECG(s_jade, 4, 'k', fs, 'JADE');
+
+PlotECG(s_sobi, 4, 'b', fs, 'SOBI');
+

Ex04_testEOGArtifactRemoval.m

@@ -1,5 +1,4 @@
-% Removing EOG artifacts from EEG signals
-%
+% Boyang Bao (bbao5@emory.edu)
 % BMI500 Course
 % Lecture:  An Introduction to Blind Source Separation and Independent Component Analysis
 %           By: R. Sameni

Ex05_testFetalECGExtraction.m

@@ -1,5 +1,4 @@
-% Extracting fetal ECG signals using various ICA algorithms
-%
+% Boyang Bao (bbao5@emory.edu)
 % BMI500 Course
 % Lecture:  An Introduction to Blind Source Separation and Independent Component Analysis
 %           By: R. Sameni