bachelor_thesis/DataAnalysisToolbox/removehumnoise.m

function [z,y] = removehumnoise(x,sr)
%-------------------------------------------------------------------
% removal of recording electric noise with filters 
%-------------------------------------------------------------------
% Input
%  x  : data to denoise (one dimensional vector)
%  dt : time step of the data in [ms] 
% Output
%  z  : filtered data (hum noise and high-frequency noise removed)
%  y  : filtered data (hum noise removed)
%-------------------------------------------------------------------
% Technical notes 
%  In Matlab, you normally have to take two steps to apply a filter to 
%  data. In the first step, you construct a filter (fir1, cheby2, butter, 
%  etc.), whose properties can be visualized with the 'freqz' function. 
%  In the second step, you apply the filter to the data using the 
%  'filtfilt' function. In the code below, a bandstop filter is used 
%  to remove the hum (50Hz) noise (output y). And then a lowpass filter 
%  is additionally applied to remove high-frequency noise (output z). 
%-------------------------------------------------------------------
% Go Ashida (go.ashida@uni-oldenburg.de) Mar 2020 
% modified by Bjarne Schultze (Dec 2021)
%-------------------------------------------------------------------

%% flag for plotting the results 
plotflag = 1; % 1: plot, 0: no plot

%% time/frequency settings
fs = sr;
% 1000/dt; % [Hz] sampling rate
nyq = fs/2; % Nyquist frequency (used for constructing a filter)

%% defining filters 
% Here examples with the Chebyshev Type II filter are shown. 
% In Matlab, however, many more filters (such as FIR, IIR, 
% Butterworth, Bessel, etc.) are available. (see help) 

% filter settings
bf1 = 50; % [Hz] center frequency of bandstop filter
bw1 = 5; % [Hz] bandwidth of bandstop filter                | adapdet by BS
bdim = 2; % filter dimension
bdb = 10; % [dB] attenuation level (10 db?)
bf2 = 10; % [dB] cutoff frequency of lowpass filter 
 
% constructing cheby2 filter 
[b_bandstop, a_bandstop] = cheby2(bdim, bdb, [(bf1-bw1/2)/nyq, (bf1+bw1)/nyq], 'stop'); 
[b_lowpass, a_lowpass]  = cheby2(bdim, bdb, bf2/nyq, 'low');  % lowpass

% showing filter properties - for more detail, see help for freqz
%figure; freqz(b_bandstop, a_bandstop, 512, fs); 
%figure; freqz(b_lowpass, a_lowpass, 512, fs); 

%% applying filter
y = filtfilt(b_bandstop, a_bandstop, x); 
z = filtfilt(b_lowpass, a_lowpass, y); 

%% plotting
if(plotflag)
  tv = (0:length(x)-1)/sr; % time vector                    | adapted by BS
  figure; hold on;
  plot(tv,x,'r-'); 
  plot(tv,y,'b-');
  plot(tv,z,'g-');
  % labelling (added by BS)
  xlabel("time [s]"); ylabel("membrane potential [mV]");
  legend("original voltage trace", "hum noise filtered", ...
      "hum + high freq. noise filtered")
end