Deleting spikes and reconstructing background noise

1 次查看（过去 30 天）
Bence Laczó
Bence Laczó 2023-6-8
评论： Star Strider 2023-6-10
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
  2 个评论
Star Strider
Star Strider 2023-6-8
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ó
Bence Laczó 2023-6-8
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?

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Star Strider
Star Strider 2023-6-8
Ran in:
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.
.
  2 个评论
Bence Laczó
Bence Laczó 2023-6-10
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
Star Strider 2023-6-10
Ran in:
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.
.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Descriptive Statistics 的更多信息

标签

产品


版本

R2019b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by