From eeb0ae7a5805c5ad403e2fdd1143e376dd3b592d Mon Sep 17 00:00:00 2001 From: t2 Date: Wed, 4 Jul 2018 16:37:16 +0200 Subject: [PATCH] Better interpolation Change-Id: I258a927fdc8912c0533df59ad832ddc739ddf146 --- brukeropusreader/opus_data.py | 51 +++++++---------------------------- 1 file changed, 9 insertions(+), 42 deletions(-) diff --git a/brukeropusreader/opus_data.py b/brukeropusreader/opus_data.py index fd3e648..6e2e13f 100644 --- a/brukeropusreader/opus_data.py +++ b/brukeropusreader/opus_data.py @@ -4,14 +4,6 @@ from numpy import linalg from scipy.interpolate import interp1d -class NoInterpolatedDataException(Exception): - pass - - -class WrongWavelengthsInConf(Exception): - pass - - class OpusData(object): def __init__(self, raw_data, interp=None, meta=None): self.raw_data = raw_data @@ -25,38 +17,13 @@ class OpusData(object): dist = linalg.norm(y_r - self.raw_data[1], ord=1) print("L1 distance between series:", dist) - def gen_interp_waves(self, conf): - wunit = (conf.wave_start - conf.wave_end) / (float(conf.iwsize - 1)) - a = conf.wave_start + wunit - iwavenumber = [0] * conf.iwsize - for i in range(len(iwavenumber)): - a -= wunit - iwavenumber[i] = a - return iwavenumber - - def interpolate_spectra(self, conf): + def interpolate(self, start, stop, num): xav, yav = self.raw_data[0], self.raw_data[1] - ab_wavenumbers = self.raw_data[0] - iwavenumber = self.gen_interp_waves(conf) - - xa_min = min(ab_wavenumbers) - xa_max = max(ab_wavenumbers) - n_interp = conf.iwsize - yi = [] - - f2 = interp1d(xav, yav, kind='cubic', assume_sorted=True) - - for k in range(n_interp): - if iwavenumber[k] < xa_min: - if k == n_interp - 1: - yi.append(yav[0]) - else: - print("Wrong wavelengths for interpolating data") - elif iwavenumber[k] > xa_max: - if k == 0: - yi.append(yav[-1]) - else: - raise WrongWavelengthsInConf() - else: - yi.append(f2(iwavenumber[k])) - return iwavenumber, yi + iwavenumber = np.linspace(start, + stop, + num) + f2 = interp1d(xav, yav, + kind='cubic', + assume_sorted=True, + fill_value='extrapolate') + return iwavenumber, f2(iwavenumber)