…Trail of papers had a recent post on standard analyses in neuroscience, which reminded me that I’ve meaning to post about signal preprocessing for a while now.

If there’s too much line noise (electrical noise at multiples of 60Hz in North America and 50Hz in Europe) in a signal, this will render the signal unusable. So it’s important to try and eliminate line noise through proper grounding at the source. Yet some line noise is more or less inevitable, and you’ll want to diminish this noise through digital filtering. This is useful, for example, as a preprocessing step for LFPs, or for cleaning up wideband signals prior to spike detection.

People sometimes use notch filters for this purpose, but this can remove too much signal, such that there’s a gap in the frequency domain around 60Hz where there’s no power at all. A better approach IMHO is to fit a function to the signal power around 60Hz and invert it to obtain a filter.

Here’s a function which does this. It uses a flexible exponential family of functions of the form for the fits. It works in little chunks of data (by default 60s) so it’s able to track nonstationarities and doesn’t take a boatload of memory. Results of a fit to peak around 60 and 180Hz are shown at the top.

function [newdata] = chunkwiseDeline(data,sr,freqs,freqrange,chunksize,showOutput) %function [newdata] = chunkwiseDeline(data,sr,freqs,chunksize) %Removes line noise from a signal in small chunks of data. %data - an n x 1 vector containing continuous data %sr - scalar, the sampling rate of the signal %freqs - a vector of frequencies around which to remove line noise, % for example [60,180] %freqrange (optional) - a scalar specifying the range of frequencies % to use for fitting a parametric function for % line noise (default: 2Hz) %chunksize (optional) - the size of a chunk in seconds (default : 60) %showOutput (opt) - if true, show plots of fits (default: false) %Algo: for every block of data, take abs(fft(dat)). Then extract the % vector corresponding to [freq(1)-freqrange, freq(1)+freqrange]. % Fit a function thefun = @(p,x) p(1)*exp(-p(2)*abs(x-p(4)).^(p(5)))+p(3); % to this range of freqs. Invert the function to obtain a delining % filter and apply it. Repeat for other frequencies. Repeat for % all blocks. Reassemble blocks. data = data(:); if nargin < 4 || isempty(freqrange) freqrange = .5; end if nargin < 5 || isempty(chunksize) chunksize = 60; end if nargin < 6 || isempty(showOutput) showOutput = false; end chunksize = chunksize*sr; %zero pad origlen = length(data); data = [data(chunksize/2:-1:1);data]; nchunks = ceil((length(data)+chunksize/2)/chunksize); padsize = nchunks*chunksize - length(data); data = [data;data(end:-1:end-padsize+1)]; newdata = zeros(size(data)); nchunks = nchunks*2-1; thewin = (0:chunksize/2-1)'/(chunksize/2-1); thewin = [thewin;thewin(end:-1:1)]; %Deline each chunk and reassemble pss = nan(5,nchunks); for ii = 1:nchunks thechunk = data((ii-1)*chunksize/2+(1:chunksize)); [delinedchunk,pss] = delineChunk(thechunk,sr,showOutput,freqs,freqrange,pss); newdata((ii-1)*chunksize/2+(1:chunksize)) = ... newdata((ii-1)*chunksize/2+(1:chunksize)) + ... delinedchunk.*thewin; end newdata = newdata(chunksize/2 + (1:origlen)); end %Deline a single chunk of data function [y pss] = delineChunk(dat,sr,showoutput,freqs,freqrange,pss) ae = []; if mod(length(dat),2) == 1 ae = dat(end); dat = dat(1:end-1); end fftdat = fft(double(dat)); a = abs(fftdat); %Now remove line noise from datlo to obtain y % %Eliminate line noise at target frequencies thefilt = ones(size(a)); winlen = round(length(dat)/sr*freqrange); %Fit a curve to this chunk of frequencies opts = optimset('Display','Off','Jacobian','on','Algorithm','levenberg-marquardt'); n = 1; for tgtr = freqs peak = tgtr/sr*length(dat); rg = round(((peak-winlen):(peak+winlen)))'; datrg = a(rg); x = (-winlen:winlen)'/winlen*freqrange; %Only adjust a few parameters at a time %convergence is better this way %everything but the exponent %Find the peak datrgsm = conv(datrg,ones(21,1),'same'); [~,peakloc] = max(datrgsm); %Set the initial parameters x0 = [max(datrg)-median(datrg),1/.2^2,median(datrg),(peakloc-1-winlen)/winlen*freqrange,1]'; if ~any(isnan(pss(:,n))) x0([2,4,5]) = pss([2,4,5],n); end [ps] = lsqcurvefit(@(x,y) thefun([x;x0(5)],y,[1;1;1;1;0]),x0(1:4),x,datrg,[],[],opts); xd = ps(4); %Everything but the center [ps] = lsqcurvefit(@(x,y) thefun([x(1:3);xd;x(4)],y,[1;1;1;0;1]),[ps(1:3);x0(5)],x,datrg,[],[],opts); %Everything [ps] = lsqcurvefit(@(x,y) thefun(x,y,[1;1;1;1;1]),[ps(1:3);xd;ps(4)],x,datrg,[],[],opts); pss(:,n) = ps; %Good, now adjust the filter in this range accordingly thefilt(rg) = ps(3)./thefun(ps,x); b = thefilt(rg); thefilt(end-rg+2) = b; if showoutput subplot(length(freqs),1,n); plot(x+tgtr,datrg,x+tgtr,thefun(ps,x)); title(sprintf('%3.1f Hz',tgtr)); drawnow; [ps(2),ps(4),ps(5)] end n=n+1; end % y is datlo with line noise removed a = fftdat.*thefilt; y = [real(ifft(a));ae]; end function [y,J] = thefun(p,x,mask) E = abs(x-p(4)).^(p(5)); M = exp(-p(2)*E); y = p(1)*M+p(3); if nargout > 1 J = [ M,... -p(1)*E.*M,... ones(size(x)),... p(1)*p(2)*p(5)*sign(x-p(4)).*abs(x-p(4)).^(p(5)-1).*M,... -p(1)*p(2)*p(5)*E.*log(abs(x-p(4))+1e-6).*M]; J = J(:,mask==1); end end

First of all thank you for this post! It’s really helpfull!

I know this post is pretty old now and you might not read this comment anymore, but I’m not sure about your Jacobian.

My math skills are a little bit out of training, but I found:

-p(1)*p(2)*E.*log(abs(x-p(4))+1e-6).*M

as a derivative according to p(5), i.e. the same as you without the p(5) factor.

Great function! Thanks for sharing.

I found that it doesn’t work with non-integer sampling rates (an annoyance with our data acquisition hardware). Anyway, just in case anyone is interested, your function works with non-integer sampling rates with the following modification to the code:

% chunksize = chunksize*sr;

33 chunksize = round(chunksize*sr);

34 if mod(chunksize,2) == 1, chunksize = chunksize-1; end

I’m guessing there’s a minimum amount of spectral information you need in order to get a good fit at 60Hz? I have 3-second snippets (not from a continuous source, but collected every 10 seconds for an hour) and the noise was still present after running this code. But I really appreciate your contributions! By contrast, I am so sick of finding papers that describe their excellent processing/pre-processing strategies but don’t give any code example of how to implement it. It really defeats the hope of standardizing our analyses (and therefore, realistic comparisons).

Your code examples are easy to understand and very useful. If I (or others) use your code in published work, how would want it cited? Thanks again!

*to be clear, I am referring to 3 second LFP snippets. So my question is better worded as, ” Is there a minimum sampling duration needed for this type of filtering to work?”

I think it should still work though you might have to do some contortions, like averaging the power spectrum of all the snippets and then fitting the 60Hz peak to a Gaussian and then apply that filter to each snippet individually…