diff --git a/Ex01_testPCA.m b/Ex01_testPCA.m index 1181070..40b29d0 100644 --- a/Ex01_testPCA.m +++ b/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 + diff --git a/Ex02_testEigenAnalysisPowerMethod.m b/Ex02_testEigenAnalysisPowerMethod.m index f1d95be..58f7b4a 100644 --- a/Ex02_testEigenAnalysisPowerMethod.m +++ b/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'); + diff --git a/Ex02_testEigenAnalysisPowerMethod.m~ b/Ex02_testEigenAnalysisPowerMethod.m~ new file mode 100644 index 0000000..3f774ed --- /dev/null +++ b/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) + diff --git a/Ex03_testICAmethods.m b/Ex03_testICAmethods.m index af123be..6888005 100644 --- a/Ex03_testICAmethods.m +++ b/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'); + diff --git a/Ex04_testEOGArtifactRemoval.m b/Ex04_testEOGArtifactRemoval.m index b6c4122..1119825 100644 --- a/Ex04_testEOGArtifactRemoval.m +++ b/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 diff --git a/Ex05_testFetalECGExtraction.m b/Ex05_testFetalECGExtraction.m index 3ac2f31..88ee48e 100644 --- a/Ex05_testFetalECGExtraction.m +++ b/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