New basis interpolation function. Needs tests.

fsl_mrs/core/basis.py

@@ -80,6 +80,9 @@ class Basis:
         # This only has bearing on the plotting currently
         self._nucleus = '1H'
 
+        # Default interpolation is Fourier Transform based.
+        self._use_fourier_interp = True
+
     @classmethod
     def from_file(cls, filepath):
         """Create a Basis object from a path
@@ -174,6 +177,17 @@ class Basis:
         """Set the nucleus string - only affects plotting"""
         self._nucleus = nucleus
 
+    @property
+    def use_fourier_interp(self):
+        """Return interpolation state"""
+        return self._use_fourier_interp
+
+    @use_fourier_interp.setter
+    def use_fourier_interp(self, true_false):
+        """Set to true to use FFT based interpolation (default)
+        Or set to False to use time domain linear interpolation."""
+        self._use_fourier_interp = true_false
+
     def save(self, out_path, overwrite=False, info_str=''):
         """Saves basis held in memory to a directory in FSL-MRS format.
 
@@ -299,10 +313,16 @@ class Basis:
            coverage than the FID.
         """
         try:
-            basis = misc.ts_to_ts(self._raw_fids,
-                                  self.original_dwell,
-                                  target_dwell,
-                                  target_points)
+            if self.use_fourier_interp:
+                basis = misc.ts_to_ts_ft(self._raw_fids,
+                                         self.original_dwell,
+                                         target_dwell,
+                                         target_points)
+            else:
+                basis = misc.ts_to_ts(self._raw_fids,
+                                      self.original_dwell,
+                                      target_dwell,
+                                      target_points)
         except misc.InsufficentTimeCoverageError:
             raise BasisHasInsufficentCoverage('The basis spectra covers too little time. '
                                               'Please reduce the dwelltime, number of points or pad this basis.')

fsl_mrs/utils/misc.py

@@ -216,6 +216,72 @@ def ts_to_ts(old_ts, old_dt, new_dt, new_n):
     return new_ts
 
 
+def ts_to_ts_ft(old_ts, old_dt, new_dt, new_n):
+    """Temporal resampling using Fourier transform based resampling
+
+    Using the method implemented in LCModel:
+    1. Data is padded or truncated in the spectral domain to match the bandwidth of the target.
+    The ifft then returns the time domain data with the right overall length.
+    2. The data is then padded or truncated in the time domain to the length of the target.
+    If the data is then FFT it return the interpolated data.
+
+    :param old_ts: Input time-domain data
+    :type old_ts: numpy.ndarray
+    :param old_dt: Input dwelltime
+    :type old_dt: float
+    :param new_dt: Target dwell time
+    :type new_dt: float
+    :param new_n: Target number of points
+    :type new_n: int
+    :rtype: numpy.ndarray
+    """
+
+    old_n = old_ts.shape[0]
+    old_t = np.linspace(old_dt, old_dt * old_n, old_n) - old_dt
+    new_t = np.linspace(new_dt, new_dt * new_n, new_n) - new_dt
+    # Round to nanoseconds
+    old_t = np.round(old_t, 9)
+    new_t = np.round(new_t, 9)
+
+    if new_t[-1] > old_t[-1]:
+        raise InsufficentTimeCoverageError('Input data covers less time than is requested by interpolation.'
+                                           ' Change interpolation points or dwell time.')
+
+    def f2s(x):
+        return np.fft.fftshift(np.fft.fft(x, axis=0), axes=0)
+
+    def s2f(x):
+        return np.fft.ifft(np.fft.ifftshift(x, axes=0), axis=0)
+
+    # Input data to frequency domain
+    old_fs = f2s(old_ts)
+
+    # Step 1: pad or truncate in the frequency domain
+    new_bw = 1 / new_dt
+    old_bw = 1 / old_dt
+    npoints_f = (new_bw - old_bw) / (old_bw / old_ts.shape[0])
+    npoints_f_half = int(np.round(npoints_f / 2))
+
+    if npoints_f_half < 0:
+        # New bandwidth is smaller than old. Truncate
+        npoints_f_half *= -1
+        step1 = s2f(old_fs[npoints_f_half:-npoints_f_half])
+    elif npoints_f_half > 0:
+        # New bandwidth is larger than old. Pad
+        npoints_f_half
+        step1 = s2f(np.pad(old_fs, (npoints_f_half, npoints_f_half), 'constant', constant_values=(0j, 0j)))
+    else:
+        step1 = s2f(old_fs)
+
+    # Step 2: pad or truncate in the temporal domain
+    if step1.shape[0] < new_n:
+        step2 = np.pad(step1, (0, new_n - step1.shape[0]), 'constant', constant_values=(0j, 0j))
+    else:
+        step2 = step1[:new_n]
+
+    return step2
+
+
 # Numerical differentiation (light)
 # Gradient Function
 def gradient(x, f):