Clean noise in neural data - MATLAB
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&gt;t0-.01),find(t&lt;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.