numpy-discussion

[Numpy-discussion] interpolating arrays (?)

[Ray S. <ra...@bl...>]

Howdy,

I'd like to resize an array (1-dimensional) a of length n into an array of 
n+x, and interpolate the values as necessary.
x might be+/-200, len(a) is usually ~500

I tried it in pure Python but it's not quite right yet. (rubberBand, below)
Someone out there must have done this... and linear interpolation for 
values is probably sufficient (instead of bi-cubic or some such)

Cheers,
Ray




==================================================================
         startSlice = 10000
         endOfSlice = 
10000+int(round(self.samplesPerRotation))    #self.samplesPerRotation=661
         ## self.dataArray is about len=200,000 of Int

         for delta in range(-2,3):
             ## interpolate the different slices to the original length
	  ##  0 should return the same array, this is a test
             newArray = self.rubberBand(
                                 self.dataArray[startSlice:endOfSlice+delta],
                                 round(self.samplesPerRotation))
             print len(self.dataArray[startSlice:endOfSlice+delta]), 
len(newArray),
             ## show the result
             print delta, Numeric.sum(self.dataArray[startSlice:endOfSlice] 
- newArray)

     def rubberBand(self, data, desiredLength):
         """ alias/interpolate the data so that it is adjusted to the 
desired array length"""
         currentLength = float(len(data))

         ## positive if the new array is to be longer
         difference = desiredLength - currentLength
         #print  difference, desiredLength, currentLength

         ## set up the desired length array
         newData = Numeric.zeros(desiredLength, Numeric.Float)

         ## accounts for binary rounding errors
         smallNum = 2**-14

         for index in range(desiredLength):
             ## find the ratio of the current index to the desired length
             ratio = index / float(desiredLength)

             ## find the same (float) position in the old data
             currIndex = int(ratio * currentLength)
             ## find the decimal part
             decimalPart = (ratio * currentLength) - currIndex
             print index, currIndex, decimalPart,
             if(decimalPart>smallNum):
                 ## interpolate
                 newData[index] = ((1 - decimalPart) * data[currIndex] + 
decimalPart * data[currIndex+1])
                 print data[currIndex], data[currIndex+1], newData[index]
             else:
                 newData[index] = data[currIndex]
                 print 'else',data[currIndex], newData[index]

         return newData

Re: [Numpy-discussion] interpolating arrays (?)

[Warren F. <fo...@sl...>]

Would you consider Fourier interpolation?

interpolated = FFT.inverse_real_fft(FFT.real_fft(data), desiredLength) * \
	desiredLength / len(data)

This works whether desiredLength is less or greater than len(data).  The
trick here is that inverse_real_fft discards data or zero-pads at high
frequencies.  I keep toying with the idea of modifying the complex
version to do the same, I cannot imagine a use for the current behavior
(truncate or pad at small negative frequencies).

There can be some pretty wierd edge effects, they can often be defeated
with padding, or, better yet, overlapping segments.  Geting that right
might be harder than fixing your approach.

It has the nice feature that it is information-preserving:  if you
interpolate to a higher number of points, then back to the original, you
should get the same numbers back, give or take a little bit of rounding
error.

Warren Focke

On Mon, 22 Mar 2004, Ray Schumacher wrote:

> Howdy,
>
> I'd like to resize an array (1-dimensional) a of length n into an array of
> n+x, and interpolate the values as necessary.
> x might be+/-200, len(a) is usually ~500
>
> I tried it in pure Python but it's not quite right yet. (rubberBand, below)
> Someone out there must have done this... and linear interpolation for
> values is probably sufficient (instead of bi-cubic or some such)
>
> Cheers,
> Ray
>
>
>
>
> ==================================================================
>          startSlice = 10000
>          endOfSlice =
> 10000+int(round(self.samplesPerRotation))    #self.samplesPerRotation=661
>          ## self.dataArray is about len=200,000 of Int
>
>          for delta in range(-2,3):
>              ## interpolate the different slices to the original length
> 	  ##  0 should return the same array, this is a test
>              newArray = self.rubberBand(
>                                  self.dataArray[startSlice:endOfSlice+delta],
>                                  round(self.samplesPerRotation))
>              print len(self.dataArray[startSlice:endOfSlice+delta]),
> len(newArray),
>              ## show the result
>              print delta, Numeric.sum(self.dataArray[startSlice:endOfSlice]
> - newArray)
>
>      def rubberBand(self, data, desiredLength):
>          """ alias/interpolate the data so that it is adjusted to the
> desired array length"""
>          currentLength = float(len(data))
>
>          ## positive if the new array is to be longer
>          difference = desiredLength - currentLength
>          #print  difference, desiredLength, currentLength
>
>          ## set up the desired length array
>          newData = Numeric.zeros(desiredLength, Numeric.Float)
>
>          ## accounts for binary rounding errors
>          smallNum = 2**-14
>
>          for index in range(desiredLength):
>              ## find the ratio of the current index to the desired length
>              ratio = index / float(desiredLength)
>
>              ## find the same (float) position in the old data
>              currIndex = int(ratio * currentLength)
>              ## find the decimal part
>              decimalPart = (ratio * currentLength) - currIndex
>              print index, currIndex, decimalPart,
>              if(decimalPart>smallNum):
>                  ## interpolate
>                  newData[index] = ((1 - decimalPart) * data[currIndex] +
> decimalPart * data[currIndex+1])
>                  print data[currIndex], data[currIndex+1], newData[index]
>              else:
>                  newData[index] = data[currIndex]
>                  print 'else',data[currIndex], newData[index]
>
>          return newData
>
>
>
> -------------------------------------------------------
> This SF.Net email is sponsored by: IBM Linux Tutorials
> Free Linux tutorial presented by Daniel Robbins, President and CEO of
> GenToo technologies. Learn everything from fundamentals to system
> administration.http://ads.osdn.com/?ad_id=1470&alloc_id=3638&op=click
> _______________________________________________
> Numpy-discussion mailing list
> Num...@li...
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>

Re: [Numpy-discussion] interpolating arrays (?)

[RayS <ra...@bl...>]

At 08:19 PM 3/22/04 -0800, Warren Focke wrote:

>Would you consider Fourier interpolation?

I'll try it out tomorrow and check the results.
Thanks!
Ray

>interpolated = FFT.inverse_real_fft(FFT.real_fft(data), desiredLength) * \
>         desiredLength / len(data)
>
>This works whether desiredLength is less or greater than len(data).  The
>trick here is that inverse_real_fft discards data or zero-pads at high
>frequencies.  I keep toying with the idea of modifying the complex
>version to do the same, I cannot imagine a use for the current behavior
>(truncate or pad at small negative frequencies).
>
>There can be some pretty wierd edge effects, they can often be defeated
>with padding, or, better yet, overlapping segments.  Geting that right
>might be harder than fixing your approach.
>
>It has the nice feature that it is information-preserving:  if you
>interpolate to a higher number of points, then back to the original, you
>should get the same numbers back, give or take a little bit of rounding
>error.
>
>Warren Focke
>
>On Mon, 22 Mar 2004, Ray Schumacher wrote:
>
> > Howdy,
> >
> > I'd like to resize an array (1-dimensional) a of length n into an array of
> > n+x, and interpolate the values as necessary.
> > x might be+/-200, len(a) is usually ~500
> >
> > I tried it in pure Python but it's not quite right yet. (rubberBand, below)
> > Someone out there must have done this... and linear interpolation for
> > values is probably sufficient (instead of bi-cubic or some such)
> >
> > Cheers,
> > Ray
> >
> >
> >
> >
> > ==================================================================
> >          startSlice = 10000
> >          endOfSlice =
> > 10000+int(round(self.samplesPerRotation))    #self.samplesPerRotation=661
> >          ## self.dataArray is about len=200,000 of Int
> >
> >          for delta in range(-2,3):
> >              ## interpolate the different slices to the original length
> >         ##  0 should return the same array, this is a test
> >              newArray = self.rubberBand(
> > 
> self.dataArray[startSlice:endOfSlice+delta],
> >                                  round(self.samplesPerRotation))
> >              print len(self.dataArray[startSlice:endOfSlice+delta]),
> > len(newArray),
> >              ## show the result
> >              print delta, Numeric.sum(self.dataArray[startSlice:endOfSlice]
> > - newArray)
> >
> >      def rubberBand(self, data, desiredLength):
> >          """ alias/interpolate the data so that it is adjusted to the
> > desired array length"""
> >          currentLength = float(len(data))
> >
> >          ## positive if the new array is to be longer
> >          difference = desiredLength - currentLength
> >          #print  difference, desiredLength, currentLength
> >
> >          ## set up the desired length array
> >          newData = Numeric.zeros(desiredLength, Numeric.Float)
> >
> >          ## accounts for binary rounding errors
> >          smallNum = 2**-14
> >
> >          for index in range(desiredLength):
> >              ## find the ratio of the current index to the desired length
> >              ratio = index / float(desiredLength)
> >
> >              ## find the same (float) position in the old data
> >              currIndex = int(ratio * currentLength)
> >              ## find the decimal part
> >              decimalPart = (ratio * currentLength) - currIndex
> >              print index, currIndex, decimalPart,
> >              if(decimalPart>smallNum):
> >                  ## interpolate
> >                  newData[index] = ((1 - decimalPart) * data[currIndex] +
> > decimalPart * data[currIndex+1])
> >                  print data[currIndex], data[currIndex+1], newData[index]
> >              else:
> >                  newData[index] = data[currIndex]
> >                  print 'else',data[currIndex], newData[index]
> >
> >          return newData
> >
> >
> >
> > -------------------------------------------------------
> > This SF.Net email is sponsored by: IBM Linux Tutorials
> > Free Linux tutorial presented by Daniel Robbins, President and CEO of
> > GenToo technologies. Learn everything from fundamentals to system
> > administration.http://ads.osdn.com/?ad_id=1470&alloc_id=3638&op=click
> > _______________________________________________
> > Numpy-discussion mailing list
> > Num...@li...
> > https://lists.sourceforge.net/lists/listinfo/numpy-discussion
> >

Re: [Numpy-discussion] interpolating arrays (?)

[Chris B. <Chr...@no...>]

RayS wrote:
> At 08:19 PM 3/22/04 -0800, Warren Focke wrote:
>> Would you consider Fourier interpolation?
> I'll try it out tomorrow and check the results.
> Thanks!

I woulnd't use Fourier interpolation unles your data support that. i.e. 
if it is periodic.

SciPy has an interpolation module. I'd download it and check it out. 
Last I tried, you can use the interpolation module without installing 
the whole SciPy package.

-Chris


-- 
Christopher Barker, Ph.D.
Oceanographer
                                     		
NOAA/OR&R/HAZMAT         (206) 526-6959   voice
7600 Sand Point Way NE   (206) 526-6329   fax
Seattle, WA  98115       (206) 526-6317   main reception

Chr...@no...

Re: [Numpy-discussion] interpolating arrays (?)

[RayS <ra...@bl...>]

Hi Tim,

I implemented the below code with good results; and if I go to a C version 
in the future I'll let you know. As the raw data is integer and a bit 
coarse at times, it will also help to have a non-linear method. I use a 
parabolic interpolation for finding centriods of clipped stellar images, 
but that's still Python too.

>With that one comes up with (untested):
>
>def rubberBand(self, y, desiredLength):
>       # Define raw so that raw[0] == 0 and raw[-1] == len(y)-1 and 
> len(raw) == desiredLength
>        raw =arange(desiredLength) * (len(y) - 1) / (float(desiredLength) - 1)
>        jVals = na.clip(na.floor(raw), 0, len(y)-2).astype('i')
>        delta = raw - jVals
>        dy = y[1:] - y[:-1]
>        return na.take(y, jVals) + delta * na.take(dy, jVals)
>
>Hope that's helpful,

quite.

>PS, I just realized, these snippets assume a 'import numarray as na' 
>somewhere above. Numeric should also work if the names are adjusted 
>appropriately.

It does. I use Numeric for when arrays are length< 2000, as someone posted 
some results to that effect a while back. numarray is certainly faster for 
images.

In general, I'm using FFT for ballpark estimation of periodicity, then 
doing time domain data comparison for more precise alignment.

Thanks to Chris, Konrad and Warren as well,

Ray