Repository: StingraySoftware/stingray
Pull request #208 (2017-07-12 17:18:21)

Bispectrum
Implementation of `Bispectrum` as a fourier transform of 3rd Order Cumulant (Triple AutoCorrelation). 
Results are verified from a matlab version of bispectrum code.

Further Addition(s):
1. lag window functions (Currently no window is applied) # #


Patch:
diff --git a/stingray/bispectrum.py b/stingray/bispectrum.py
new file mode 100644
index 00000000..06e5280c
--- /dev/null
+++ b/stingray/bispectrum.py
@@ -0,0 +1,374 @@
+from __future__ import division
+
+import numpy as np
+from scipy.linalg import toeplitz
+from scipy.fftpack import fftshift, fft2, ifftshift, fftfreq
+
+from stingray import lightcurve
+from stingray.utils import simon
+import matplotlib.pyplot as plt
+
+
+class Bispectrum(object):
+    def __init__(self, lc, maxlag=None, scale='biased'):
+        """
+                Makes a :class:`Bispectrum` object from a given :class:`Lightcurve`.
+                
+                Bispectrum is a higher order time series analysis method and is calculated by indirect method as
+                fourier transform of triple auto-correlation function also called as 3rd Order cumulant.
+                
+                Parameters
+                ----------
+                lc: lightcurve.Lightcurve object
+                    The light curve data for bispectrum calculation.
+                maxlag: int, optional, default None
+                    Maximum lag on both positive and negative sides of 
+                    3rd order cumulant (Similar to lags in correlation).
+                    if None, max lag is set to one-half of length of lightcurve.
+                scale: {'biased', 'unbiased'}, optional, default 'biased'
+                    Flag to decide biased or unbiased normalization for 3rd order cumulant function.
+
+                Attributes
+                ----------
+                lc: lightcurve.Lightcurve 
+                    The light curve data for bispectrum.
+                fs: float
+                    Sampling freq of light curve.        
+                n: int
+                    Total Number of samples of light curve observations.
+                maxlag: int
+                    Maximum lag on both positive and negative sides of 
+                    3rd order cumulant (Similar to lags in correlation)
+                signal: numpy.ndarray
+                    Row vector of lightcurve counts for matrix operations
+                scale: {'biased', 'unbiased'}
+                    Flag to decide biased or unbiased normalization for 3rd order cumulant function. 
+                lags: numpy.ndarray
+                    An array of time lags for which 3rd order cumulant is calculated
+                freq: numpy.ndarray
+                    An array of freq values for bispectrum.
+                cum3: numpy.ndarray
+                    A maxlag*2+1 x maxlag*2+1 matrix containing 3rd order cumulant data for different lags.
+                bispec: numpy.ndarray
+                    A maxlag*2+1 x maxlag*2+1 matrix containing bispectrum data for different frequencies.
+                bispec_mag: numpy.ndarray
+                    Magnitude of Bispectrum
+                bispec_phase: numpy.ndarray
+                    Phase of Bispectrum
+                References
+                ----------     
+                [1] The biphase explained: understanding the asymmetries invcoupled Fourier components of astronomical timeseries
+                by Thomas J. Maccarone Department of Physics, Box 41051, Science Building, Texas Tech University, Lubbock TX 79409-1051
+                School of Physics and Astronomy, University of Southampton, SO16 4ES
+                    
+                [2] T. S. Rao, M. M. Gabr, An Introduction to Bispectral Analysis and Bilinear Time
+                Series Models, Lecture Notes in Statistics, Volume 24, D. Brillinger, S. Fienberg,
+                J. Gani, J. Hartigan, K. Krickeberg, Editors, Springer-Verlag, New York, NY, 1984.
+                
+                [3] Matlab version of bispectrum under following link. 
+                https://www.mathworks.com/matlabcentral/fileexchange/60-bisp3cum
+                
+                Example
+                -------
+                >> from stingray.lightcurve import Lightcurve
+                >> from stingray.bispectrum import Bispectrum
+                >> lc = Lightcurve([1,2,3,4,5],[2,3,1,1,2])
+                >> bs = Bispectrum(lc,maxlag=1)
+                >> bs.lags
+                array([-1.,  0.,  1.])
+                >> bs.freq
+                array([-0.5,  0. ,  0.5])
+                >> bs.cum3
+                array([[-0.2976,  0.1024,  0.1408],
+                    [ 0.1024,  0.144 , -0.2976],
+                    [ 0.1408, -0.2976,  0.1024]])
+                >> bs.bispec_mag
+                array([[ 1.26336794,  0.0032    ,  0.0032    ],
+                    [ 0.0032    ,  0.16      ,  0.0032    ],
+                    [ 0.0032    ,  0.0032    ,  1.26336794]])
+                >> bs.bispec_phase
+                array([[ -9.65946229e-01,   2.25347190e-14,   3.46944695e-14],
+                    [  0.00000000e+00,   3.14159265e+00,   0.00000000e+00],
+                    [ -3.46944695e-14,  -2.25347190e-14,   9.65946229e-01]])                  
+        """
+
+        # Function call to create Bispectrum Object
+        self._make_bispectrum(lc, maxlag, scale)
+
+    def _make_bispectrum(self, lc, maxlag, scale):
+        """
+            Makes a Bispectrum Object with given lighcurve, maxlag and scale. 
+        """
+
+        if not isinstance(lc, lightcurve.Lightcurve):
+            raise TypeError('lc must be a lightcurve.ightcurve object')
+
+        self.lc = lc
+        self.fs = 1 / lc.dt
+        self.n = self.lc.n
+
+        if maxlag is None:
+            # if maxlag is not specified, it is set to half of length of lightcurve
+            self.maxlag = np.int(self.lc.n / 2)
+        else:
+            if not (isinstance(maxlag, int)):
+                raise ValueError('maxlag must be an integer')
+
+            # if negative maxlag is entered, convert it to +ve
+            if maxlag < 0:
+                self.maxlag = -maxlag
+            else:
+                self.maxlag = maxlag
+
+        if isinstance(scale, str) is False:
+            raise TypeError("scale must be a string")
+
+        if scale.lower() not in ["biased", "unbiased"]:
+            raise ValueError("scale can only be either 'biased' or 'unbiased'.")
+        self.scale = scale.lower()
+
+        # Other Atributes
+        self.lags = None
+        self.cum3 = None
+        self.freq = None
+        self.bispec = None
+        self.bispec_mag = None
+        self.bispec_phase = None
+
+        # converting to a row vector to apply matrix operations
+        self.signal = np.reshape(lc, (1, len(self.lc.counts)))
+
+        # Mean subtraction before bispecrum calculation
+        self.signal = self.signal - np.mean(lc.counts)
+
+        self._cumulant3()
+        self._normalize_cumulant3()
+        self._cal_bispec()
+
+    def _cumulant3(self):
+        """
+            Calculates the 3rd Order cummulant of the lightcurve.
+            Assigns: 
+            self.cum3, 
+            self.lags
+        """
+        # Initialize square cumulant matrix if zeros
+        cum3_dim = 2 * self.maxlag + 1
+        self.cum3 = np.zeros((cum3_dim, cum3_dim))
+
+        # calculate lags for different values of 3rd order cumulant
+        lagindex = np.arange(-self.maxlag, self.maxlag + 1)
+        self.lags = lagindex * self.lc.dt
+
+        # Defines indices for matrices
+        ind = np.arange((self.n - self.maxlag) - 1, self.n)
+        ind_t = np.arange(self.maxlag, self.n)
+        zero_maxlag = np.zeros((1, self.maxlag))
+        zero_maxlag_t = zero_maxlag.transpose()
+
+        sig = self.signal.transpose()
+
+        rev_signal = np.array([self.signal[0][::-1]])
+        col = np.concatenate((sig[ind], zero_maxlag_t), axis=0)
+        row = np.concatenate((rev_signal[0][ind_t], zero_maxlag[0]), axis=0)
+
+        # converts row and column into a toeplitz matrix
+        toep = toeplitz(col, row)
+        rev_signal = np.repeat(rev_signal, [2 * self.maxlag + 1], axis=0)
+
+        # Calulates Cummulant of 1D signal i.e. Lightcurve counts
+        self.cum3 = self.cum3 + np.matmul(np.multiply(toep, rev_signal), toep.transpose())
+
+    def _normalize_cumulant3(self):
+        """
+        Scales (biased or ubiased) the 3rd Order cumulant of the lightcurve .
+        Updates: 
+        self.cum3
+        """
+
+        # Biased scaling of cummulant
+        if self.scale == 'biased':
+            self.cum3 = self.cum3 / self.n
+        else:
+            # unbiased Scaling of cummulant
+            maxlag1 = self.maxlag + 1
+
+            # Scaling matrix initialized used to do unbiased normalization of cumulant
+            scal_matrix = np.zeros((maxlag1, maxlag1), dtype='int64')
+
+            # Calculate scaling matrix for unbiased normalization
+            for k in range(maxlag1):
+                maxlag1k = (maxlag1 - (k + 1))
+                scal_matrix[k, k:maxlag1] = np.tile(self.n - maxlag1k, (1, maxlag1k + 1))
+            scal_matrix += np.triu(scal_matrix, k=1).transpose()
+
+            maxlag1ind = np.arange(self.maxlag - 1, -1, -1)
+            lagdiff = self.n - maxlag1
+
+            # Rows and columns for Toeplitz matrix
+            col = np.arange(lagdiff, self.n - 1)
+            col = np.reshape(col, (1, len(col))).transpose()
+            row = np.arange(lagdiff, (self.n - 2 * self.maxlag) - 1, -1)
+            row = np.reshape(row, (1, len(row)))
+
+            # Toeplitz matrix
+            toep_matrix = toeplitz(col, row)
+            # Matrix used to concatenate with scaling matrix
+            conc_mat = np.array([scal_matrix[self.maxlag, maxlag1ind]])
+            join_matrix = np.concatenate((toep_matrix, conc_mat), axis=0)
+            scal_matrix = np.concatenate((scal_matrix, join_matrix), axis=1)
+            co_mat = scal_matrix[maxlag1ind, :]
+            co_mat = co_mat[:, np.arange(2 * self.maxlag, -1, -1)]
+
+            # Scaling matrix calculated
+            scal_matrix = np.concatenate((scal_matrix, co_mat), axis=0)
+            # Set numbers less than 1 to be equal to 1
+            scal_matrix[scal_matrix < 1] = 1
+            self.cum3 = np.divide(self.cum3, scal_matrix)
+
+    def _cal_bispec(self):
+        """
+            Calculates bispectrum as a fourier transform of 3rd Order Cumulant.
+            Assigns: 
+            self.freq
+            self.bispec
+            self.bispec_mag
+            self.bispec_phase
+        """
+        self.freq = (1 / 2) * self.fs * (self.lags / self.lc.dt) / self.maxlag
+        self.bispec = fftshift(fft2(ifftshift(self.cum3)))
+        self.bispec_mag = np.abs(self.bispec)
+        self.bispec_phase = np.angle((self.bispec))
+
+    def plot_cum3(self, axis=None, save=False, filename=None):
+
+        """
+            Plot the 3rd order cumulant as function of time lags using Matplotlib.
+            Plot the self.cum3 on a graph with ``self.lags`` on x-axis and y-axis and
+            ``self.cum3`` on z-axis
+            Parameters
+            ----------
+            axis : list, tuple, string, default None
+                Parameter to set axis properties of Matplotlib figure. For example
+                it can be a list like ``[xmin, xmax, ymin, ymax]`` or any other
+                acceptable argument for `matplotlib.pyplot.axis()` function.
+            save : boolean, optional (default=False)
+                If True, save the figure with specified filename.
+            filename : str
+                File name of the image to save. Depends on the boolean ``save``.
+            Returns
+ +          ---------
+            plt : matplotlib.pyplot object
+                Reference to plot, call show() to display it
+        """
+
+        try:
+            import matplotlib.pyplot as plt
+        except ImportError:
+            raise ImportError("Matplotlib required for plot()")
+
+        cont = plt.contourf(self.lags, self.lags, self.cum3, 100, cmap=plt.cm.Spectral_r)
+        plt.colorbar(cont)
+        plt.title('3rd Order Cumulant')
+        plt.xlabel('lags 1')
+        plt.ylabel('lags 2')
+
+        if axis is not None:
+            plt.axis(axis)
+
+        if save:
+            if filename is None:
+                plt.savefig('bispec_cum3.png')
+            else:
+                plt.savefig(filename)
+        return plt
+
+    def plot_mag(self, axis=None, save=False, filename=None):
+
+        """
+            Plot the Magnitude of Bispectrum as function of freq using Matplotlib.
+            Plot the self.bispec_mag on a graph with ``self.freq`` on x-axis and y-axis and
+            ``self.bispec_mag`` on z-axis
+            Parameters
+            ----------
+            axis : list, tuple, string, default None
+                Parameter to set axis properties of Matplotlib figure. For example
+                it can be a list like ``[xmin, xmax, ymin, ymax]`` or any other
+                acceptable argument for `matplotlib.pyplot.axis()` function.
+            save : boolean, optional (default=False)
+                If True, save the figure with specified filename.
+            filename : str
+                File name of the image to save. Depends on the boolean ``save``.
+            Returns
+            ---------
+            plt : matplotlib.pyplot object
+                Reference to plot, call show() to display it
+        """
+        try:
+            import matplotlib.pyplot as plt
+        except ImportError:
+            raise ImportError("Matplotlib required for plot()")
+
+        cont = plt.contourf(self.freq, self.freq, self.bispec_mag, 100, cmap=plt.cm.Spectral_r)
+        plt.colorbar(cont)
+        plt.title('Bispectrum Magnitude')
+        plt.xlabel('freq 1')
+        plt.ylabel('freq 2')
+
+        if axis is not None:
+            plt.axis(axis)
+
+        if save:
+            if filename is None:
+                plt.savefig('bispec_mag.png')
+            else:
+                plt.savefig(filename)
+        return plt
+
+    def plot_phase(self, axis=None, save=False, filename=None):
+
+        """
+            Plot the Phase of Bispectrum as function of freq using Matplotlib.
+            Plot the self.bispec_phase on a graph with ``self.phase`` on x-axis and y-axis and
+            ``self.bispec_phase`` on z-axis
+            Parameters
+            ----------
+            axis : list, tuple, string, default None
+                Parameter to set axis properties of Matplotlib figure. For example
+                it can be a list like ``[xmin, xmax, ymin, ymax]`` or any other
+                acceptable argument for `matplotlib.pyplot.axis()` function.
+            save : boolean, optional (default=False)
+                If True, save the figure with specified filename.
+            filename : str
+                File name of the image to save. Depends on the boolean ``save``.
+            Returns
+ +          ---------
+            plt : matplotlib.pyplot object
+                Reference to plot, call show() to display it
+        """
+        try:
+            import matplotlib.pyplot as plt
+        except ImportError:
+            raise ImportError("Matplotlib required for plot()")
+
+        cont = plt.contourf(self.freq, self.freq, self.bispec_phase, 100, cmap=plt.cm.Spectral_r)
+        plt.colorbar(cont)
+        plt.title('Bispectrum Phase')
+        plt.xlabel('freq 1')
+        plt.ylabel('freq 2')
+
+        if axis is not None:
+            plt.axis(axis)
+
+        # Save figure
+        if save:
+            if filename is None:
+                plt.savefig('bispec_phase.png')
+            else:
+                plt.savefig(filename)
+        return plt
+
+
+lc = lightcurve.Lightcurve([1, 2, 3, 4, 5], [2, 3, 1, 1, 2])
+bs = Bispectrum(lc)
+print(bs.bispec_mag)


Test patch:
diff --git a/stingray/tests/test_bispectrum.py b/stingray/tests/test_bispectrum.py
new file mode 100644
index 00000000..814cef7f
--- /dev/null
+++ b/stingray/tests/test_bispectrum.py
@@ -0,0 +1,325 @@
+import numpy as np
+
+from astropy.tests.helper import pytest
+import warnings
+import os
+
+from stingray import Lightcurve
+from stingray.bispectrum import Bispectrum
+from stingray.exceptions import StingrayError
+
+try:
+    import matplotlib.pyplot as plt
+    HAS_MPL = True
+except ImportError:
+    HAS_MPL = False
+
+
+class TestBispectrum(object):
+    @classmethod
+    def setup_class(cls):
+        cls.lc = Lightcurve([1, 2, 3, 4, 5], [2, 3, 2, 4, 1])
+        cls.lc1 = Lightcurve([0.5, 1.0, 1.5, 2.0, 2.5, 3.0], [2, 1, 3, 1, 4, 2])
+
+    def test_create_bispectrum(self):
+        """
+        Demonstrate that we can create a Bispectrum object.
+        """
+        Bispectrum(self.lc)
+
+    def test_wrong_lc(self):
+        lc = [1, 2, 3, 4]
+        with pytest.raises(TypeError):
+            Bispectrum(lc)
+
+    def test_wrong_maxlag(self):
+        with pytest.raises(ValueError):
+            Bispectrum(self.lc, maxlag='123')
+
+    def test_maxlag_none(self):
+        bs = Bispectrum(self.lc)
+        assert bs.maxlag == np.int(self.lc.n / 2)
+
+    def test_neg_maxlag(self):
+        bs = Bispectrum(self.lc, maxlag=-2)
+        assert bs.maxlag == 2
+
+    def test_bispectrum_with_none_maxlag(self):
+        bs = Bispectrum(self.lc)
+        lags = np.array([-2, -1, 0, 1, 2])
+        freq = np.array([-0.5, -0.25, 0., 0.25, 0.5])
+        cum3 = np.array([[0.0576, 0.1216, 0.1376, -0.2176, -0.0448],
+                         [0.1216, -0.6752, 0.4128, 0.1216, -0.0096],
+                         [0.1376, 0.4128, 0.288, -0.6752, 0.0576],
+                         [-0.2176, 0.1216, -0.6752, 0.4128, 0.1216],
+                         [-0.0448, -0.0096, 0.0576, 0.1216, 0.1376]])
+
+        bispec = np.array([[1.26572936 - 1.96410531e+00j, 1.01680863 - 2.33088331e+00j,
+                            0.53357150 - 3.20079088e-01j, -0.39280863 + 1.79793854e+00j,
+                            0.37973002 + 0.00000000e+00j],
+                           [1.01680863 - 2.33088331e+00j, -0.29772936 - 1.20447928e+00j,
+                            -0.11757150 - 7.55604229e-02j, 0.03626998 + 1.11022302e-16j,
+                            -0.39280863 - 1.79793854e+00j],
+                           [0.53357150 - 3.20079088e-01j, -0.11757150 - 7.55604229e-02j,
+                            0.27200000 + 0.00000000e+00j, -0.11757150 + 7.55604229e-02j,
+                            0.53357150 + 3.20079088e-01j],
+                           [-0.39280863 + 1.79793854e+00j, 0.03626998 - 1.11022302e-16j,
+                            -0.11757150 + 7.55604229e-02j, -0.29772936 + 1.20447928e+00j,
+                            1.01680863 + 2.33088331e+00j],
+                           [0.37973002 + 0.00000000e+00j, -0.39280863 - 1.79793854e+00j,
+                            0.53357150 + 3.20079088e-01j, 1.01680863 + 2.33088331e+00j,
+                            1.26572936 + 1.96410531e+00j]])
+
+        bispec_mag = np.array([[2.33661732, 2.54301333, 0.62221312, 1.84034823, 0.37973002],
+                               [2.54301333, 1.24073087, 0.13975849, 0.03626998, 1.84034823],
+                               [0.62221312, 0.13975849, 0.272, 0.13975849, 0.62221312],
+                               [1.84034823, 0.03626998, 0.13975849, 1.24073087, 2.54301333],
+                               [0.37973002, 1.84034823, 0.62221312, 2.54301333, 2.33661732]])
+
+        bispec_phase = np.array([[-9.98346367e-01, -1.15944968e+00, -5.40331561e-01, 1.78589369e+00,
+                                  0.00000000e+00],
+                                 [-1.15944968e+00, -1.81312395e+00, -2.57038310e+00, 3.06099713e-15,
+                                  -1.78589369e+00],
+                                 [-5.40331561e-01, -2.57038310e+00, 0.00000000e+00, 2.57038310e+00, 5.40331561e-01],
+                                 [1.78589369e+00, -3.06099713e-15, 2.57038310e+00, 1.81312395e+00,
+                                  1.15944968e+00],
+                                 [0.00000000e+00, -1.78589369e+00, 5.40331561e-01, 1.15944968e+00,
+                                  9.98346367e-01]])
+        assert bs.lc == self.lc
+        assert np.isclose(bs.fs, 1)
+        assert bs.maxlag == 2
+        assert bs.n == 5
+        assert bs.scale == 'biased'
+        assert np.allclose(bs.lags, lags)
+        assert np.allclose(bs.freq, freq)
+        assert np.allclose(bs.cum3, cum3)
+        assert np.allclose(bs.bispec, bispec)
+        assert np.allclose(bs.bispec_mag, bispec_mag)
+        assert np.allclose(bs.bispec_phase, bispec_phase)
+
+    def test_wrong_scale_type(self):
+        with pytest.raises(TypeError):
+            Bispectrum(self.lc, scale=1)
+
+    def test_wrong_scale_value(self):
+        with pytest.raises(ValueError):
+            Bispectrum(self.lc, scale='non-biased')
+
+    def test_bispectrum_unbiased_scale(self):
+        bs = Bispectrum(self.lc, scale='unbiased')
+
+        lags = np.array([-2, -1, 0, 1, 2])
+        freq = np.array([-0.5, -0.25, 0., 0.25, 0.5])
+
+        cum3 = np.array([[0.096, 0.20266667, 0.22933333, - 0.544, - 0.224],
+                         [0.20266667, - 0.844, 0.516, 0.20266667, - 0.024],
+                         [0.22933333, 0.516, 0.288, - 0.844, 0.096],
+                         [-0.544, 0.20266667, - 0.844, 0.516, 0.20266667],
+                         [-0.224, - 0.024, 0.096, 0.20266667, 0.22933333]])
+
+        bispec = np.array([[1.78211085 - 2.10567238e+00j, 1.56502881 - 3.02261738e+00j,
+                            0.63208770 - 8.00197720e-01j, - 1.20769548 + 2.56558544e+00j,
+                            0.49792363 + 2.22044605e-16j],
+                           [1.56502881 - 3.02261738e+00j, - 1.12477752 - 1.08193727e+00j,
+                            0.12524563 - 1.88901057e-01j, 0.25940971 + 2.22044605e-16j,
+                            - 1.20769548 - 2.56558544e+00j],
+                           [0.63208770 - 8.00197720e-01j,
+                            0.12524563 - 1.88901057e-01j,
+                            - 0.08800000 + 0.00000000e+00j,
+                            0.12524563 + 1.88901057e-01j,
+                            0.63208770 + 8.00197720e-01j],
+                           [-1.20769548 + 2.56558544e+00j, 0.25940971 - 2.22044605e-16j,
+                            0.12524563 + 1.88901057e-01j, - 1.12477752 + 1.08193727e+00j,
+                            1.56502881 + 3.02261738e+00j],
+                           [0.49792363 - 2.22044605e-16j, - 1.20769548 - 2.56558544e+00j,
+                            0.63208770 + 8.00197720e-01j,
+                            1.56502881 + 3.02261738e+00j,
+                            1.78211085 + 2.10567238e+00j]])
+
+        bispec_mag = np.array([[2.75858211, 3.40375249, 1.01973097, 2.83562286, 0.49792363],
+                               [3.40375249, 1.56067701, 0.22664968, 0.25940971, 2.83562286],
+                               [1.01973097, 0.22664968, 0.088, 0.22664968, 1.01973097],
+                               [2.83562286, 0.25940971, 0.22664968, 1.56067701, 3.40375249],
+                               [0.49792363, 2.83562286, 1.01973097, 3.40375249, 2.75858211]])
+
+        bispec_phase = np.array([[-8.68432005e-01, -1.09303185e+00, -9.02235465e-01, 2.01075415e+00,
+                                  4.45941091e-16],
+                                 [-1.09303185e+00, -2.37560564e+00, -9.85320934e-01, 8.55961046e-16,
+                                  -2.01075415e+00],
+                                 [-9.02235465e-01, -9.85320934e-01, 3.14159265e+00, 9.85320934e-01,
+                                  9.02235465e-01],
+                                 [2.01075415e+00, -8.55961046e-16, 9.85320934e-01, 2.37560564e+00,
+                                  1.09303185e+00],
+                                 [-4.45941091e-16, -2.01075415e+00, 9.02235465e-01, 1.09303185e+00,
+                                  8.68432005e-01]])
+
+        assert bs.lc == self.lc
+        assert np.isclose(bs.fs, 1)
+        assert bs.maxlag == 2
+        assert bs.n == 5
+        assert bs.scale == 'unbiased'
+        assert np.allclose(bs.lags, lags)
+        assert np.allclose(bs.freq, freq)
+        assert np.allclose(bs.cum3, cum3)
+        assert np.allclose(bs.bispec, bispec)
+        assert np.allclose(bs.bispec_mag, bispec_mag)
+        assert np.allclose(bs.bispec_phase, bispec_phase)
+
+    def test_lc1_with_diff_lag(self):
+        bs = Bispectrum(self.lc1, maxlag=1)
+        lags = np.array([-0.5, 0, 0.5])
+        freq = np.array([-1, 0., 1])
+        cum3 = np.array([[0.37114198, -0.62885802, -0.02160494],
+                         [-0.62885802, 0.59259259, 0.37114198],
+                         [-0.02160494, 0.37114198, -0.62885802]])
+
+        bispec = np.array([[0.93595679 + 2.59807621e+00j, 0.61419753 + 0j,
+                            0.61419753 + 2.22044605e-16j],
+                           [0.61419753 + 1.11022302e-16j, -0.22376543 + 0j,
+                            0.61419753 - 1.11022302e-16j],
+                           [0.61419753 - 2.22044605e-16j, 0.61419753 + 0j,
+                            0.93595679 - 2.59807621e+00j]])
+
+        bispec_mag = np.array([[2.76152406, 0.61419753, 0.61419753],
+                               [0.61419753, 0.22376543, 0.61419753],
+                               [0.61419753, 0.61419753, 2.76152406]])
+
+        bispec_phase = np.array([[1.22501950e+00, 0, 3.61519859e-16],
+                                 [1.80759930e-16, 3.14159265e+00, -1.80759930e-16],
+                                 [-3.61519859e-16, 0, -1.22501950e+00]])
+
+        assert bs.lc == self.lc1
+        assert np.isclose(bs.fs, 2)
+        assert bs.maxlag == 1
+        assert bs.n == 6
+        assert bs.scale == 'biased'
+        assert np.allclose(bs.lags, lags)
+        assert np.allclose(bs.freq, freq)
+        assert np.allclose(bs.cum3, cum3)
+        assert np.allclose(bs.bispec, bispec)
+        assert np.allclose(bs.bispec_mag, bispec_mag)
+        assert np.allclose(bs.bispec_phase, bispec_phase)
+
+    def test_lc1_unbiased_scale(self):
+        bs = Bispectrum(self.lc1, maxlag=1, scale='unbiased')
+        lags = np.array([-0.5, 0, 0.5])
+        freq = np.array([-1, 0., 1])
+
+        cum3 = np.array([[0.44537037, -0.75462963, -0.03240741],
+                         [-0.75462963, 0.59259259, 0.44537037],
+                         [-0.03240741, 0.44537037, -0.75462963]])
+
+        bispec = np.array([[0.99166667 + 3.11769145e+00j, 0.62500000 + 0j,
+                            0.62500000 + 0j],
+                           [0.62500000 + 2.22044605e-16j, -0.40000000 + 0j,
+                            0.62500000 - 2.22044605e-16j],
+                           [0.62500000 + 0j, 0.62500000 + 0j,
+                            0.99166667 - 3.11769145e+00j]])
+
+        bispec_mag = np.array([[3.27160554, 0.625, 0.625],
+                               [0.625, 0.4, 0.625],
+                               [0.625, 0.625, 3.27160554]])
+
+        bispec_phase = np.array([[1.26283852, 0, 0],
+                                 [0, 3.14159265, 0],
+                                 [0, 0, -1.26283852]])
+
+        assert bs.lc == self.lc1
+        assert np.isclose(bs.fs, 2)
+        assert bs.maxlag == 1
+        assert bs.n == 6
+        assert bs.scale == 'unbiased'
+        assert np.allclose(bs.lags, lags)
+        assert np.allclose(bs.freq, freq)
+        assert np.allclose(bs.cum3, cum3)
+        assert np.allclose(bs.bispec, bispec)
+        assert np.allclose(bs.bispec_mag, bispec_mag)
+        assert np.allclose(bs.bispec_phase, bispec_phase)
+
+    @pytest.mark.skipif(HAS_MPL, reason='Matplotlib is already installed if condition is met')
+    def test_plot_matplotlib_not_installed(self):
+        bs = Bispectrum(self.lc)
+        with pytest.raises(ImportError) as excinfo:
+            bs.plot_cum3()
+            bs.plot_mag()
+            bs.plot_phase()
+        message = str(excinfo.value)
+        assert "Matplotlib required for plot()" in message
+
+    @pytest.mark.skipif(not HAS_MPL, reason='Matplotlib is not installed')
+    def test_plot_cum3(self):
+        bs = Bispectrum(self.lc)
+        bs.plot_cum3()
+        assert plt.fignum_exists(1)
+
+    @pytest.mark.skipif(not HAS_MPL, reason='Matplotlib is not installed')
+    def test_plot_mag(self):
+        bs = Bispectrum(self.lc)
+        bs.plot_mag()
+        assert plt.fignum_exists(1)
+
+    @pytest.mark.skipif(not HAS_MPL, reason='Matplotlib is not installed')
+    def test_plot_phase(self):
+        bs = Bispectrum(self.lc)
+        bs.plot_phase()
+        assert plt.fignum_exists(1)
+
+    @pytest.mark.skipif(not HAS_MPL, reason='Matplotlib is not installed')
+    def test_plot_cum3_axis(self):
+        bs = Bispectrum(self.lc)
+        bs.plot_cum3(axis=[0, 1, 0, 100])
+        assert plt.fignum_exists(1)
+
+    @pytest.mark.skipif(not HAS_MPL, reason='Matplotlib is not installed')
+    def test_plot_mag_axis(self):
+        bs = Bispectrum(self.lc)
+        bs.plot_mag(axis=[0, 1, 0, 100])
+        assert plt.fignum_exists(1)
+
+    @pytest.mark.skipif(not HAS_MPL, reason='Matplotlib is not installed')
+    def test_plot_phase_axis(self):
+        bs = Bispectrum(self.lc)
+        bs.plot_phase(axis=[0, 1, 0, 100])
+        assert plt.fignum_exists(1)
+
+    @pytest.mark.skipif(not HAS_MPL, reason='Matplotlib is not installed')
+    def test_plot_cum3_default_filename(self):
+        bs = Bispectrum(self.lc)
+        bs.plot_cum3(save=True)
+        assert os.path.isfile('bispec_cum3.png')
+        os.unlink('bispec_cum3.png')
+
+    @pytest.mark.skipif(not HAS_MPL, reason='Matplotlib is not installed')
+    def test_plot_mag_default_filename(self):
+        bs = Bispectrum(self.lc)
+        bs.plot_mag(save=True)
+        assert os.path.isfile('bispec_mag.png')
+        os.unlink('bispec_mag.png')
+
+    @pytest.mark.skipif(not HAS_MPL, reason='Matplotlib is not installed')
+    def test_plot_phase_default_filename(self):
+        bs = Bispectrum(self.lc)
+        bs.plot_phase(save=True)
+        assert os.path.isfile('bispec_phase.png')
+        os.unlink('bispec_phase.png')
+
+    @pytest.mark.skipif(not HAS_MPL, reason='Matplotlib is not installed')
+    def test_plot_cum3_custom_filename(self):
+        bs = Bispectrum(self.lc)
+        bs.plot_cum3(save=True, filename='cum3.png')
+        assert os.path.isfile('cum3.png')
+        os.unlink('cum3.png')
+
+    @pytest.mark.skipif(not HAS_MPL, reason='Matplotlib is not installed')
+    def test_plot_mag_custom_filename(self):
+        bs = Bispectrum(self.lc)
+        bs.plot_phase(save=True, filename='mag.png')
+        assert os.path.isfile('mag.png')
+        os.unlink('mag.png')
+
+    @pytest.mark.skipif(not HAS_MPL, reason='Matplotlib is not installed')
+    def test_plot_phase_custom_filename(self):
+        bs = Bispectrum(self.lc)
+        bs.plot_phase(save=True, filename='phase.png')
+        assert os.path.isfile('phase.png')
+        os.unlink('phase.png')
\ No newline at end of file