MATLAB: How can ı do spike sorting to filtered EEG signal.

eegspikethreshold

I have this sample and ı need to do spike sorting. I start with finding thresholding value and want to see changes with this code. I have trouble with plotting in for loop and also further steps. I am very new in Matlab. If you can help it would be great thank you.

for k = 1:0.1:5
    
omg = median(abs(samp)/0.6745);
Thr =k * omg;
[pks, locs] = findpeaks(samp);
findpeaks(samp, 'Threshold', Thr);
[r,NumPeaks] = size(pks);
NumPeaks
end
hold on

Best Answer

Try this:

D = load('samp.mat');
samp = D.samp;
data = 1E-3;                                                    % What Are ‘data’?
omg = median(abs(data)/0.6745);
% th = 1:0.1:5;
th = 5 : -0.1 : 1;                                              % Easier To Count Down Here
for k = 1:numel(th)
    Thr = th(k) * omg;
    [pks{k}, locs{k}] = findpeaks(samp, 'Threshold', Thr);      % Peaks Above This ‘Thr’ Value
    kprv = max(k-1,1);                                          % Previous Index Value
    lv = ~ismember(locs{k},locs{kprv});                         % Spikes Above Current ‘Thr’ & Below Previous ‘Thr’
    nspikes(k) = nnz(lv);                                       % Number Of Spikes Above Current ‘Thr’ & Below Previous ‘Thr’
    ThrPks{k} = pks{k}(lv);                                     % Peak Amplitudes Above Current ‘Thr’ & Below Previous ‘Thr’
    ThrLcs{k} = locs{k}(lv);                                    % Locations Of Peak Amplitudes Above Current ‘Thr’ & Below Previous ‘Thr’
end
figure
bar(th, nspikes)
xlabel('Threshold')
ylabel('# Spikes')

I am not certain what you want, so I saved everything about your data in my code.

Related Solutions

MATLAB: How to select a specific peak in a repeating pattern

Try this:

clc;  % Clear command window.
clear;  % Delete all variables.
close all;  % Close all figure windows except those created by imtool.
workspace;  % Make sure the workspace panel is showing.
fontSize = 24;
s = load('dataa.mat')
z = s.z;
plot(z, 'b-', 'LineWidth', 2);
grid on;
hold on;
title('z vs. Index', 'FontSize', fontSize);
xlabel('Index', 'FontSize', fontSize);
ylabel('z', 'FontSize', fontSize);
% Label different squarish regions.
[labeledSignal, numRegions] = bwlabel(z > 200);
% Find indexes of each group
props = regionprops(labeledSignal, 'PixelIdxList');
% Examine each group to find out where the signal decreases
for k = 1 : numRegions
  % Get indexes of only this group.
  theseIndexes = props(k).PixelIdxList;
  % Plot a red line where the group starts.
  line([theseIndexes(1), theseIndexes(1)], ylim, 'Color', 'r', 'LineWidth', 2);
  % Get differences.  If it's going up, differences will be positive. 
  % When it turns around, the difference values will be zero or negative.
  diffs = diff(z(theseIndexes));
  % Find out where diffs <= 0
  peakIndex(k) = find(diffs <= 0, 1, 'first') + theseIndexes(1) - 1;
  % Tell user what the index is.
  fprintf('The peak for group #%d is at index %d.\n', k, peakIndex(k));
  % Plot a magenta line where the first peak in the group is.
  line([peakIndex(k), peakIndex(k)], ylim, 'Color', 'm', 'LineWidth', 2);
end
  
%------------------------------------------------------------------------------
% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0, 0.04, 1, 0.96]);
% Get rid of tool bar and pulldown menus that are along top of figure.
% set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off')

Each group starts at the red line, and the very first peak has its highest point at the magenta line.

MATLAB: Finding points of signal

Well you could start at a peak and go forward and backward until the signal begins to turn upward.

Or you could just use the widths returned by findpeaks(), which is about half way up. I have a demo for this:

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 22;
fprintf('Beginning to run %s.m ...\n', mfilename);
EOG_Data = load('EOG_Data.mat');
EOG_Data_v = double(EOG_Data.EOG_Data_4ch(1,:)) - double(EOG_Data.EOG_Data_4ch(2,:));
fs = 256;
N = length(EOG_Data_v);
T = 1/fs;
EOG_Time = [0:(N - 1)] * T;
subplot(2, 1, 1);
plot(EOG_Time, EOG_Data_v);
title('x_v[n]', 'FontSize', fontSize);
grid on; 
Blink = EOG_Data_v(500:720); %template
[R, lags] = xcorr(EOG_Data_v, Blink); %cross-correlation
R = R(lags>=0);
lags = lags(lags>=0);
[peakValues, indexesOfPeaks, widths, proms] = findpeaks(R, 'MinPeakHeight',10e6); %finding peaks
subplot(2, 1, 2);
plot(EOG_Time, R);
title('Occurance of Blinks from Cross-Correlation plot', 'FontSize', fontSize);
xlabel('Time in seconds', 'FontSize', fontSize); 
ylabel('Voltage in mV', 'FontSize', fontSize); 
grid on; 
axis tight;
hold on;
plot(EOG_Time(indexesOfPeaks), peakValues, '*r')
% Plot widths
valleyIndexes = zeros(length(widths), 2);
for k = 1 : length(widths)
	% Get the indexes.
	valleyIndexes(k, 1) = floor(indexesOfPeaks(k) - widths(k)/2);
	valleyIndexes(k, 2) = ceil(indexesOfPeaks(k) + widths(k)/2);
	% Get the times.
	t1 = EOG_Time(valleyIndexes(k, 1));
	t2 = EOG_Time(valleyIndexes(k, 2));
	% Draw vertical lines on the plot.
	xline(t1, 'Color', 'm', 'LineWidth', 2);
	xline(t2, 'Color', 'm', 'LineWidth', 2);
end
hold off;
% Maximize figure window.
g = gcf;
g.WindowState = 'maximized'
fprintf('Done running %s.m.\n', mfilename);

Best Answer

Related Solutions

MATLAB: How to select a specific peak in a repeating pattern

MATLAB: Finding points of signal

Related Question