%%
clc;
clear
close all;
load B02T.mat;
fs=250;
% t = 0.004:1/fs:1;
eeg1 = data{1,1}.X;
channel_1= eeg1(:,1);
ch_1=channel_1(1:2000);
figure;
subplot(2,1,1);
plot(ch_1);
xlabel('Time (s)');
ylabel('Amplitude');
legend('ch1');
title("Plot channel "+1+ " for Structure 1X1");
hold on;
grid on;
eeg_signal = eeg1(1:2000); % first 2000 samples
subplot(2,1,2);
plot(eeg_signal,'black');
hold on;
grid on;
legend('EEG Signal');
title('Raw EEG Signal');
%%
[IMF, residual, info] = emd(eeg_signal);
imf_count = max(info.NumIMF);
figure;
for i=1:imf_count
subplot(6,2,i)
plot(IMF(:,i))
title("IMF"+i);
end
subplot(6,2,imf_count+1)
plot(residual)
title("Residue");
%%
thresholds = thselect(IMF, 'rigrsure'); % determine threshold for each IMF using universal thresholding
for i = 1:size(IMF, 2)
IMF(:, i) = wthresh(IMF(:, i), 's', thresholds(i));
end
% Reconstruct denoised signal from thresholded IMFs
denoised_signal = sum(IMF, 2) + residual;
clean_EEG = denoised_signal(1:2000);
% Plot results
figure;
subplot(2,1,1);
plot(eeg_signal);
title('Original signal');
subplot(2,1,2);
plot(clean_EEG);
legend( 'Denoised signal');
title('Denoised signal');
%%
original_signal = eeg_signal;
reconstructed_signal =clean_EEG;
% Calculate SNR
noise = original_signal - reconstructed_signal;
SNR = 20 * log10(norm(original_signal) / norm(noise));
% Calculate MSE
% MSE = mean((norm(original_signal) - norm(reconstructed_signal)).^2);
MSE = immse(norm(original_signal), norm(reconstructed_signal));
disp(['SNR: ', num2str(SNR), ' dB']);
disp(['MSE: ', num2str(MSE)]);
