I'd first try to display the running average over a number of points, like 5 or 10. This way, a single discrepancy in the values only have a little impact on the graph. Of course, it depends on how accurate you need the graph to be.
A simple (ad hoc) way is to just take a weighted average (tunable by alpha) at each point with its neighbors:
data(2:n-1) = alpha*data(2:n-1) + (1-alpha)*0.5*(data(1:n-2)+data(3:n))
or some variation thereof. Yes, to be more sophisticated you can Fourier transform your data first, then cut off the high frequencies. Something like:
f = fft(data)
f(n/2+1-20:n/2+20) = zeros(40,1)
smoothed = real(ifft(f))
This cuts out the highest 20 frequencies. Be careful to cut them out symmetrically otherwise the inverse transform is no longer real. You need to carefully choose the cutoff frequency for the right level of smoothing. This is a very simple kind of filtering (box filtering in frequency domain), so you can try gently attenuating high order frequencies if the distortion is unacceptable.
FFT isn't a bad idea, but it's probably overkill here. Running or moving averages give generally poor results and should be avoided for anything besides late homework (and white noise).
I'd use Savitzky-Golay filtering (in Matlab sgolayfilt(...)). This will give you the best results for what you are looking for - some local smoothing while maintaining the shape of the curve.
-Paul