Merge pull request #6 from myokit/spectral

Added two spectral analysis methods

Author: MichaelClerx

README.md

@@ -3,7 +3,34 @@
 
 # Datkit
 
-This module contains some simple methods that are useful when analysing time
+This module contains some simple methods that are useful when analysing time 
 series data.
 
-The code is tested on a recent version of Ubuntu & Python 3, but is so simple that it should work everywhere else too.
+The code is tested on a recent version of Ubuntu & Python 3, but is so simple
+that it should work everywhere else too.
+
+## Installation
+
+To install the latest release from PyPI, use
+
+```
+pip install datkit
+```
+
+## Installation for development
+
+To install from the repo, use e.g.
+```
+python setup.py install -e .
+```
+
+Tests can then be run with
+```
+python -m unittest
+```
+
+And docs can be built with
+```
+cd docs
+make clean html
+```

datkit/__init__.py

@@ -72,3 +72,9 @@
     moving_average,
     window_size,
 )
+
+from ._spectral import (  # noqa
+    amplitude_spectrum,
+    power_spectral_density,
+)
+

datkit/_spectral.py

@@ -0,0 +1,80 @@
+#
+# Methods to inspect time series in the frequency domain.
+# Note that much more advanced methods are available e.g. in scipy.
+#
+# This file is part of Datkit.
+# For copyright, sharing, and licensing, see https://github.com/myokit/datkit/
+#
+import numpy as np
+
+import datkit
+
+
+def amplitude_spectrum(times, values):
+    """
+    Calculates the amplitude spectrum of a regularly spaced time series
+    ``(times, current)``, returning a tuple ``(frequency, amplitude)``.
+
+    Example::
+
+        t = np.linspace(0, 10, 1000)
+        v = 3 * np.sin(t * (2 * np.pi * 2)) + 10 * np.cos(t * (2 * np.pi * 3))
+        f, a = amplitude_spectrum(t, v)
+
+        fig = plt.figure()
+        ax = fig.add_subplot(2, 1, 1)
+        ax.plot(t, v)
+        ax = fig.add_subplot(2, 1, 2)
+        ax.set_xlim(0, 10)
+        ax.stem(f, a)
+        plt.show()
+
+    """
+    dt = datkit.sampling_interval(times)
+    n = len(times)
+    # Select only the positive frequencies
+    m = n // 2
+    # FFT, normalised by number of point
+    a = np.absolute(np.fft.fft(values)[:m]) / n
+    # Only using one half, so multiply values of mirrored points by 2
+    a[1:-1] *= 2
+    # Get corresponding frequencies and return
+    f = np.fft.fftfreq(n, dt)[:m]
+    return f, a
+
+
+def power_spectral_density(times, values):
+    """
+    Estimates the power spectral density of a regularly spaced time series
+    ``(times, current)``, returning a tuple ``(frequency, density)``.
+
+    For times in units "T" and values in units "V", the returned frequencies
+    will be in units "F=1/T" and the densities in "V^2/F".
+
+    For a less noisy version, use ``scipy.signal.welch(values,
+    1 / datkit.sampling_interval(times))`` instead.
+
+    Example::
+
+        t = np.linspace(0, 10, 1000)
+        v = 5 * np.sin(t * (2 * np.pi * 2)) + 10 * np.cos(t * (2 * np.pi * 3))
+        f, psd = power_spectral_density(t, v)
+
+        fig = plt.figure()
+        ax = fig.add_subplot(2, 1, 1)
+        ax.plot(t, v)
+        ax = fig.add_subplot(2, 1, 2)
+        ax.set_yscale('log')
+        ax.set_xlim(0, 10)
+        ax.plot(f, psd)
+        plt.show()
+
+    """
+    dt = datkit.sampling_interval(times)
+    n = len(times)
+    m = n // 2
+    p = np.absolute(np.fft.fft(values)[:m])**2 * dt / n
+    p[1:-1] *= 2
+    f = np.fft.fftfreq(n, dt)[:m]
+    return f, p
+

datkit/tests/test_spectral.py

@@ -0,0 +1,61 @@
+#!/usr/bin/env python3
+#
+# Tests the datkit spectral analysis methods.
+#
+# This file is part of Datkit.
+# For copyright, sharing, and licensing, see https://github.com/myokit/datkit/
+#
+import unittest
+
+import numpy as np
+
+import datkit as d
+
+
+class SpectralTest(unittest.TestCase):
+    """ Tests methods from the hidden _spectral module. """
+
+    def test_amplitude_spectrum(self):
+
+        t = np.linspace(0, 10, 1000)
+        v = 3 * np.sin(t * (2 * np.pi * 3)) + 10 * np.cos(t * (2 * np.pi * 7))
+        f, a = d.amplitude_spectrum(t, v)
+        self.assertAlmostEqual(f[30], 2.997)
+        self.assertAlmostEqual(f[70], 6.993)
+        self.assertLess(a[29], 0.1)
+        self.assertAlmostEqual(a[30], 3, 0)
+        self.assertLess(a[31], 0.1)
+        self.assertLess(a[69], 1)
+        self.assertAlmostEqual(a[70], 10, 0)
+        self.assertLess(a[71], 1)
+
+        t = np.linspace(0, 10, 1235)
+        v = 6 * np.sin(t * (2 * np.pi * 2)) + 4 * np.cos(t * (2 * np.pi * 5))
+        f, a = d.amplitude_spectrum(t, v)
+        self.assertAlmostEqual(f[20], 1.998, 3)
+        self.assertAlmostEqual(f[50], 4.996, 3)
+        self.assertLess(a[19], 0.11)
+        self.assertAlmostEqual(a[20], 6, 0)
+        self.assertLess(a[21], 0.11)
+        self.assertLess(a[49], 0.2)
+        self.assertAlmostEqual(a[50], 4, 0)
+        self.assertLess(a[51], 0.2)
+
+    def test_power_spectral_density(self):
+
+        t = np.linspace(0, 10, 1000)
+        v = 3 * np.sin(t * (2 * np.pi)) + 7 * np.cos(t * (2 * np.pi * 3))
+        f, psd = d.power_spectral_density(t, v)
+        self.assertAlmostEqual(f[10], 0.999, 3)
+        self.assertAlmostEqual(f[30], 2.997, 3)
+
+        self.assertLess(psd[9], 0.01)
+        self.assertAlmostEqual(psd[10], 45, 1)
+        self.assertLess(psd[11], 0.01)
+        self.assertLess(psd[29], 0.3)
+        self.assertAlmostEqual(psd[30], 245, 0)
+        self.assertLess(psd[31], 0.3)
+
+
+if __name__ == '__main__':
+    unittest.main()

docs/source/checking_time_vectors.rst

@@ -2,6 +2,9 @@
 Checking time vectors
 *********************
 
+The methods below can be used to test if sequences (typically numpy arrays)
+contain regularly sampled times.
+
 .. currentmodule:: datkit
 
 .. autofunction:: is_increasing

docs/source/finding_points.rst

@@ -2,6 +2,10 @@
 Finding points
 **************
 
+The methods listed here let you find indices corresponding to points in time,
+or perform aggregate functions (e.g. mean, min, max) on intervals within a time
+series.
+
 .. currentmodule:: datkit
 
 Indices and values at specific times

docs/source/index.rst

@@ -10,7 +10,8 @@ series data.
 ..  toctree::
 
     self
+    checking_time_vectors
     finding_points
     smoothing
-    checking_time_vectors
+    spectral_analysis
 

docs/source/smoothing.rst

@@ -2,6 +2,8 @@
 Smoothing
 *********
 
+The methods below can perform basic smoothing of noisy time series.
+
 .. currentmodule:: datkit
 
 .. autofunction:: moving_average

docs/source/spectral_analysis.rst

@@ -0,0 +1,14 @@
+*****************
+Spectral analysis
+*****************
+
+The methods listed below can be used to perform very rudimentary spectral (or
+frequency domain) analysis. They are most useful to make quick plots of the
+frequency content of a time series.
+
+.. currentmodule:: datkit
+
+.. autofunction:: amplitude_spectrum
+
+.. autofunction:: power_spectral_density
+