[SciPy-user] Sample and hold

[Robert Clewley]

Hi Martin,

On Sat, 3 Jun 2006, Martin Spacek wrote:
>
> I'd like the output to have regularly spaced data points. This means
> that the exact times of the original data points will be lost, but
> that's OK.

FWIW, you might be interested in a slightly more sophisticated solution: 
abstracting the data and treating it like a true curve. You can have 
arbitrary access to the curve for convenient resampling at your desired 
interval length.

For this you could use the interp0d class that I adapted from scipy's 
interp1d. You could also easily modify it to use the input points as the 
left or right-hand side of the constant intervals (mine uses them as right 
points, and averages the y values from the interval endpoints). Naturally, 
this solution involves more overhead than the methods previously posted.

Below is the modified __call__ method of interp1d (from scipy 0.3.2) that 
you need to create interp0d. Alternatively, the whole thing (including a 
curve class over these interpolators) can be found in the PyDSTool 
package. Note that the original data points are retained inside the class. 
A demo of it working with independent time and data arrays in this package 
also appears below.

-Rob


     def __call__(self,x_new):
         """Find piecewise-constant interpolated y_new = <name>(x_new).

         Inputs:
           x_new -- New independent variables.

         Outputs:
           y_new -- Piecewise-constant interpolated values corresponding to x_new.
         """
         # 1. Handle values in x_new that are outside of x.  Throw error,
         #    or return a list of mask array indicating the outofbounds values.
         #    The behavior is set by the bounds_error variable.
         x_new_1d = atleast_1d(x_new)
         out_of_bounds = self._check_bounds(x_new_1d)
         # 2. Find where in the orignal data, the values to interpolate
         #    would be inserted.
         #    Note: If x_new[n] = x[m], then m is returned by searchsorted.
         x_new_indices = searchsorted(self.x,x_new_1d)
         # 3. Clip x_new_indices so that they are within the range of
         #    self.x indices and at least 1.  Removes mis-interpolation
         #    of x_new[n] = x[0]
         x_new_indices = clip(x_new_indices,1,len(self.x)-1).astype(Int)
         # 4. Calculate the region that each x_new value falls in.
         lo = x_new_indices - 1; hi = x_new_indices

         # !! take() should default to the last axis (IMHO) and remove
         # !! the extra argument.
         # 5. Calculate the actual value for each entry in x_new.
         y_lo = take(self.y,lo,axis=self.interp_axis)
         y_hi = take(self.y,hi,axis=self.interp_axis)
         y_new = (y_lo+y_hi)/2.
         # 6. Fill any values that were out of bounds with NaN
         # !! Need to think about how to do this efficiently for
         # !! mutli-dimensional Cases.
         yshape = y_new.shape
         y_new = y_new.flat
         new_shape = list(yshape)
         new_shape[self.interp_axis] = 1
         sec_shape = [1]*len(new_shape)
         sec_shape[self.interp_axis] = len(out_of_bounds)
         out_of_bounds.shape = sec_shape
         new_out = ones(new_shape)*out_of_bounds
         putmask(y_new, new_out.flat, self.fill_value)
         y_new.shape = yshape
         # Rotate the values of y_new back so that they correspond to the
         # correct x_new values.
         result = swapaxes(y_new,self.interp_axis,self.axis)
         try:
             len(x_new)
             return result
         except TypeError:
             return result[0]
         return result


Example of use in PyDSTool, using the InterpolateTable class:

>> timeData = linspace(0, 10, 30)
>> sindata = sin(timeData)
>> xData = {'sinx': sindata}
>> sin_pcwc_interp = InterpolateTable({'tdata': timeData,
                               'ics': xData,
                               'name': 'interp0d',
                               'method': 'constant',
                               'checklevel': 1,
                               'abseps': 1e-5
                               }).compute('sin_interp')
>> sin_pcwc_interp(2*pi)  # should be close to zero!
0.094554071417120258
>> t = linspace(0,10,300) # new time array, much finer resolution
>> y = sin_pcwc_interp(t)  # new data (as pointset, so use array(y) to get array)