Deleting spikes and reconstructing background noise

Hi!

I would like to analyse local field potentials recorded by a microelectrode. To do this I need to delete the spikes and replace them with background data from a random location that does not contain any spikes. I have the exact locations of the spikes in a variable named spikeMark. I made the following code which seems working, but I would like to ask your opinion about it, and your suggestions how to make it more efficient. The raw data is in the data variable.

Thanks for your help in advance!

for i=1:length(data)           % deleting spikes
    if spikeMark(i) == 1
        data(i-12:i+60) = 0;
    end
end

y = buffer(data,73,72,'nodelay');  %replacing spikes with background noise
for i=1:length(y)
    B(i) = any(y(:,i) == 0);
end

y(:, B) = [];
availableToUse = y;

for k = 1:length(data)                                              
    if data(k) == 0
        randomIndex = randi(size(availableToUse,2), 1, 1);
        data(k:k+72) = availableToUse(:,randomIndex);
    end
end

[Star Strider]

Instead of setting the spikes (and apparently their surrounding data) to 0, it might be easier to set them to NaN and then use the fillmissing function (introduced in R2016b) to interpolate them.

I’m not posting this as an answer because I don’t have your data to experiment with.

[Bence Laczó]

I have added one data file (sampling rate is 24000Hz), and added the positions of the detected spikes (this data is in msec). Can you post me how the fillmissing function can be used?

Answer 1

[Star Strider]

While it is possible to replace the spikes with broadband noise, here I just replaced them with a short sine segment —

LD1 = load('central+4.mat');
data = LD1.data;
L = numel(data);
Fs = 2.4E+5;
t = linspace(0, L-1, L)/Fs;

LD2 = load('timestamp.mat');
spikeMark = LD2.timestamp;

figure
plot(t, data)
hold on
plot(t(spikeMark), data(spikeMark), 'r+')
hold off
grid

figure
plot(t, data)
hold on
plot(t(spikeMark), data(spikeMark), 'r+')
hold off
grid
xlim([0 0.001])

figure
plot(t, data)
hold on
plot(t(spikeMark), data(spikeMark), 'r+')
hold off
grid
xlim([0 0.001]+0.3)
ylim([-1 1]*80)

Fn = Fs/2;
NFFT = 2^nextpow2(L)

NFFT = 262144

FTdata = fft((data-mean(data)).*hann(L), NFFT)/L;
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);

figure
plot(Fv, abs(FTdata(Iv))*2)
grid
xlim([0 6]*1E+4)

for k = 1:numel(spikeMark)
    idxrng = max(1,spikeMark(k)-12) : min(L, spikeMark(k)+20);      % Changed Range, Included Limits Check
    % data(idxrng) = NaN;
    data(idxrng) = 10*sin(t(1:numel(idxrng))*2*pi*2E+4);
end

figure
plot(t, data)
hold on
plot(t(spikeMark), data(spikeMark), 'r+')
hold off
grid
xlim([0 0.001])
ylim([-1 1]*80)

figure
plot(t, data)
hold on
plot(t(spikeMark), data(spikeMark), 'r+')
hold off
grid
xlim([0 0.001]+0.3)
ylim([-1 1]*80)

The fillmissing function would work here (I tried it first), however it would do a simpler type of interpolation (linear or other methods), over the gap created by the NaN values, rather than replacing them with something more in keeping with what you want. Here, I did everything in one loop.

.

[Bence Laczó]

First I would like to thank you for your help! That code is really nice I really like how simply you solved the problem. But I really need to use the background noise to replace the spikes as I would like to reproduce the data from a former study.

[Star Strider]

It seems that you are using the randi function to create the background noise, however this is actually a version of ‘white’ noise. It will not have the same spectral characteristics as tthe rest of your signal.

One option would be to simply duplicate the segment of the signal immediately following the deleted sections, providing that does not index beyond the end of the vector, otherwise use the preceding section instead —

LD1 = load('central+4.mat');
data = LD1.data;
L = numel(data);
Fs = 2.4E+5;
t = linspace(0, L-1, L)/Fs;

LD2 = load('timestamp.mat');
spikeMark = LD2.timestamp;

figure
plot(t, data)
hold on
plot(t(spikeMark), data(spikeMark), 'r+')
hold off
grid

figure
plot(t, data)
hold on
plot(t(spikeMark), data(spikeMark), 'r+')
hold off
grid
xlim([0 0.001])

figure
plot(t, data)
hold on
plot(t(spikeMark), data(spikeMark), 'r+')
hold off
grid
xlim([0 0.001]+0.3)
ylim([-1 1]*80)

Fn = Fs/2;
NFFT = 2^nextpow2(L)

NFFT = 262144

FTdata = fft((data-mean(data)).*hann(L), NFFT)/L;
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);

figure
plot(Fv, abs(FTdata(Iv))*2)
grid
xlim([0 6]*1E+4)

for k = 1:numel(spikeMark)
    idxrng = max(1,spikeMark(k)-12) : min(L, spikeMark(k)+20);      % Changed Range, Included Limits Check
    % data(idxrng) = NaN;
    % data(idxrng) = 10*sin(t(1:numel(idxrng))*2*pi*2E+4 + 2*pi*rand(size(idxrng))-pi);
    if all((idxrng+numel(idxrng)) < L)
        data(idxrng) = data(idxrng+numel(idxrng));
    else
        data(idxrng) = data(idxrng-numel(idxrng));
    end
end

figure
plot(t, data)
hold on
plot(t(spikeMark), data(spikeMark), 'r+')
hold off
grid
xlim([0 0.001])
ylim([-1 1]*80)

figure
plot(t, data)
hold on
plot(t(spikeMark), data(spikeMark), 'r+')
hold off
grid
xlim([0 0.001]+0.3)
ylim([-1 1]*80)

The advantage of this (as I see it) is that it retains the spectral characteristics of the vector. From a signal processing perspective, this is important.

.