[Numpy-discussion] interpolating arrays (?)
RayS
rays at blue-cove.com
Mon Mar 22 21:39:00 EST 2004
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
> > Numpy-discussion at lists.sourceforge.net
> > https://lists.sourceforge.net/lists/listinfo/numpy-discussion
> >
More information about the NumPy-Discussion
mailing list