# Clean noise in neural data - MATLAB

I'm using Matlab to process neurophysiological data and i'm currently dealing with this sort of noise in my signal:

Close up

EDIT:

For clarification, the clean signal should look like this:

EDIT 2:

Clean

Noise contaminated

These spectra were taken from 10s segments as the whole block is considerably bigger. Both spectra are from segments that include the data in the first and third figures (with and without noise bits of data). The noise events are sparse within each data block.

This happens somewhat frequently across datasets and due to it, i'm having some trouble detecting events of interest (i.e. neural spike detection). As i'm trying to detect small amplitude spikes in the signal, it's critical to preserve as much of the original signal as possible.

Any suggestion on how to better approach this problem? What would be te best way to correct these shifts/spikes in the signal, keeping as much of the signal as possible?

Any suggestion/help will be immensely appreciated!

Cheers, M

• It is unclear to me (I do not have experience with neurophysiological data) what you consider as noise and what you want to detect. Can you clarify it? – m7913d Nov 27 '17 at 16:15
• @m7913d: the noise in this sample are the high amplitude peaks and upward shifts. The signal should (ideally) be in the range of the signal before the first big shift in the first figure. – Oiko Nov 27 '17 at 16:19
• Just added a figure of a clean segment of the same dataset. – Oiko Nov 27 '17 at 16:28
• To me it's obvious that all of the spikes are signals. The question is, which ones are you interested in? I think an identification of the neural signals is crucial. From my limited knowledge, I think that neural signals are short discharge pulses from the neurons. These seem to be fairly well categorized by experiments. Having that in hand, you could sort through the data looking for the characteristic spikes. Cern has done a lot of work with this sort of thing but from my work in cytometry, scatterplots can be useful for classification; lacking analytic criteria. – rrogers Nov 29 '17 at 14:02
• @rrogers: Thanks for the input! In fact what you see here is noise, probably due to some mechanical artefacts on the drive/electrode/cable interface. Real spikes (action potentials) are significantly smaller and do not produce the kind of baseline shift we see here. The problem is that this noise increases the amplitude standard deviation, which i use to detect real spikes (+/-100uV). Thus it is critical to be able to remove these artefacts. Filtering does not seem to solve this... How could I go about correcting these shifts while keeping the signal on the downward slope after the big peaks? – Oiko Nov 29 '17 at 18:32

Following, a way to remove the 3 exponential pulses.

1.- Rebuilding signal + interferer from graph provided

clear all;clc;close all

A1=A(:,:,1);
figure(1);imshow(A1)

A2=~imbinarize(A1);
figure(2);imshow(A2)
% 001


problem: bwmorph breaks causality and splits signal

A3=bwmorph(A2,'skel',Inf);
figure(3);imshow(A3);
% 002


[sz1 sz2]=size(A2)

for k=1:1:sz2
L1=A2(:,k);
n1=find(L1);

A2(:,k)=zeros(1,sz1);
if ~isempty(n1)
A2(floor(mean(n1)),k)=1;
end
end
figure(4);imshow(A2)
% 003


Now signal is correctly mapped.

Since the contaminated signal is available to the question originator, the original data should be inserted here:

s=0
for k=1:1:sz2
L1=A2(:,k)
s=[s find(L1)];
end
s(1)=[]


2.- where is mean(s)?

don't do just s=s-mean(s)

Left behind the y=0 axis tick on purpose

% 004

mean_s=287

s=-(s-mean_s)


3.- Rebuilding time reference

t=linspace(0,10,sz2)
dt=mean(diff(t))


4.- Plotting signal + interference

figure(5);plot(t,s);grid on
xlabel('t[seconds]');ylabel('s+i')  % signal+interferer


The VLF bending within the initial 5 seconds may come from conducted interference; contact, applying an ECG sensor on a mobile phone?

The base station or mobile seeking answer, starts with certain power and then as the packets are exchanged Base Station Mobile Terminal tell one another to increase or reduce power, usually the Base Station starts higher than needed to make sure reaching cell edge, and then reduces power until MS says level ok.

The negative sharp spikes show there may be something metallic reflecting the interfering signal.

You have provided a reference signal, what it should be like, but there's no time reference, just the amount of samples, without further information, it's not sure the reference signal duration is also 10 seconds, is it?

[pks,locs]=findpeaks(s,'MinPeakHeight',50,'MinPeakDistance',50)
s_peaks=s(locs)


5.- When do the 3 exponential decay interfering pulses start?

t_peaks=t(locs)


6.- What is the mean delay between those 3 pulses measured in time?

diff(t_peaks)


7.- What is the mean delay between those 3 pulses measured in amount samples?

diff(locs)


8.- What is the real amplitude of the pulses? approximately the 1st sample

ds_peaks=s(locs)-s(locs-10)


9.- what is the amount of samples per interfering pulse?

nT1=floor(mean(diff(locs)))


10.- What are the respective attenuation constant of the 3 exponential pulses?

a=-1./(dt*nT1)*log(1./ds_peaks)


11.- 1st interfering pulse start time measured directly on the graph

t0=.688


12.- What is the sample numeral of the start moment of the 1st interfering pulse?

nt0=[1:1:floor(mean(intersect(find(t>t0-.01),find(t<t0+.01))))]


13.- Values of each interfering pulse:

e1=ds_peaks(1)*exp(-a(1)*dt*([1:1:nT1]+t0))
e2=ds_peaks(2)*exp(-a(2)*dt*([1:1:nT1]+t0))
e3=ds_peaks(3)*exp(-a(3)*dt*([1:1:nT1]+t0))
e_interf=[zeros(1,numel(nt0)) e1 e2 e3 zeros(1,sz2-3*nT1-numel(nt0))]
hold all
plot(t,e_interf)


s2=s-e_interf
figure(6);plot(t,s2);grid on
xlabel('t[seconds]');ylabel('s+i-3pulses')  % signal+interferrer


If interested I can show how to use findpeaks to remove the sharp negative pulses.

• What is the relevance of the "Mobile Base Stations" paragraph to the original question? – A_A Jul 1 '18 at 6:32
• I am guessing what is producing such bursts, and the ones left behind? the negative pulses, may be same base station / phone / whatever using higher band, like the 3 removed pulses could be 2.4GHz and the negative ones same data but at a higher frequency therefore shorter. I just mentioned it in case the question originator would like to comment and add more details on the test set-up, does it makes sense? – John Bofarull Jul 3 '18 at 1:38

You can reduce the large decaying deviations by using a high pass filter. You might also reduce the strong spikes with a low pass filter, thus forming in combination, a band pass filter. The hard part is determining the low and high frequencies for the band pass and to some extent, the stop band attenuations, while at the same time making a judgement about what is an acceptable result.

If you added spectrums of both the noise corrupted and noise-free signals, more specific and more helpful advice could be offered. It would appear that your signals could be improved but more detail is needed to make specific suggestions.

response to edit 2

There is a lot of spectral overlap so there will be a limit to the degree of improvement. I would definitely first high pass filter everything above 3 Hz, to reduce the large exponential deviations. If that produces something promising, you could proceed with the harder large amplitude spikes. They are very detectable so you could ignore them or blank out those intervals but that may not be suitable to your application.

You might also try to reduce the large spikes with a soft limit in the time domain followed by low pass filtering the results below about 49Hz, or first try the low pass filter.

• Thanks for your help! I'll try to post the spectra as soon as possible. – Oiko Nov 28 '17 at 10:24
• I just added the figures of the spectra. They correspond to 10s of data encompassing the signals in my 1st and 3rd figures (above). The noise events (groups of 5 to 10 such shifts/peaks) are scattered randomly through the data blocks. Thanks in advance! – Oiko Nov 28 '17 at 13:47