It's impossible to say what the best method is without a model of spike formation and a function that quantifies the cost of residual errors in the smoothed data. Without these, you could try some techniques and try to assess the performance.
There are many different possibilities. Every method that reduces the noise will do some damage to the data, and the problem is to find which is the best trade-off.
A reasonable starting point might be median filtering, using medfilt2() and medfilt1(). You could also look at more general nonlinear filters, which you can implement using nlfilter().
Another possibility is linear filtering. For example, you can smooth with a Gaussian kernel using fspecial() and conv2(). There's an easy interface in Gradients with Gaussian Smoothing in the File Exchange.
There are many more advanced methods if these aren't adequate, but it will involve experiment and research to get it right.
