Processing and identifying noises in EEG signal
Show older comments
I'm working on a final project of processing the EEG signal from a 24-day-old infant. The data file is in .edf format, and I've extracted and plotted it as a time series. You can find the data file here: https://file.io/3fWPetZYclsH.
The project's requirements are:
- Analyze the data to identify signal and noises.
- Design filters to eliminate the DC component, Electricity noise, and others (if contaminated).
My first approach was to use the power spectral density (PSD), but this idea was rejected as my professor wanted us to use more fundamental algorithms like DFT of Hibert-Huang Transform (HHT). Another problem with my prior work was that my filtered data lost all of the data after a a specific frequency, which my professor attributed to using a 'second-order' filter. Below is my attempt on creating the filters. The two tasks i need help is:
- Use DFT and/or HHT to identify signal and noises. It is strictly required to identify where the noises are.
- Designing filters that effectively eliminate noise without losing essential data points.
I'd appreciate any guidance from anyone. Thank you in advance for your assistance!
def highpass_filter(data, sfreq, cutoff=0.5):
nyquist = 0.5 * sfreq
normal_cutoff = cutoff / nyquist
b, a = butter(1, normal_cutoff, btype='high', analog=False)
return filtfilt(b, a, data)
def notch_filter(data, sfreq, freq=50.0, quality=30.0):
nyquist = 0.5 * sfreq
notch_freq = freq / nyquist
b, a = iirnotch(notch_freq, quality)
return filtfilt(b, a, data)
def bandpass_filter(data, sfreq, lowcut=0.5, highcut=8.0):
nyquist = 0.5 * sfreq
low = lowcut / nyquist
high = highcut / nyquist
b, a = butter(1, [low, high], btype='band')
return filtfilt(b, a, data)
Accepted Answer
Star Strider
on 25 Jun 2024
0 votes
The data file doesn’t exist. At least it’s not available to me. It’s bettter to attach it here, anyway. If it has an excluded exttension, use the zip function to enclose it in a .zip file that you can upload here (providing it’s within the size limitts).
You can make things easier for yourself (since you obviously have the Signal Processing Toolbox) by using the highpass, bandstop, and bandpass functions. Use the 'ImpulseResponse','iir' name-value pair in each function call to force them to design a very efficientt elliptic filter. There are other name-value pair arguments to those functions that can allow you to modify the filter design, however I usually go with the default values. (An order 1 Butterworth filter isn’t going to do much actual filtering. Use the freqz function to view their Bode plots.)
Implementing the fft function is straightforward, and there are several examples in the documentation and in Answers. (I wrote several of them here, using essentially the same code, and most recently a function I created.)
Your approach is correct. The Fourier transform will demonstrate the frequency components of the signal. Frequency-selective filters will be useful for band-limited noise, including mains/powerline noise, however broadband noise will be difficult to deal with. There is no best solution for that, however wavelet denoising and the Savitzky-Golay filter are commonly used to reduce or eliminate broadband noise. The tradeoff is in reducing the noise without compromising the signal morphology. Badly designed filters that eliminate most of the noise can render the resulting filtered signal clinically useless.
5 Comments
Image Analyst
on 25 Jun 2024
What is his code? It does not look like MATLAB code.
David
on 25 Jun 2024
It took a bit of experimenting to get the edfread function to work (I don’tt have much experience with it), however after that, the rest is straightforward.
I’m not certain what you want to do, so I suggested some optionos here, including taking the Fourier transforms, filtering the signals, and calculating the band power in each lead (although you will need to check toe band power calculations to be certain they conform to the literature calculations).
A prototype analysis to get you started —
zf = dir('*.zip');
Uz = unzip(zf.name);
Lv1 = contains(Uz,'.edf');
EDF = Uz{Lv1};
Data = edfread(EDF)
VN = Data.Properties.VariableNames;
Fs = 1600/8;
tv = linspace(0, 1599, 1600).'/Fs;
Ofst = seconds(Data.('Record Time'));
tm = repmat(tv,1,size(Data,1)) + Ofst.';
t = tm(:);
% Data.Fp1{1}
Datac = table2array(Data);
DataMtx = cell2mat(cat(1,Datac));
Fp1v = DataMtx(:,1);
figure
plot(t, DataMtx)
grid
xlabel('Time (s)')
ylabel('Amplitude')
[FT_EEG,Fv] = FFT1(DataMtx,t);
figure
plot(Fv, abs(FT_EEG)*2)
grid
Ax = gca;
Ax.XScale = 'log';
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Fourier Transforms of EEG Data')
legend(VN, 'Location','best', 'NumColumns',2)
xlim([0 50])
figure
tiledlayout(5,4)
for k = 1:size(FT_EEG,2)
nexttile
semilogx(Fv, abs(FT_EEG(:,k))*2)
grid
title(VN{k})
ylim([0 30])
end
sgtitle('EEG Lead Fourier Transform Ensemble')
[Alpha,fdbp1] = bandpass(DataMtx, [8 13], Fs, 'ImpulseResponse','iir'); % Filter Leads Using Common EEG Band Definitions
[Beta,fdbp2] = bandpass(DataMtx, [13 30], Fs, 'ImpulseResponse','iir');
[Delta,fdbp3] = bandpass(DataMtx, [0.5 4], Fs, 'ImpulseResponse','iir');
[Theta,fdbp4] = bandpass(DataMtx, [4 7], Fs, 'ImpulseResponse','iir');
for k = 1:size(DataMtx,2) % Calculate Band Power In Each LEad
AlphaPwr(:,k) = trapz(t,Alpha(:,k).^2);
BetaPwr(:,k) = trapz(t,Beta(:,k).^2);
DeltaPwr(:,k) = trapz(t,Delta(:,k).^2);
ThetaPwr(:,k) = trapz(t,Theta(:,k).^2);
end
RowNames = {'Alpha Power','Beta Power','Delta Power','Theta Power'};
Band_Power_Mtx = cat(1,AlphaPwr,BetaPwr,DeltaPwr,ThetaPwr);
Band_PowerTable = array2table(Band_Power_Mtx, 'RowNames',RowNames, 'VariableNAmes',VN)
function [FTs1,Fv] = FFT1(s,t)
t = t(:);
L = numel(t);
if size(s,2) == L
s = s.';
end
Fs = 1/mean(diff(t));
Fn = Fs/2;
NFFT = 2^nextpow2(L);
FTs = fft((s - mean(s)) .* hann(L).*ones(1,size(s,2)), NFFT)/sum(hann(L));
Fv = Fs*(0:(NFFT/2))/NFFT;
% Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
FTs1 = FTs(Iv,:);
end
I will of course help with anything else you want to do with these signals.
.
David
on 25 Jun 2024
Star Strider
on 25 Jun 2024
Edited: Star Strider
on 25 Jun 2024
Thank you!
My apologies for the delay. Away for a bit.
There is not much high-frequency noise, and I didn’t see any mains/powerline noise at all. This indicates that a good reference electtrode was employed during the recording. There is a snr function I would recommend to calculate that. You can eliminate the D-C offset with a highpass filter. Experiment with the cutoff frrequency,m probably beginning at 3E-3 Hz, then re-calculate the fft and see the result. Another option, witthout using any filter, is simply to subtract the mean of each channel from that channel to eliminate the D-C offset, although my ‘FFT1’ function does that automatically, so there is actually no constant basellilne offset in the Fourier transformed signals, although there could be some baseline variation due to a ‘wandering’ non-constant baselline.
That would go something like this —
zf = dir('*.zip');
Uz = unzip(zf.name);
Lv1 = contains(Uz,'.edf');
EDF = Uz{Lv1};
Data = edfread(EDF)
VN = Data.Properties.VariableNames;
Fs = 1600/8;
tv = linspace(0, 1599, 1600).'/Fs;
Ofst = seconds(Data.('Record Time'));
tm = repmat(tv,1,size(Data,1)) + Ofst.';
t = tm(:);
% Data.Fp1{1}
Datac = table2array(Data);
DataMtx = cell2mat(cat(1,Datac));
DataMtxHPF = highpass(DataMtx, 3E-3, Fs, 'ImpulseResponse','iir');
figure
plot(t, DataMtx)
grid
xlabel('Time (s)')
ylabel('Amplitude')
[FT_EEG,Fv] = FFT1(DataMtxHPF,t);
figure
plot(Fv, abs(FT_EEG)*2)
grid
Ax = gca;
Ax.XScale = 'log';
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Fourier Transforms of EEG Data')
legend(VN, 'Location','best', 'NumColumns',2)
xlim([0 50])
figure
tiledlayout(5,4)
for k = 1:size(FT_EEG,2)
nexttile
semilogx(Fv, abs(FT_EEG(:,k))*2)
grid
title(VN{k})
ylim([0 30])
end
sgtitle('EEG Lead Fourier Transform Ensemble')
[Alpha,fdbp1] = bandpass(DataMtx, [8 13], Fs, 'ImpulseResponse','iir'); % Filter Leads Using Common EEG Band Definitions
[Beta,fdbp2] = bandpass(DataMtx, [13 30], Fs, 'ImpulseResponse','iir');
[Delta,fdbp3] = bandpass(DataMtx, [0.5 4], Fs, 'ImpulseResponse','iir');
[Theta,fdbp4] = bandpass(DataMtx, [4 7], Fs, 'ImpulseResponse','iir');
for k = 1:size(DataMtx,2) % Calculate Band Power In Each LEad
AlphaPwr(:,k) = trapz(t,Alpha(:,k).^2);
BetaPwr(:,k) = trapz(t,Beta(:,k).^2);
DeltaPwr(:,k) = trapz(t,Delta(:,k).^2);
ThetaPwr(:,k) = trapz(t,Theta(:,k).^2);
end
RowNames = {'Alpha Power','Beta Power','Delta Power','Theta Power'};
Band_Power_Mtx = cat(1,AlphaPwr,BetaPwr,DeltaPwr,ThetaPwr);
Band_PowerTable = array2table(Band_Power_Mtx, 'RowNames',RowNames, 'VariableNAmes',VN)
function [FTs1,Fv] = FFT1(s,t)
t = t(:);
L = numel(t);
if size(s,2) == L
s = s.';
end
Fs = 1/mean(diff(t));
Fn = Fs/2;
NFFT = 2^nextpow2(L);
FTs = fft((s - mean(s)) .* hann(L).*ones(1,size(s,2)), NFFT)/sum(hann(L));
Fv = Fs*(0:(NFFT/2))/NFFT;
% Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
FTs1 = FTs(Iv,:);
end
Make appropriate changes (including to the highpass stopband) to get the result you want.
EDIT — (25 Jun 2024 at 23:55)
Minor tweak, and typographcal error correction.
.
More Answers (1)
Diego Caro
on 26 Jun 2024
0 votes
If allowed by your instructors, you could use the EEGLAB toolbox, which is an environment bulit to analyze EEG signals. The GUI is very intuitive and easy to follow. To identify artifacts on your signal, you could run the ICA algorithm, which effectively classifies the orignal signal components into brain, muscle, eye, heart, line noise or other.
After performing the ICA algorith on an 8 electrode EEG signal, this is what comes out. In simple terms, it tells you the probability of each component to come from certain source.
As you can see, the second component has a 64.7% chance of being Line Noise. To identify wheter this is true or not, you could take a look at the extended analysis of this component.
From the power spectrum, it is obvious that this component is corrupted with Line Noise (because of the notorious spike at 60 Hz). After that, you could use the EEGLAB built-in functions, or desing your own to eliminate unwanted artifacts.
Regards,
Diego.
Categories
Find more on EEG/MEG/ECoG in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
