From 961bc1affea79a70994bcad7e2fef8d63c237bf4 Mon Sep 17 00:00:00 2001 From: Scott Ransom Date: Sun, 17 Nov 2019 13:27:41 -0500 Subject: [PATCH] Fixed integer division and bad sinc implementation issues --- python/presto/sinc_interp.py | 35 +++++++++++++---------------------- 1 file changed, 13 insertions(+), 22 deletions(-) diff --git a/python/presto/sinc_interp.py b/python/presto/sinc_interp.py index 58c97810e..9eb22ee39 100644 --- a/python/presto/sinc_interp.py +++ b/python/presto/sinc_interp.py @@ -3,15 +3,6 @@ import numpy as Num import numpy.fft as FFT -def sinc(xs): - """ - sinc(xs): - Return the sinc function [i.e. sin(pi * xs)/(pi * xs)] - for the values xs. - """ - pxs = Num.pi*xs - return Num.where(Num.fabs(pxs)<1e-3, 1.0-pxs*pxs/6.0, Num.sin(pxs)/pxs) - def kaiser_window(xs, halfwidth, alpha): """ kaiser_window(xs, halfwidth, alpha): @@ -94,14 +85,14 @@ def windowed_sinc_interp(data, newx, halfwidth=None, if hi_pt >= len(data): hi_pt = len(data)-1 print("Warning: trying to access above the highest index!") - halfwidth = (hi_pt-lo_pt)/2 + halfwidth = (hi_pt-lo_pt)//2 pts = Num.arange(2*halfwidth)+lo_pt xs = newx - pts if window.lower() is "kaiser": - win = _window_function[window](xs, len(data)/2, alpha) + win = _window_function[window](xs, len(data)//2, alpha) else: - win = _window_function[window](xs, len(data)/2) - return Num.add.reduce(Num.take(data, pts) * win * sinc(xs)) + win = _window_function[window](xs, len(data)//2) + return Num.add.reduce(Num.take(data, pts) * win * Num.sinc(xs)) def periodic_interp(data, zoomfact, window='hanning', alpha=6.0): """ @@ -122,16 +113,16 @@ def periodic_interp(data, zoomfact, window='hanning', alpha=6.0): comb = Num.reshape(Num.transpose(comb), (newN,)) # Compute the offsets xs = Num.zeros(newN, dtype='d') - xs[:newN/2+1] = Num.arange(newN/2+1, dtype='d')/zoomfact - xs[-newN/2:] = xs[::-1][newN/2-1:-1] + xs[:newN//2+1] = Num.arange(newN//2+1, dtype='d')/zoomfact + xs[-newN//2:] = xs[::-1][newN//2-1:-1] # Calculate the sinc times window for the kernel if window.lower()=="kaiser": - win = _window_function[window](xs, len(data)/2, alpha) + win = _window_function[window](xs, len(data)//2, alpha) else: - win = _window_function[window](xs, len(data)/2) - kernel = win * sinc(xs) + win = _window_function[window](xs, len(data)//2) + kernel = win * Num.sinc(xs) if (0): - plotxy(sinc(xs), color='yellow') + plotxy(Num.sinc(xs), color='yellow') plotxy(win) plotxy(kernel, color='red') closeplot() @@ -155,8 +146,8 @@ def periodic_interp(data, zoomfact, window='hanning', alpha=6.0): # The "sampled" data Ndata = 100 - data = theo[::Ntheo/Ndata] - data_phases = theo_phases[::Ntheo/Ndata] + data = theo[::Ntheo//Ndata] + data_phases = theo_phases[::Ntheo//Ndata] # The values to interpolate Ncalc = 30 @@ -167,7 +158,7 @@ def periodic_interp(data, zoomfact, window='hanning', alpha=6.0): plotxy(data, data_phases, line=None, symbol=3, color='green') # Do the interpolation one point at a time - halfwidth = Ndata/2-5 + halfwidth = Ndata//2-5 calc_vals = [] for phs in calc_phases: calc_vals.append(windowed_sinc_interp(data, phs*len(data), halfwidth))