Hi Matthew,
On Tue, Feb 22, 2011 at 9:42 PM, Matthew Marcus
Should the line out[i] = (x[i+1] - x[i])/2.0 be out[i] = (x[i+1] - x[i-1])/2.0 ?
That would account for the 2.0. That is, in fact, exactly what I proposed was meant by 'numerical differential'.
Yes, sorry for the typo, it should b def deriv(x): npts = len(x) out = array(npts) out[0] = x[1] - x[0] out[npts-1] = x[npts-1] - x[npts-2] for i in range(1, npts-2): # (1, 2, 3, ... npts-2) out[i] = (x[i+1] - x[i-1])/2.0 return out
I don't actually resample. I do a cubic spline interpolation onto the given grid and pull the derivatives out that way. Of course, you have to test for grids that aren't strictly monotonic. The smooth() function as given would definitely cause problems with non-uniform tabulation unless resampled. My program (2-column Editor on the ALS Beamline 10.3.2 website) offers resampling and tensioned spline fitting but does not smooth the derivative before or after - the user has to take control of that using the tensioned spline fit. I thought about smoothing but never got around to figuring out what makes a decent smooth algorithm for general data. A possibility would be a slow version of Savitzky-Golay in which you do the polynomial fits explicitly instead of being able to use pre-computed coefficients as in the actual S-G method. Similarly, a VSFT (Very Slow Fourier Transform) in which you fit the data to a sum of sines and cosines (linear fit - the frequencies are fixed), filter using your favorite low-pass filter, then inverse-transform, might be good. Or, there's always reading the literature to see what actual numerical mathematicians have done :-)
Smooth() is definitely over-simplistic and will work poorly for non-uniform grids. Using a Savitzky-Golay smoothing function is probably a good idea. Having the abilty to use more sophisticated numerical algorithms on demand would be very nice.... --Matt