If you have the Signal Processing Toolbox, the findpeaks function (with appropriate name-value pair argument choices and parameters) could help you identify the spikes and their locations.
Keep your original independent variable vector, since you’ll need it for the interpolation. Then use the location index vector provided by findpeaks to eliminate them. Then interpolate.
To illustrate:
x = 1:100; % Create Data
y = 1 + 0.1*x.*rand(1,100); % Create Data
spikes = randi(100, 1, 5); % Create Data
y(spikes) = randi(50, 1, 5); % Create Data
[spk,loc] = findpeaks(y,'MinPeakHeight',10); % Find ‘spikes’
xorg = x; % Save Original ‘x’
yorg = y; % Save Original ‘y’
x(loc) = []; % Eliminate ‘spikes’ From Data
y(loc) = []; % Eliminate ‘spikes’ From Data
yint = interp1(x, y, xorg); % Interpolate Missing Data
figure(1)
subplot(2,1,1)
plot(xorg, yorg)
grid
title('Original Data')
yl = get(gca, 'YLim')
subplot(2,1,2)
plot(xorg, yint)
grid
title('Interpolated Data Without Spikes')
set(gca, 'YLim',yl)