From cbefb9efc091f8b4c9f6fb7c0a7fda10b9e4b248 Mon Sep 17 00:00:00 2001
From: MarvinQuiet <mawenjing1993@126.com>
Date: Fri, 30 Oct 2020 11:59:56 -0400
Subject: [PATCH] lab during class

---
 Ex01_testPCA.m                      | 65 +++++++++++++++++++++++++++++
 Ex02_testEigenAnalysisPowerMethod.m |  8 ++++
 Ex03_testICAmethods.m               | 38 ++++++++++++++++-
 Ex04_testEOGArtifactRemoval.m       |  4 +-
 Ex05_testFetalECGExtraction.m       |  2 +-
 5 files changed, 114 insertions(+), 3 deletions(-)

diff --git a/Ex01_testPCA.m b/Ex01_testPCA.m
index 1181070..00078f5 100644
--- a/Ex01_testPCA.m
+++ b/Ex01_testPCA.m
@@ -1,4 +1,5 @@
 % Principal component analysis
+% Wenjing Ma (wenjing.ma@emory.edu)
 %
 % BMI500 Course
 % Lecture:  An Introduction to Blind Source Separation and Independent Component Analysis
@@ -36,3 +37,67 @@
 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');
+
+% remove mean
+x_demeaned = x - mean(x, 2) * ones(1, size(x, 2));
+
+% covariance matrix of the input
+Cx = cov(x_demeaned');
+
+% eigen value decomposition
+[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');
+
+% use eigenvalues to observe the energy (-> subset of eigenvalues)
+
+% check variance
+x_var = var(x_demeaned, [], 2);
+x_var2 = diag(Cx);  % the variance
+
+% Decorrelate the channels
+y = V' * x_demeaned;
+Cy = cov(y');
+y_var = diag(Cy);
+
+% check total engegy match
+x_total_energy = sum(x_var);
+Cx_trace = trace(Cx);
+eigenvalue_sum = sum(D);
+Cy_trace = trace(Cy);
+
+% did not scale, just rotating, so the variances do not change
+
+% partial energy
+x_partial_energy = 100.0 * cumsum(D(end:-1:1))./x_total_energy;
+
+% set a cut off 
+th = 99.9;
+N_eigs_to_keep = find(x_partial_energy <= th, 1, 'last');
+
+% find a compressed version of x
+x_compressed = V(:, N-N_eigs_to_keep+1:N) * y(N-N_eigs_to_keep+1:N, :);
+% eig(cov(x_compressed'))
+
+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
\ No newline at end of file
diff --git a/Ex02_testEigenAnalysisPowerMethod.m b/Ex02_testEigenAnalysisPowerMethod.m
index f1d95be..0df77eb 100644
--- a/Ex02_testEigenAnalysisPowerMethod.m
+++ b/Ex02_testEigenAnalysisPowerMethod.m
@@ -1,4 +1,5 @@
 % The power method for eigenvalue decomposition
+% Wenjing Ma (wenjing.ma@emory.edu)
 %
 % BMI500 Course
 % Lecture:  An Introduction to Blind Source Separation and Independent Component Analysis
@@ -24,3 +25,10 @@
 % Cx = x * x';
 Cx = cov(x');
 
+% HW, power method to calculate eigenvalue decomposition
+% eigs function, for high-dimensional data
+% the algorithm is in the slides (power method for EVD)
+
+
+
+
diff --git a/Ex03_testICAmethods.m b/Ex03_testICAmethods.m
index af123be..abb9dda 100644
--- a/Ex03_testICAmethods.m
+++ b/Ex03_testICAmethods.m
@@ -16,7 +16,7 @@
 clear
 close all
 
-example = 3;
+example = 1;
 switch example
     case 1 % A sample EEG from the OSET package
         load EEGdata textdata data % Load a sample EEG signal
@@ -46,3 +46,39 @@
 N = size(x, 1); % The number of channels
 T = size(x, 2); % The number of samples per channel
 
+% test-ICA procedure, fastica, JADE, SOBI
+%PlotECG(x, 4, 'b', fs, 'Raw data channels');
+
+% fastica
+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);
+    
+% 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);
+
+% Plot
+PlotECG(s_fastica, 4, 'r', fs, 'fastica');
+PlotECG(s_jade, 4, 'k', fs, 'JADE');
+PlotECG(s_sobi, 4, 'b', fs, 'SOBI');
+
+% need domain knowledge to tell the channel and whether to trust the result
+
+
+
+
+
+
+
diff --git a/Ex04_testEOGArtifactRemoval.m b/Ex04_testEOGArtifactRemoval.m
index b6c4122..c7a5052 100644
--- a/Ex04_testEOGArtifactRemoval.m
+++ b/Ex04_testEOGArtifactRemoval.m
@@ -1,5 +1,5 @@
 % Removing EOG artifacts from EEG signals
-%
+% Wenjing Ma (wenjing.ma@emory.edu)
 % BMI500 Course
 % Lecture:  An Introduction to Blind Source Separation and Independent Component Analysis
 %           By: R. Sameni
@@ -32,3 +32,5 @@
 T = size(x, 2); % The number of samples per channel
 t = (0 : T - 1)/fs;
 
+
+
diff --git a/Ex05_testFetalECGExtraction.m b/Ex05_testFetalECGExtraction.m
index 3ac2f31..4af763e 100644
--- a/Ex05_testFetalECGExtraction.m
+++ b/Ex05_testFetalECGExtraction.m
@@ -1,5 +1,5 @@
 % Extracting fetal ECG signals using various ICA algorithms
-%
+% Wenjing Ma (wenjing.ma@emory.edu)
 % BMI500 Course
 % Lecture:  An Introduction to Blind Source Separation and Independent Component Analysis
 %           By: R. Sameni