obspy.signal.interpolation.lanczos_interpolation
- lanczos_interpolation(data, old_start, old_dt, new_start, new_dt, new_npts, a, window='lanczos', *args, **kwargs)[source]
Function performing Lanczos resampling, see https://en.wikipedia.org/wiki/Lanczos_resampling for details. Essentially a finite support version of sinc resampling (the ideal reconstruction filter). For large values of
ait converges towards sinc resampling. If used for downsampling, make sure to apply an appropriate anti-aliasing lowpass filter first.Note
In most cases you do not want to call this method directly but invoke it via either the
obspy.core.stream.Stream.interpolate()orobspy.core.trace.Trace.interpolate()method. These offer a nicer API that naturally integrates with the rest of ObsPy. Usemethod="lanczos"to use this interpolation method. In that case the only additional parameters of interest areaandwindow.- Parameters:
data (array_like) – Array to interpolate.
old_start (float) – The start of the array as a number.
old_dt – The time delta of the current array.
new_start (float) – The start of the interpolated array. Must be greater or equal to the current start of the array.
new_dt (float) – The desired new time delta.
new_npts (int) – The new number of samples.
a (int) – The width of the window in samples on either side. Runtime scales linearly with the value of
abut the interpolation also gets better.window (str) – The window used to taper the sinc function. One of
"lanczos","hanning","blackman". The window determines the trade-off between “sharpness” and the amplitude of the wiggles in the pass and stop band. Please use theplot_lanczos_windows()function to judge these for any given application.
Values of
a>= 20 show good results even for data that has energy close to the Nyquist frequency. If your data is extremely oversampled you can get away with much smallera’s.To get an idea of the response of the filter and the effect of the different windows, please use the
plot_lanczos_windows()function.Also be aware of any boundary effects. All values outside the data range are assumed to be zero which matters when calculating interpolated values at the boundaries. At each side the area with potential boundary effects is
a*old_dt. If you want to avoid any boundary effects you will have to remove these values.Mathematical Details:
The \(\operatorname{sinc}\) function is defined as
\[\operatorname{sinc}(t) = \frac{\sin(\pi t)}{\pi t}.\]The Lanczos kernel is then given by a multiplication of the \(\operatorname{sinc}\) function with an additional window function resulting in a finite support kernel.
\[\begin{split}\begin{align} L(t) = \begin{cases} \operatorname{sinc}(t)\, \cdot \operatorname{sinc}(t/a) & \text{if } t \in [-a, a] \text{ and } \texttt{window} = \texttt{lanczos}\\ \operatorname{sinc}(t)\, \cdot \frac{1}{2} (1 + \cos(\pi\, t/a)) & \text{if } t \in [-a, a] \text{ and } \texttt{window} = \texttt{hanning}\\ \operatorname{sinc}(t)\, \cdot \left( \frac{21}{50} + \frac{1}{2} \cos(\pi\, t/a) + \frac{2}{25} \cos (2\pi\, t/a) \right) & \text{if } t \in [-a, a] \text{ and } \texttt{window} = \texttt{blackman}\\ 0 & \text{else} \end{cases} \end{align}\end{split}\]Finally interpolation is performed by convolving the discrete signal \(s_i\) with that kernel and evaluating it at the new time samples \(t_j\):
\[\begin{align} S(t_j) = \sum_{i = \left \lfloor{t_j / \Delta t}\right \rfloor -a + 1} ^{\left \lfloor{t_j / \Delta t}\right \rfloor + a} s_i L(t_j/\Delta t - i), \end{align}\]where \(\lfloor \cdot \rfloor\) denotes the floor function. For more details and justification please see [Burger2009] and [vanDriel2015].