Clean noise in neural data - MATLAB

Question

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

Answer 1

Following, a way to remove the 3 exponential pulses.

1.- Rebuilding signal + interferer from graph provided

clear all;clc;close all

A=imread('001.jpg');
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

Comments:

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.

Answer 2

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.