Downsampling / Resampling using scipy.signal.resample

[driel]

Downsampling by now is only done by lowpass filtering and integer decimation. Might be convinient to wrap scipy.signal.resample in Trace and Stream. As far as I see, scipy.signal.resample is equivalent to the exact reconstruction but done in frequency domain:

 http://en.wikipedia.org/wiki/Sampling_theorem#Reconstruction

 http://en.wikipedia.org/wiki/Whittaker%E2%80%93Shannon_interpolation_formula

Opinions, comments, experience with this?

Martin

[barsch]

+1 - but should be wrapped in obspy.signal to make sure that scipy as dependency is installed

[beyreuth]

yes, resample is a correct functions for this. We discussed this before but did not find time to implement it..., so +1 from too

[driel]

I read a bit more into the theory, and there is one difference to resampling using the sampling theoreme: in fourier domain, it is assumed, that the function is periodic. This leads to 'ringing', if the function does not go to zero at the boundaries. Zero padding might help, but I guess FFT makes detrending/tapering necessary to be stable in that case.

Also see  http://stackoverflow.com/questions/1851384/resampling-interpolating-matrix

On the other hand, in timedomain it is quite expensive (N² compared to N log(N)).

So, I suggest to wrap the scipy function and add some comment in the docstring...

[megies]

would be very handy to have for a lot of people i imagine. i would also go for the scipy frequency domain routine. if tapering/detrending is necessary i would include it controlled via kwargs with the defaults being True.

[barsch]

In [3442/obspy]:

• adding resample to Trace and Stream - fixes #314

[beyreuth]

The new trace method now holds:

1264            from scipy.signal import resample 
1265            factor = float(self.stats.npts) / num 
1266            self.data = resample(self.data, num, window=window) 
1267            self.stats.delta *= factor 

• I think we should use a hanning window here as default, using the default window None will show severe 'ringing' for non periodic signals as seismological data, as mentioned by driel above.
• There is no lowpassing applied here, this will result in severe aliasing. We must use a similar approach like ../../../obspy.core/obspy/core/trace.py line 1345 here.
• I'll try to look into at the end of the week.

[anonymous]

I would suggest to make the lowpass filter optional (maybe default), because in some applications it is necessary to use some specific filter and then it is better to do that manually before the resampling. I had some issues with the old decimate routine because of that.