diff --git a/requirements.txt b/requirements.txt
index 4ae0d82..31df5ab 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -2,4 +2,5 @@ nose
 python-dotenv
 pytest
 requests
+numpy
 
diff --git a/setup.py b/setup.py
index 242b9a2..1b9ff1e 100644
--- a/setup.py
+++ b/setup.py
@@ -6,11 +6,12 @@
 
 setuptools.setup(
     name='pysofar',
-    version='0.1.10',
+    version='0.1.11',
     license='Apache 2 Licesnse',
     install_requires=[
         'requests',
-        'python-dotenv'
+        'python-dotenv',
+        'numpy'
     ],
     description='Python client for interfacing with the Sofar Wavefleet API to access Spotter Data',
     long_description=readme_contents,
diff --git a/src/pysofar/spotter.py b/src/pysofar/spotter.py
index 62e45e7..febf7de 100644
--- a/src/pysofar/spotter.py
+++ b/src/pysofar/spotter.py
@@ -9,7 +9,9 @@
 Authors: Mike Sosa
 """
 from pysofar.sofar import SofarApi, WaveDataQuery
-
+from pysofar.wavespectra.spectrum1D import WaveSpectrum1D, WaveSpectrum1DInput
+from pysofar.wavefleet_exceptions import ExceptionNoFrequencyData
+import typing
 
 # --------------------- Devices ----------------------------------------------#
 class Spotter:
@@ -279,3 +281,42 @@ def grab_data(self, limit: int = 20,
         _data = _query.execute()
 
         return _data
+
+    def spectra(
+            self,
+            limit: int = 20,
+            start_date: str = None, end_date: str = None,
+            ) -> typing.List[WaveSpectrum1D]:
+        """
+        Grabs the requested spectra for this spotter based on the given keyword arguments
+
+        :param limit: The limit for data to grab. Defaults to 20, For frequency data max of 100 samples at a time.
+        :param start_date: ISO 8601 formatted date string. If not included defaults to beginning of spotters history
+        :param end_date: ISO 8601 formatted date string. If not included defaults to end of spotter history
+
+        :return: Data as a FrequencyDataList Object
+        """
+
+        json_data = self.grab_data(
+            limit=limit,
+            start_date=start_date,
+            end_date=end_date,
+            include_frequency_data=True,
+            include_directional_moments=True)
+
+        if not json_data['frequencyData']:
+            raise ExceptionNoFrequencyData(
+                f'Spotter {self.id} has no spectral data for the requested time range')
+
+        out = []
+        for spectrum in json_data['frequencyData']:
+            # Load the data into the input object- this is a one-to-one
+            # mapping of keys between dictionaries (we would not _need_ to
+            # map onto a input object-> we could pass spectrum directly to
+            # the constructor. But it is clearer to fail here if something
+            # is amiss with the api JSON response format)
+            wave_spectrum_input = WaveSpectrum1DInput(**spectrum)
+
+            out.append(WaveSpectrum1D(wave_spectrum_input))
+
+        return out
\ No newline at end of file
diff --git a/src/pysofar/wavefleet_exceptions.py b/src/pysofar/wavefleet_exceptions.py
index 8d02bc3..ea5c1c4 100644
--- a/src/pysofar/wavefleet_exceptions.py
+++ b/src/pysofar/wavefleet_exceptions.py
@@ -21,4 +21,10 @@ class CouldNotRetrieveFile(Exception):
     """
     Query raised when a requested datafile could not be retrieved.
     """
+    pass
+
+class ExceptionNoFrequencyData(Exception):
+    """
+    Query raised when no frequency data is available for the spotter
+    """
     pass
\ No newline at end of file
diff --git a/src/pysofar/wavespectra/__init__.py b/src/pysofar/wavespectra/__init__.py
new file mode 100644
index 0000000..acb9bd5
--- /dev/null
+++ b/src/pysofar/wavespectra/__init__.py
@@ -0,0 +1,2 @@
+from .spectrum2D import WaveSpectrum2D
+from .spectrum1D import WaveSpectrum1D
\ No newline at end of file
diff --git a/src/pysofar/wavespectra/mem.py b/src/pysofar/wavespectra/mem.py
new file mode 100644
index 0000000..7c2ed24
--- /dev/null
+++ b/src/pysofar/wavespectra/mem.py
@@ -0,0 +1,53 @@
+import numpy
+
+def mem(directions_radians: numpy.ndarray, a1, b1, a2, b2) -> numpy.ndarray:
+    """
+    This function uses the maximum entropy method by Lygre and Krogstadt (1986,JPO)
+    to estimate the directional shape of the spectrum.
+
+    Lygre, A., & Krogstad, H. E. (1986). Maximum entropy estimation of the directional
+    distribution in ocean wave spectra. Journal of Physical Oceanography, 16(12), 2052-2060.
+
+    :param direction_radians: radian directions (going to, anti-clockswise from east) want to
+                              evaluate the spectrum on
+    :param a1: 1st cosine Fourier coefficient of the directional distribution
+    :param b1: 1st sine Fourier coefficient of the directional distribution
+    :param a2: 2nd cosine Fourier coefficient of the directional distribution
+    :param b2: 2nd sine Fourier coefficient of the directional distribution
+
+    :return: Directional distribution as a numpy array
+
+    Note that:
+    d1 = a1; d2 =b1; d3 = a2 and d4=b2 in the defining equations 10.
+    """
+
+    # Ensure that these are numpy arrays
+    a1 = numpy.atleast_1d(a1)
+    b1 = numpy.atleast_1d(b1)
+    a2 = numpy.atleast_1d(a2)
+    b2 = numpy.atleast_1d(b2)
+
+    number_of_directions = len(directions_radians)
+
+    c1 = a1 + 1j * b1
+    c2 = a2 + 1j * b2
+    #
+    # Eq. 13 L&K86
+    #
+    Phi1 = (c1 - c2 * numpy.conj(c1)) / (1 - c1 * numpy.conj(c1))
+    Phi2 = c2 - Phi1 * c1
+    #
+    e1 = numpy.exp(-directions_radians * 1j)
+    e2 = numpy.exp(-directions_radians * 2j)
+
+    numerator = (1 - Phi1 * numpy.conj(c1) - Phi2 * numpy.conj(c2))
+    denominator = numpy.abs(1 - Phi1[:,None] * e1[None,:]
+                              - Phi2[:,None] * e2[None,:]) ** 2
+
+    D = numpy.real( numerator[:,None] / denominator  ) / numpy.pi / 2
+
+    # Normalize to 1. in discrete sense
+    integralApprox = numpy.sum(D,axis=-1) * numpy.pi * 2. / number_of_directions
+    D = D / integralApprox[:,None]
+
+    return numpy.squeeze(D)
\ No newline at end of file
diff --git a/src/pysofar/wavespectra/spectrum1D.py b/src/pysofar/wavespectra/spectrum1D.py
new file mode 100644
index 0000000..ef952aa
--- /dev/null
+++ b/src/pysofar/wavespectra/spectrum1D.py
@@ -0,0 +1,83 @@
+import numpy
+from .spectrum2D import WaveSpectrum2D, WaveSpectrum2DInput
+from .wavespectrum import WaveSpectrum, WaveSpectrumInput
+from .tools import to_datetime, datetime_to_iso_time_string
+from .mem import mem
+from typing import List, Union
+
+class WaveSpectrum1DInput(WaveSpectrumInput):
+    a1: Union[List[float], numpy.ndarray]
+    b1: Union[List[float], numpy.ndarray]
+    a2: Union[List[float], numpy.ndarray]
+    b2: Union[List[float], numpy.ndarray]
+
+class WaveSpectrum1D(WaveSpectrum):
+    spectral_density_units = 'm**2/Hertz'
+
+    def __init__(self,
+                 wave_spectrum1D_input:WaveSpectrum1DInput
+                 ):
+        super().__init__(wave_spectrum1D_input)
+        self._a1 = numpy.array(wave_spectrum1D_input['a1'])
+        self._b1 = numpy.array(wave_spectrum1D_input['b1'])
+        self._b2 = numpy.array(wave_spectrum1D_input['b2'])
+        self._a2 = numpy.array(wave_spectrum1D_input['a2'])
+        self._e = numpy.array(wave_spectrum1D_input['varianceDensity'])
+
+    def frequency_moment(self, power: int, fmin=0, fmax=numpy.inf) -> float:
+        range = self._range(fmin,fmax)
+
+        return numpy.trapz(
+            self.variance_density[range] * self.frequency[range] ** power,
+            self.frequency[range])
+
+    def _create_wave_spectrum_input(self)->WaveSpectrum1DInput:
+        return WaveSpectrum1DInput(
+            frequency=list(self.frequency),
+            varianceDensity=list(self.variance_density),
+            timestamp=datetime_to_iso_time_string(self.timestamp),
+            latitude=self.latitude,
+            longitude=self.longitude,
+            a1=list(self.a1),
+            b1=list(self.b1),
+            a2=list(self.a2),
+            b2=list(self.b2)
+            )
+
+    def spectrum2d(self, number_of_directions: int,
+                   method:str='maximum_entropy_method')->WaveSpectrum2D:
+        """
+        Construct a 2D spectrum based on the 1.5D spectrum and a spectral
+        reconstruction method.
+
+        :param number_of_directions: length of the directional vector for the
+        2D spectrum. Directions returned are in degrees
+        """
+        direction = numpy.linspace(0, 360, number_of_directions,
+                                   endpoint=False)
+
+        # Jacobian to transform distribution as function of radian angles into
+        # degrees.
+        Jacobian = numpy.pi / 180
+
+        if method.lower() in ['maximum_entropy_method', 'mem']:
+            # reconstruct the directional distribution using the maximum entropy
+            # method.
+            D = mem(direction * numpy.pi / 180, self.a1, self.b1, self.a2,
+                    self.b2) * Jacobian
+        else:
+            raise Exception(f'unsupported spectral estimator method: {method}')
+
+        wave_spectrum2D_input = WaveSpectrum2DInput(
+            frequency=self.frequency,
+            directions=direction,
+            varianceDensity=self.variance_density[:, None] * D,
+            timestamp=self.timestamp,
+            longitude=self.longitude,
+            latitude=self.latitude
+        )
+
+        # We return a 2D wave spectrum object.
+        return WaveSpectrum2D(wave_spectrum2D_input)
+
+
diff --git a/src/pysofar/wavespectra/spectrum2D.py b/src/pysofar/wavespectra/spectrum2D.py
new file mode 100644
index 0000000..0cf2015
--- /dev/null
+++ b/src/pysofar/wavespectra/spectrum2D.py
@@ -0,0 +1,63 @@
+import numpy
+from .wavespectrum import WaveSpectrum, WaveSpectrumInput
+from typing import List, Union
+from .tools import datetime_to_iso_time_string
+
+
+class WaveSpectrum2DInput(WaveSpectrumInput):
+    directions: Union[List[float], numpy.ndarray]
+
+class WaveSpectrum2D(WaveSpectrum):
+    def __init__(self,
+                 wave_spectrum2D_input:WaveSpectrum2DInput
+                 ):
+
+        super().__init__(wave_spectrum2D_input)
+        self.direction = wave_spectrum2D_input['directions']
+        self._a1 = self._directional_moment('a', 1, normalized=True)
+        self._b1 = self._directional_moment('b', 1, normalized=True)
+        self._a2 = self._directional_moment('a', 2, normalized=True)
+        self._b2 = self._directional_moment('b', 2, normalized=True)
+        self._e = self._directional_moment('zero', 0, normalized=False)
+
+    def _delta(self):
+        angles = self.direction
+        forward_diff = (numpy.diff(angles, append=angles[0]) + 180) % 360 - 180
+        backward_diff = (numpy.diff(angles,
+                                    prepend=angles[-1]) + 180) % 360 - 180
+        return (forward_diff + backward_diff) / 2
+
+    def _directional_moment(self, kind='zero', order=0,
+                            normalized=True) -> numpy.array:
+        delta = self._delta()
+        if kind == 'a':
+            harmonic = numpy.cos(self.radian_direction * order) * delta
+        elif kind == 'b':
+            harmonic = numpy.sin(self.radian_direction * order) * delta
+        elif kind == 'zero':
+            harmonic = delta
+        else:
+            raise Exception('Unknown moment')
+        values = numpy.sum(self.variance_density * harmonic[None, :], axis=-1)
+
+        if normalized:
+            scale = numpy.sum(self.variance_density * delta[None, :], axis=-1)
+        else:
+            scale = 1
+
+        return values / scale
+
+    def frequency_moment(self, power: int, fmin=0, fmax=numpy.inf) -> float:
+        range = self._range(fmin, fmax)
+        return numpy.trapz(self.e[range] * self.frequency[range] ** power,
+                           self.frequency[range])
+
+    def _create_wave_spectrum_input(self)->WaveSpectrum2DInput:
+        return WaveSpectrum2DInput(
+            frequency=list(self.frequency),
+            directions=list(self.direction),
+            varianceDensity=list(self.variance_density),
+            timestamp=datetime_to_iso_time_string(self.timestamp),
+            latitude=self.latitude,
+            longitude=self.longitude,
+        )
\ No newline at end of file
diff --git a/src/pysofar/wavespectra/tools.py b/src/pysofar/wavespectra/tools.py
new file mode 100644
index 0000000..bbc1908
--- /dev/null
+++ b/src/pysofar/wavespectra/tools.py
@@ -0,0 +1,23 @@
+from datetime import datetime, timezone
+import typing
+
+def datetime_to_iso_time_string(time: datetime):
+    return time.strftime("%Y-%m-%dT%H:%M:%S.%fZ")
+
+def to_datetime(time: typing.Union[float, int, datetime, str]):
+    if time is None:
+        return None
+
+    if isinstance(time, datetime):
+        return time.astimezone(timezone.utc)
+    elif isinstance(time, str):
+        if time[-1] == "Z":
+            # From isoformat does not parse "Z" as a valid timezone designator
+            time = time[:-1] + "+00:00"
+        try:
+            return datetime.fromisoformat(time)
+        except ValueError as e:
+            return datetime.strptime(time, "%Y-%m-%dT%H:%M:%S.%f%z")
+    else:
+        return datetime.fromtimestamp(time, tz=timezone.utc)
+
diff --git a/src/pysofar/wavespectra/wavespectrum.py b/src/pysofar/wavespectra/wavespectrum.py
new file mode 100644
index 0000000..f9ebd47
--- /dev/null
+++ b/src/pysofar/wavespectra/wavespectrum.py
@@ -0,0 +1,177 @@
+import numpy
+from pysofar.wavespectra.tools import to_datetime, datetime_to_iso_time_string
+from typing import TypedDict, List
+
+class WaveSpectrumInput(TypedDict):
+    frequency: List[float]
+    varianceDensity: List
+    timestamp: str
+    latitude: float
+    longitude: float
+
+class WaveSpectrum():
+    """
+    Base spectral class.
+    """
+    frequency_units = 'Hertz'
+    angular_units = 'Degrees'
+    spectral_density_units = 'm**2/Hertz'
+    angular_convention = 'Wave travel direction (going-to), measured anti-clockwise from East'
+
+    def __init__(self,
+                 wave_spectrum_input: WaveSpectrumInput
+                 ):
+
+        self.frequency = numpy.array(wave_spectrum_input['frequency'])
+        self.variance_density = numpy.array(wave_spectrum_input['varianceDensity'])
+        self.direction = None
+        self._a1 = None
+        self._b1 = None
+        self._a2 = None
+        self._b2 = None
+        self._e = None
+        self.timestamp = to_datetime(wave_spectrum_input['timestamp'])
+        self.longitude = wave_spectrum_input['longitude']
+        self.latitude = wave_spectrum_input['latitude']
+
+    def frequency_moment(self, power: int, fmin=0, fmax=numpy.inf) -> float:
+        pass
+
+    def _create_wave_spectrum_input(self)->WaveSpectrumInput:
+        return WaveSpectrumInput(
+            frequency=list(self.frequency),
+            varianceDensity=list(self.variance_density),
+            timestamp=datetime_to_iso_time_string(self.timestamp),
+            latitude=self.latitude,
+            longitude=self.longitude
+        )
+
+    def _range(self, fmin=0.0, fmax=numpy.inf)->numpy.ndarray:
+        return (self.frequency >= fmin) & (self.frequency < fmax)
+
+    @property
+    def radian_direction(self) -> numpy.ndarray:
+        return self.direction * numpy.pi / 180
+
+    @property
+    def radian_frequency(self) -> numpy.ndarray:
+        return self.frequency * 2 * numpy.pi
+
+    @property
+    def e(self) -> numpy.array:
+        return self._e
+
+    @property
+    def a1(self) -> numpy.array:
+        return self._a1
+
+    @property
+    def b1(self) -> numpy.array:
+        return self._b1
+
+    @property
+    def a2(self) -> numpy.array:
+        return self._a2
+
+    @property
+    def b2(self) -> numpy.array:
+        return self._b2
+
+    @property
+    def A1(self) -> numpy.array:
+        return self.a1 * self.e
+
+    @property
+    def B1(self) -> numpy.array:
+        return self.b1 * self.e
+
+    @property
+    def A2(self) -> numpy.array:
+        return self.a2 * self.e
+
+    @property
+    def B2(self) -> numpy.array:
+        return self.b2 * self.e
+
+    def m0(self, fmin=0, fmax=numpy.inf) -> float:
+        return self.frequency_moment(0, fmin, fmax)
+
+    def m1(self, fmin=0, fmax=numpy.inf) -> float:
+        return self.frequency_moment(1, fmin, fmax)
+
+    def m2(self, fmin=0, fmax=numpy.inf) -> float:
+        return self.frequency_moment(2, fmin, fmax)
+
+    def hm0(self, fmin=0, fmax=numpy.inf) -> float:
+        return 4 * numpy.sqrt(self.m0(fmin, fmax))
+
+    def tm01(self, fmin=0, fmax=numpy.inf) -> float:
+        return self.m0(fmin, fmax) / self.m1(fmin, fmax)
+
+    def tm02(self, fmin=0, fmax=numpy.inf) -> float:
+        return numpy.sqrt(self.m0(fmin, fmax) / self.m2(fmin, fmax))
+
+    def peak_index(self, fmin=0, fmax=numpy.inf) -> float:
+        range = self._range(fmin, fmax)
+
+        return numpy.argmax(self.e[range])
+
+    def peak_frequency(self, fmin=0, fmax=numpy.inf) -> float:
+        return self.frequency[self.peak_index(fmin, fmax)]
+
+    def peak_period(self, fmin=0, fmax=numpy.inf) -> float:
+        return 1 / self.peak_frequency(fmin, fmax)
+
+    def peak_direction(self, fmin=0, fmax=numpy.inf):
+        index = self.peak_index(fmin, fmax)
+        a1 = self.a1[index]
+        b1 = self.b1[index]
+        return self._mean_direction(a1, b1)
+
+    def peak_spread(self, fmin=0, fmax=numpy.inf):
+        index = self.peak_index(fmin, fmax)
+        a1 = self.a1[index]
+        b1 = self.b1[index]
+        return self._spread(a1, b1)
+
+    @staticmethod
+    def _mean_direction(a1, b1):
+        return numpy.arctan2(b1, a1) * 180 / numpy.pi
+
+    @staticmethod
+    def _spread(a1, b1):
+        return numpy.sqrt(
+            2 - 2 * numpy.sqrt(a1 ** 2 + b1 ** 2)) * 180 / numpy.pi
+
+    @property
+    def mean_direction(self):
+        return self._mean_direction(self.a1, self.b1)
+
+    @property
+    def mean_spread(self):
+        return self._spread(self.a1, self.b1)
+
+    def _spectral_weighted(self, property, fmin=0, fmax=numpy.inf):
+        range = (self._range(fmin,fmax)) & numpy.isfinite( property )
+
+        return numpy.trapz(property[range] * self.e[range],
+                           self.frequency[range]) / self.m0(fmin, fmax)
+
+    def bulk_direction(self, fmin=0, fmax=numpy.inf):
+        return self._mean_direction(self.bulk_a1(fmin, fmax),
+                                    self.bulk_b1(fmin, fmax))
+
+    def bulk_spread(self, fmin=0, fmax=numpy.inf):
+        return self._spread(self.bulk_a1(fmin, fmax), self.bulk_b1(fmin, fmax))
+
+    def bulk_a1(self, fmin=0, fmax=numpy.inf):
+        return self._spectral_weighted(self.a1, fmin, fmax)
+
+    def bulk_b1(self, fmin=0, fmax=numpy.inf):
+        return self._spectral_weighted(self.b1, fmin, fmax)
+
+    def bulk_a2(self, fmin=0, fmax=numpy.inf):
+        return self._spectral_weighted(self.a2, fmin, fmax)
+
+    def bulk_b2(self, fmin=0, fmax=numpy.inf):
+        return self._spectral_weighted(self.b2, fmin, fmax)
\ No newline at end of file