# Removing line noise from LFPs, wideband signals

…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 $a\exp(-b|x-c|^d)+e$ 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

```

## 5 comments

1. Florian says:

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.

2. Daniel says:

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

3. Brian says:

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!

• Brian says:

*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?”

• xcorr says:

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…