Repository: openmrslab/suspect

Issue:
Spectra should be first class data objects
Currently MRSData objects contain FID data and provide various convenience methods and properties, while spectra are bare ndarrays without any additional information. It would be better if the two classes inherited from (let's say) MRSBase giving some core methods, and then provided their own methods on top of those. Calling .spectrum() on an MRSData would return an MRSSpectrum object and calling .fid() on an MRSSpectrum object would return an MRSData object. Then users can work with whichever domain they prefer, without having to go back and forth all the time.

Discussion:
coveralls: 
[![Coverage Status](https://coveralls.io/builds/10772873/badge)](https://coveralls.io/builds/10772873)

Coverage increased (+0.2%) to 75.135% when pulling **a63bfd373567d86cbfe31d411b34652e2f3e284c on mrsobjects_inheritance** into **522105c8852802054a8849405e74c8cb8fa2bc60 on master**.


Patch:
diff --git a/requirements.txt b/requirements.txt
index 5af82e3..76bf75c 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -4,3 +4,4 @@ nbsphinx
 ipykernel
 numpy
 lmfit
+pydicom
diff --git a/suspect/__init__.py b/suspect/__init__.py
index 53dcc1a..5a58e22 100644
--- a/suspect/__init__.py
+++ b/suspect/__init__.py
@@ -1,4 +1,4 @@
-from suspect.mrsdata import MRSData
+from suspect.mrsobjects import MRSBase, MRSData, MRSSpectrum
 from suspect._transforms import *
 from . import processing
 from . import io
diff --git a/suspect/core.py b/suspect/core.py
index e667fe1..d34ae77 100644
--- a/suspect/core.py
+++ b/suspect/core.py
@@ -1,14 +1,11 @@
-from .mrsdata import MRSData
-
-
 def adjust_phase(data, zero_phase, first_phase=0, fixed_frequency=0):
     """
-    Adjust the phase of an MRSData object
+    Adjust the phase of an MRSBase object
 
     Parameters
     ----------
-    data : MRSData
-        The MRSData object to be phased
+    data : MRSSpectrum
+        The MRSSpectrum object to be phased
     zero_phase : scalar
         The change to the zero order phase, in radians
     first_phase : scalar, optional
@@ -19,15 +16,15 @@ def adjust_phase(data, zero_phase, first_phase=0, fixed_frequency=0):
 
     Returns
     -------
-    out : MRSData
-        A new MRSData object with adjusted phase.
+    out : MRSSpectrum
+        A new MRSSpectrum object with adjusted phase.
     """
     return data.adjust_phase(zero_phase, first_phase, fixed_frequency)
 
 
 def adjust_frequency(data, frequency_shift):
     """
-    Adjust the centre frequency of an MRSData object.
+    Adjust the centre frequency of an MRSBase object.
 
     Parameters
     ----------
diff --git a/suspect/mrsdata.py b/suspect/mrsobjects.py
similarity index 67%
rename from suspect/mrsdata.py
rename to suspect/mrsobjects.py
index ebe8a9f..a1fa3a9 100644
--- a/suspect/mrsdata.py
+++ b/suspect/mrsobjects.py
@@ -1,7 +1,7 @@
 import numpy
 
 
-class MRSData(numpy.ndarray):
+class MRSBase(numpy.ndarray):
     """
     numpy.ndarray subclass with additional metadata like sampling rate and echo
     time.
@@ -22,7 +22,7 @@ class MRSData(numpy.ndarray):
         return obj
 
     def __array_finalize__(self, obj):
-        # if this instance is being created by slicing from another MRSData, copy the parameters across
+        # if this instance is being created by slicing from another MRSBase, copy the parameters across
         self._dt = getattr(obj, 'dt', None)
         self._f0 = getattr(obj, 'f0', None)
         self._te = getattr(obj, 'te', 30)
@@ -38,25 +38,25 @@ class MRSData(numpy.ndarray):
             return numpy.ndarray.__array_wrap__(self, obj)
 
     def __str__(self):
-        return "<MRSData instance f0={0}MHz TE={1}ms dt={2}ms>".format(self.f0, self.te, self.dt * 1e3)
+        return "<MRSBase instance f0={0}MHz TE={1}ms dt={2}ms>".format(self.f0, self.te, self.dt * 1e3)
 
     def inherit(self, new_array):
-        """Converts a generic numpy ndarray into an MRSData instance by copying its own MRS specific parameters.
+        """Converts a generic numpy ndarray into an MRSBase instance by copying its own MRS specific parameters.
 
-        This is useful when performing some processing on the MRSData object gives a bare ndarray result.
+        This is useful when performing some processing on the MRSBase object gives a bare ndarray result.
 
         Parameters
         ----------
         new_array : numpy ndarray
-            Generic ndarray to be converted to MRSData.
+            Generic ndarray to be converted to MRSBase.
 
         Returns
         -------
-        cast_array : MRSData
-            New MRSData instance with data from new_array and parameters from self.
+        cast_array : MRSBase
+            New MRSBase instance with data from new_array and parameters from self.
 
         """
-        cast_array = new_array.view(MRSData)
+        cast_array = new_array.view(type(self))
         cast_array._dt = self.dt
         cast_array._f0 = self.f0
         cast_array._te = self.te
@@ -110,16 +110,6 @@ class MRSData(numpy.ndarray):
         """
         return self._f0
 
-    def spectrum(self):
-        """
-        Returns
-        -------
-        ndarray
-            The Fourier-transformed and shifted data
-
-        """
-        return numpy.fft.fftshift(numpy.fft.fft(self, axis=-1), axes=-1)
-
     def hertz_to_ppm(self, frequency):
         """Converts a frequency in Hertz to the corresponding PPM for this dataset.
 
@@ -199,7 +189,7 @@ class MRSData(numpy.ndarray):
         return numpy.prod(self.voxel_dimensions)
 
     def to_scanner(self, x, y, z):
-        """Converts a 3d position in MRSData space to the scanner reference frame
+        """Converts a 3d position in MRSBase space to the scanner reference frame
 
         Parameters
         ----------
@@ -217,14 +207,14 @@ class MRSData(numpy.ndarray):
 
         """
         if self.transform is None:
-            raise ValueError("No transform set for MRSData object {}".format(self))
+            raise ValueError("No transform set for MRSBase object {}".format(self))
 
         transformed_point = self.transform * numpy.matrix([x, y, z, 1]).T
 
         return numpy.squeeze(numpy.asarray(transformed_point))[0:3]
 
     def from_scanner(self, x, y, z):
-        """Converts a 3d position in the scanner reference frame to the MRSData space
+        """Converts a 3d position in the scanner reference frame to the MRSBase space
 
         Parameters
         ----------
@@ -238,37 +228,131 @@ class MRSData(numpy.ndarray):
         Returns
         -------
         ndarray
-            Squeezed ndarray representing a point in 3d MRSData space
+            Squeezed ndarray representing a point in 3d MRSBase space
 
         """
         if self.transform is None:
-            raise ValueError("No transform set for MRSData object {}".format(self))
+            raise ValueError("No transform set for MRSBase object {}".format(self))
 
         transformed_point = numpy.linalg.inv(self.transform) * numpy.matrix([x, y, z, 1]).T
 
         return numpy.squeeze(numpy.asarray(transformed_point))[0:3]
 
-    def adjust_phase(self, zero_phase, first_phase=0, fixed_frequency=0):
+
+class MRSData(MRSBase):
+    """
+    MRS data in the time domain.
+    """
+
+    def spectrum(self):
+        """
+        Returns
+        -------
+        MRSSpectrum
+            The Fourier-transformed and shifted data, represented as a spectrum
+
+        """
+        spectrum = self.inherit(numpy.fft.fftshift(numpy.fft.fft(self, axis=-1), axes=-1)).view(MRSSpectrum)
+        return spectrum
+
+    def adjust_phase(self, zero_phase, first_phase=0., fixed_frequency=0.):
         """
         Adjust the phases of the signal.
 
         Refer to suspect.adjust_phase for full documentation.
 
+        Parameters
+        ----------
+        zero_phase: float
+            The zero order phase shift in radians
+        first_phase: float
+            The first order phase shift in radians per Hertz
+        fixed_frequency: float
+            The frequency at which the first order phase shift is zero
+
+        Returns
+        -------
+        out : MRSData
+            Phase adjusted FID
+
         See Also
         --------
         suspect.adjust_phase : equivalent function
         """
-        # easiest to apply the phase shift in the frequency domain
-        # TODO when MRSSpectrum is a real class, this function can delegate
-        # to that one.
+        # easiest to do this in the spectral domain
         spectrum = self.spectrum()
+        return spectrum.adjust_phase(zero_phase, first_phase, fixed_frequency).fid()
+
+    def adjust_frequency(self, frequency_shift):
+        """
+        Adjust the centre frequency of the signal.
+
+        Refer to suspect.adjust_frequency for full documentation.
+
+        Parameters
+        ----------
+        frequency_shift: float
+            The amount to shift the frequency, in Hertz.
+
+        Returns
+        -------
+        out : MRSData
+            Frequency adjusted FID
+
+        See Also
+        --------
+        suspect.adjust_frequency : equivalent function
+        """
+        correction = numpy.exp(2j * numpy.pi * (frequency_shift * self.time_axis()))
+        return self.inherit(numpy.multiply(self, correction))
+
+
+class MRSSpectrum(MRSBase):
+    """
+    MRS data in the frequency domain
+    """
+
+    def fid(self):
+        """
+        Returns
+        -------
+        MRSData
+            The inverse-Fourier-shifted and inverse-Fourier-transformed data, represented as a FID
+        """
+        fid = self.inherit(numpy.fft.ifft(numpy.fft.ifftshift(self, axes=-1), axis=-1)).view(MRSData)
+        return fid
+
+    def adjust_phase(self, zero_phase, first_phase=0., fixed_frequency=0.):
+        """
+        Adjust the phases of the signal.
+
+        Refer to suspect.adjust_phase for full documentation.
+
+        Parameters
+        ----------
+        zero_phase: float
+            The zero order phase shift in radians
+        first_phase: float
+            The first order phase shift in radians per Hertz
+        fixed_frequency: float
+            The frequency at which the first order phase shift is zero
+
+        Returns
+        -------
+        out : MRSSpectrum
+            Phase adjusted spectrum
+
+        See Also
+        --------
+        suspect.adjust_phase : equivalent function
+        """
         phase_ramp = numpy.linspace(-self.sw / 2,
                                     self.sw / 2,
                                     self.np,
                                     endpoint=False)
         phase_shift = zero_phase + first_phase * (fixed_frequency + phase_ramp)
-        phased_spectrum = spectrum * numpy.exp(1j * phase_shift)
-        return self.inherit(numpy.fft.ifft(numpy.fft.ifftshift(phased_spectrum, axes=-1), axis=-1))
+        phased_spectrum = self * numpy.exp(1j * phase_shift)
+        return phased_spectrum
 
     def adjust_frequency(self, frequency_shift):
         """
@@ -283,12 +367,11 @@ class MRSData(numpy.ndarray):
 
         Returns
         -------
-        out : MRSData
-            Frequency adjusted FID
+        out : MRSSpectrum
+            Frequency adjusted spectrum
 
         See Also
         --------
         suspect.adjust_frequency : equivalent function
         """
-        correction = numpy.exp(2j * numpy.pi * (frequency_shift * self.time_axis()))
-        return self.inherit(numpy.multiply(self, correction))
+        return self.fid().adjust_frequency(frequency_shift).spectrum()


Test patch:
diff --git a/tests/test_mrs/test_mrsdata.py b/tests/test_mrs/test_mrsdata.py
index 49c19bf..cb34c84 100644
--- a/tests/test_mrs/test_mrsdata.py
+++ b/tests/test_mrs/test_mrsdata.py
@@ -31,3 +31,12 @@ def test_zero_rank_to_scalar():
     data = suspect.MRSData(numpy.ones(1024, 'complex'), 5e-4, 123)
     sum = numpy.sum(data)
     assert numpy.isscalar(sum)
+
+
+def test_spectrum_2_fid():
+    data = suspect.MRSData(numpy.ones(1024, 'complex'), 5e-4, 123)
+    spectrum = data.spectrum()
+    assert type(spectrum) == suspect.MRSSpectrum
+    numpy.testing.assert_equal(spectrum, numpy.fft.fftshift(numpy.fft.fft(data)))
+    fid = spectrum.fid()
+    numpy.testing.assert_equal(data, fid)