DC offset removal filter (in Matlab)
DC offset removal filter
Sometimes slow drifts or DC drifts mask low frequencies that could be of interest, especially when the analysis included fourier-based methods. For example the Contingent Negative Variation (CNV) signal is a very low frequency signal that should be preserved while avoiding to include in the analysis the slow drifts inherent sometimes with EEG recordings. These filters were made for a CNV EEG analysis. Sometimes, in other kinds of recordings an initial drift from baseline exists at the beginning of recordings. For those cases, a specialized DC offset removal filter is necessary for not cutting out parts of data that could be of interest.
Two filters, two solutions are proposed below:
(the examples assume your data are already loaded in EEGLAB)
Solution a: Using an IIR low-pass filter
General idea: Application of a low pass FIR filter (Roll off 0.001 to 0.05 Hz, butter, (IIR), 6 dB attenuation in the stop band), then removal from the signal the filtered lowpass signal
Filtered_low_pass=filterLP(signal)
Final_signal=signal –filterLP(signal).
This is depicted in the following figure for a representative channel.
Figure 1: channel 1 of dataset 101b, raw dataset (blue), the Filtered_low_pass(red) and the final signal Filtered_signal (green). A sinusoid of 0.01 Hz is also plotted for comparison.
If we apply this procedure to 8 selected channels we get the following picture:
Figure 2: filtered (blue) and raw (green).
The lowpass filter characteristics are
h=fdesign.lowpass('Fp,Fst,Ap,Ast',0.001,0.05,1,6,250);
d=design(h,'butter');
fvtool(d)
The program:
DC_offset_removal_21_10_2011_a.m
%% DC offset removal. 21.10.2011 Maria L. Stavrinou
%a)solution 1
%% 1. load the EEG.data and save it as another set
data= EEG.data; % nchan x timeduration_all
%% 2 Low pass filter
h=fdesign.lowpass('Fp,Fst,Ap,Ast',0.001,0.05,1,6,250);
d=design(h,'butter');
fvtool(d)
% check on a single channel
ch1=double(data(1,:));
y=filter(d, ch1);%y=filtfilt(d.Numerator,1,ch1);
figure; plot(ch1); hold on; plot(y, 'r')
hold on; plot((ch1-y), 'g')
%% remove DC from all channels
nchan=size(data,1);
for k=1:nchan
temp=double(data(k,:));
%temp_filt=filtfilt(d.Numerator,1,temp);
temp_filt=filter(d,temp);
data_filt(k,:)=-data(k,:)+temp_filt;
%data_filt(k,:)=(data_filt(k,:)-mean(data_filt(k,:)))
clear temp temp_filt
end
%% 4 check on one file & decide how to cut the file.
% %% see an example file
ch12=data_filt([5 6 7 12 62 72 106 129],:);
figure(100);
nn=size(ch12,1);
for kk=1:nn
a=squeeze(ch12(kk,:));
plot(a);hold on;
clear a
end
hold on;
raw_ch12=data([5 6 7 12 62 72 106 129],:);
for kk=1:nn
a=squeeze(raw_ch12(kk,:));
plot(a, 'g');hold on;legen1=legend('green=raw blue=filtered')
clear a
end
disp('Write down how to cut from the start in datapoints')
%% 5 reload the filtered back to the EEGLAB structure
EEG.data=data_filt;
Solution b: Using an FIR filter without phase shift (time delays)
General idea: Application of a low pass FIR filter (Roll off 0.001 to 0.05 Hz, Equiripple, (FIR), 8 dB attenuation in the stop band), then removal from the signal the filtered lowpass signal
Filtered_low_pass=filterLP(signal)
Final_signal=signal –filterLP(signal).
Filtered_signal=Filtered_signal-mean(Filtered_signal); (for every channel, its own mean is subtracted)
This is depicted in the following figure for a representative channel.
Figure 1: channel 1 of dataset 101b, raw dataset (blue), the Filtered_low_pass(red) and the final signal Filtered_signal (green).
If we apply this procedure to 8 selected channels we get the following picture:
Figure 2: filtered (blue) and raw (green).
The lowpass filter characteristics are
h=fdesign.lowpass('Fp,Fst,Ap,Ast',0.001,0.05,1,8,250);
d=design(h,'equiripple');
fvtool(d)
The program:
DC_offset_removal_21_10_2011_b.m
%% DC offset removal. 21.10.2011 Maria L. Stavrinou
%b)solution 2 %%
%% load eeg data
data= EEG.data; % nchan x timeduration_all
%% 2 Low pass filter
h=fdesign.lowpass('Fp,Fst,Ap,Ast',0.001,0.05,1,8,250);
d=design(h,'equiripple');
fvtool(d)
% check on a single channel
ch1=double(data(1,:));
y=filtfilt(d.Numerator,1,ch1);
figure; plot(ch1); hold on; plot(y, 'r')
hold on; plot((ch1-y), 'g')
%% 3 remove DC from all channels
nchan=size(data,1);
for k=1:nchan
temp=double(data(k,:));
temp_filt=filtfilt(d.Numerator,1,temp);
%temp_filt=filter(d,temp);
data_filt(k,:)=data(k,:)-temp_filt;
data_filt(k,:)=(data_filt(k,:)-mean(squeeze(data_filt(k,:))));
clear temp temp_filt
end
%% 4 check on channels of interest
ch12=data_filt([5 6 7 12 62 72 106 129],:);
figure(100);
nn=size(ch12,1);
for kk=1:nn
a=squeeze(ch12(kk,:));
plot(a);hold on;
clear a
end
hold on;
raw_ch12=data([5 6 7 12 62 72 106 129],:);
for kk=1:nn
a=squeeze(raw_ch12(kk,:));
plot(a, 'g');hold on;legen1=legend('green=raw blue=filtered')
clear a
end
clear nn
disp('Write down how to cut from the start in datapoints')
%% 5 reload the filtered back to the EEGLAB structure
EEG.data=data_filt;