Repository: StingraySoftware/stingray
Pull request #242 (2017-08-13 10:50:25)

Window Functions and 2D windows for Bispectrum


Patch:
diff --git a/stingray/bispectrum.py b/stingray/bispectrum.py
index 06e5280c..0a9fdd61 100644
--- a/stingray/bispectrum.py
+++ b/stingray/bispectrum.py
@@ -2,15 +2,15 @@ from __future__ import division
 
 import numpy as np
 from scipy.linalg import toeplitz
-from scipy.fftpack import fftshift, fft2, ifftshift, fftfreq
+from scipy.fftpack import fftshift, fft2, ifftshift, fft
+from  scipy.linalg import hankel
 
 from stingray import lightcurve
-from stingray.utils import simon
-import matplotlib.pyplot as plt
+import stingray.utils as utils
 
 
 class Bispectrum(object):
-    def __init__(self, lc, maxlag=None, scale='biased'):
+    def __init__(self, lc, maxlag=None, window=None, scale='biased'):
         """
                 Makes a :class:`Bispectrum` object from a given :class:`Lightcurve`.
                 
@@ -25,8 +25,11 @@ class Bispectrum(object):
                     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.
+                window : {'uniform', 'parzen', 'hamming', 'hanning', 'traingular', 'welch', 'blackmann', 'flat-top'}, optional, default 'uniform'
+                    Type of Window to apply for Bispectrum.
                 scale: {'biased', 'unbiased'}, optional, default 'biased'
                     Flag to decide biased or unbiased normalization for 3rd order cumulant function.
+                    
 
                 Attributes
                 ----------
@@ -93,9 +96,9 @@ class Bispectrum(object):
         """
 
         # Function call to create Bispectrum Object
-        self._make_bispectrum(lc, maxlag, scale)
+        self._make_bispetrum(lc, maxlag, window, scale)
 
-    def _make_bispectrum(self, lc, maxlag, scale):
+    def _make_bispetrum(self, lc, maxlag, window, scale):
         """
             Makes a Bispectrum Object with given lighcurve, maxlag and scale. 
         """
@@ -103,6 +106,16 @@ class Bispectrum(object):
         if not isinstance(lc, lightcurve.Lightcurve):
             raise TypeError('lc must be a lightcurve.ightcurve object')
 
+        # Available Windows. Used to resolve window paramneter
+        WINDOWS = ['uniform', 'parzen', 'hamming', 'hanning', 'triangular', 'welch', 'blackmann', 'flat-top']
+
+        if window:
+            if not isinstance(window, str):
+                raise TypeError('Window must be specified as string!')
+            window = window.lower()
+            if window not in WINDOWS:
+                raise ValueError("Wrong window specified or window function is not available")
+
         self.lc = lc
         self.fs = 1 / lc.dt
         self.n = self.lc.n
@@ -127,6 +140,13 @@ class Bispectrum(object):
             raise ValueError("scale can only be either 'biased' or 'unbiased'.")
         self.scale = scale.lower()
 
+        if window is None:
+            self.window_name = 'No Window'
+            self.window = None
+        else:
+            self.window_name = window
+            self.window = self._get_window()
+
         # Other Atributes
         self.lags = None
         self.cum3 = None
@@ -145,6 +165,28 @@ class Bispectrum(object):
         self._normalize_cumulant3()
         self._cal_bispec()
 
+    def _get_window(self):
+        """
+            Returns a window function of self.window_name type
+        """
+        N = 2 * self.maxlag + 1
+        window_even = utils.create_window(N, self.window_name)
+
+        # 2d even window
+        window2d = np.array([window_even, ] * N)
+
+        ## One-sided window with zero padding
+        window = np.zeros(N)
+        window[:self.maxlag + 1] = window_even[self.maxlag:]
+        window[self.maxlag:] = 0
+
+        # 2d window function to apply to bispectrum
+        row = np.concatenate(([window[0]], np.zeros(2 * self.maxlag)))
+        toep_matrix = toeplitz(window, row)
+        toep_matrix += np.tril(toep_matrix, -1).transpose()
+        window = toep_matrix[..., ::-1] * window2d * window2d.transpose()
+        return window
+
     def _cumulant3(self):
         """
             Calculates the 3rd Order cummulant of the lightcurve.
@@ -236,7 +278,13 @@ class Bispectrum(object):
             self.bispec_phase
         """
         self.freq = (1 / 2) * self.fs * (self.lags / self.lc.dt) / self.maxlag
-        self.bispec = fftshift(fft2(ifftshift(self.cum3)))
+
+        # Apply window if specified otherwise calculate with applying window
+        if self.window is None:
+            self.bispec = fftshift(fft2(ifftshift(self.cum3)))
+        else:
+            self.bispec = fftshift(fft2(ifftshift(self.cum3 * self.window)))
+
         self.bispec_mag = np.abs(self.bispec)
         self.bispec_phase = np.angle((self.bispec))
 
@@ -367,8 +415,3 @@ class Bispectrum(object):
             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)
diff --git a/stingray/utils.py b/stingray/utils.py
index a9695b33..8eaebc1d 100644
--- a/stingray/utils.py
+++ b/stingray/utils.py
@@ -3,6 +3,7 @@ from __future__ import (absolute_import, unicode_literals, division,
 import sys
 import collections
 import numbers
+from six import string_types
 
 import warnings
 import numpy as np
@@ -13,10 +14,12 @@ import numpy as np
 HAS_NUMBA = False
 try:
     from numba import jit
+
     HAS_NUMBA = True
 except ImportError:
     warnings.warn("Numba not installed. Faking it")
 
+
     class jit(object):
 
         def __init__(self, *args, **kwargs):
@@ -30,7 +33,7 @@ except ImportError:
 
 
 def _root_squared_mean(array):
-    return np.sqrt(np.sum(array**2)) / len(array)
+    return np.sqrt(np.sum(array ** 2)) / len(array)
 
 
 def simon(message, **kwargs):
@@ -52,7 +55,6 @@ def simon(message, **kwargs):
 
 
 def rebin_data(x, y, dx_new, yerr=None, method='sum', dx=None):
-
     """Rebin some data to an arbitrary new data resolution. Either sum
     the data points in the new bins or average them.
 
@@ -116,17 +118,17 @@ def rebin_data(x, y, dx_new, yerr=None, method='sum', dx=None):
         prev_frac = int_i + 1 - i
         prev_bin = int_i
         total += prev_frac * y[prev_bin]
-        totalerr += prev_frac * (yerr[prev_bin]**2)
+        totalerr += prev_frac * (yerr[prev_bin] ** 2)
 
         if i + step_size < len(x):
             # Fractional part of next bin:
             next_frac = i + step_size - int(i + step_size)
             next_bin = int(i + step_size)
             total += next_frac * y[next_bin]
-            totalerr += next_frac * (yerr[next_bin]**2)
+            totalerr += next_frac * (yerr[next_bin] ** 2)
 
-        total += sum(y[int(i+1):int(i+step_size)])
-        totalerr += sum(yerr[int(i+1):int(step_size)]**2)
+        total += sum(y[int(i + 1):int(i + step_size)])
+        totalerr += sum(yerr[int(i + 1):int(step_size)] ** 2)
         output.append(total)
         outputerr.append(np.sqrt(totalerr))
 
@@ -151,7 +153,7 @@ def rebin_data(x, y, dx_new, yerr=None, method='sum', dx=None):
         ybin = ybin[:-1]
         ybinerr = ybinerr[:-1]
 
-    new_x0 = (x[0] - (0.5*dx_old)) + (0.5*dx_new)
+    new_x0 = (x[0] - (0.5 * dx_old)) + (0.5 * dx_new)
     xbin = np.arange(ybin.shape[0]) * dx_new + new_x0
 
     return xbin, ybin, ybinerr, step_size
@@ -214,7 +216,7 @@ def rebin_data_log(x, y, f, y_err=None, dx=None):
     # until we reach the maximum frequency, increase the width of each
     # frequency bin by f
     while binx[-1] <= maxx:
-        binx.append(binx[-1] + dx*(1.0+f))
+        binx.append(binx[-1] + dx * (1.0 + f))
         dx = binx[-1] - binx[-2]
 
     # compute the mean of the ys that fall into each new frequency bin.
@@ -324,7 +326,7 @@ def is_int(obj):
     return isinstance(obj, (numbers.Integral, np.integer))
 
 
-def get_random_state(random_state = None):
+def get_random_state(random_state=None):
     if not random_state:
         random_state = np.random.mtrand._rand
     else:
@@ -332,7 +334,7 @@ def get_random_state(random_state = None):
             random_state = np.random.RandomState(random_state)
         elif not isinstance(random_state, np.random.RandomState):
             raise ValueError("{value} can't be used to generate a numpy.random.RandomState".format(
-                value = random_state
+                value=random_state
             ))
 
     return random_state
@@ -345,7 +347,7 @@ def baseline_als(y, lam, p, niter=10):
     https://www.researchgate.net/publication/228961729_Technical_Report_Baseline_Correction_with_Asymmetric_Least_Squares_Smoothing
     The Python translation is partly from
     http://stackoverflow.com/questions/29156532/python-baseline-correction-library
-    
+
     Parameters
     ----------
     y : array of floats
@@ -364,8 +366,8 @@ def baseline_als(y, lam, p, niter=10):
     for i in range(niter):
         W = sparse.spdiags(w, 0, L, L)
         Z = W + lam * D.dot(D.transpose())
-        z = sparse.linalg.spsolve(Z, w*y)
-        w = p * (y > z) + (1-p) * (y < z)
+        z = sparse.linalg.spsolve(Z, w * y)
+        w = p * (y > z) + (1 - p) * (y < z)
     return z
 
 
@@ -406,3 +408,95 @@ def excess_variance(lc, normalization='fvar'):
     elif normalization == 'none' or normalization is None:
         return var_xs, var_xs_err
 
+
+def create_window(N, window_type='uniform'):
+    """ A method to create window functions commonly used in signal processing.
+        Windows supported are:
+        Hamming, Hanning, uniform(rectangular window), triangular window, blackmann window among others.
+
+        Parameters
+        ----------
+        N : int
+            Total number of data points in window. If negative, abs is taken.
+        window_type : {'uniform', 'parzen', 'hamming', 'hanning', 'traingular', 'welch', 'blackmann', 'flat-top'}, optional, default 'uniform'
+            Type of window to create.
+         Returns
+        -------
+        window: numpy.ndarray
+            Window function of length N.
+    """
+
+    if not isinstance(N, int):
+        raise TypeError('N (window length) must be an integer')
+
+    WINDOWS = ['uniform', 'parzen', 'hamming', 'hanning', 'triangular', 'welch', 'blackmann', 'flat-top']
+
+    if not isinstance(window_type, string_types):
+        raise TypeError('type of window must be specified as string!')
+
+    window_type = window_type.lower()
+    if window_type not in WINDOWS:
+        raise ValueError("Wrong window type specified or window function is not available")
+
+    # Return empty array as window if N = 0
+    if N == 0:
+        return np.array([])
+
+    window = None
+    N = abs(N)
+
+    # Window samples index
+    n = np.arange(N)
+
+    # Constants
+    N_minus_1 = N - 1
+    N_by_2 = np.int((np.floor((N_minus_1) / 2)))
+
+    # Create Windows
+    if window_type == 'uniform':
+        window = np.ones(N)
+
+    if window_type == 'parzen':
+        N_parzen = np.int(np.ceil((N + 1) / 2))
+        N2_plus_1 = np.int(np.floor((N_parzen / 2))) + 1
+
+        window = np.zeros(N_parzen)
+        windlag0 = np.arange(0, N2_plus_1) / (N_parzen - 1)
+        windlag1 = 1 - np.arange(N2_plus_1, N_parzen) / (N_parzen - 1)
+        window[:N2_plus_1] = 1 - (1 - windlag0) * windlag0 * windlag0 * 6
+        window[N2_plus_1:] = windlag1 * windlag1 * windlag1 * 2
+        lagindex = np.arange(N_parzen - 1, 0, -1)
+        window = np.concatenate((window[lagindex], window))
+        window = window[:N]
+
+    if window_type == 'hamming':
+        window = 0.54 - 0.46 * np.cos((2 * np.pi * n) / N_minus_1)
+
+    if window_type == 'hanning':
+        window = 0.5 * (1 - np.cos(2 * np.pi * n / N_minus_1))
+
+    if window_type == 'triangular':
+        window = 1 - np.abs((n - (N_by_2)) / N)
+
+    if window_type == 'welch':
+        N_minus_1_by_2 = N_minus_1 / 2
+        window = 1 - np.square((n - N_minus_1_by_2) / N_minus_1_by_2)
+        
+    if window_type == 'blackmann':
+        a0 = 0.42659
+        a1 = 0.49656
+        a2 = 0.076849
+        window = a0 - a1 * np.cos((2 * np.pi * n) / N_minus_1) + a2 * np.cos((4 * np.pi * n) / N_minus_1)
+
+    if window_type == 'flat-top':
+        a0 = 1
+        a1 = 1.93
+        a2 = 1.29
+        a3 = 0.388
+        a4 = 0.028
+        window = a0 - a1 * np.cos((2 * np.pi * n) / N_minus_1) + \
+                 a2 * np.cos((4 * np.pi * n) / N_minus_1) - \
+                 a3 * np.cos((6 * np.pi * n) / N_minus_1) + \
+                 a4 * np.cos((8 * np.pi * n) / N_minus_1)
+
+    return window
\ No newline at end of file


Test patch:
diff --git a/stingray/tests/test_bispectrum.py b/stingray/tests/test_bispectrum.py
index 814cef7f..789e0c17 100644
--- a/stingray/tests/test_bispectrum.py
+++ b/stingray/tests/test_bispectrum.py
@@ -236,6 +236,289 @@ class TestBispectrum(object):
         assert np.allclose(bs.bispec_mag, bispec_mag)
         assert np.allclose(bs.bispec_phase, bispec_phase)
 
+    def test_bispectrum_window_none(self):
+        bs = Bispectrum(self.lc, scale='unbiased')
+        assert bs.window is None
+        assert bs.window_name == 'No Window'
+
+    def test_bispectrum_window_uniform(self):
+        bs = Bispectrum(self.lc, maxlag=2, window='uniform')
+        lags = np.array([-2, -1, 0, 1, 2])
+        freq = np.array([-0.5, -0.25, 0., 0.25, 0.5])
+        bispec = np.array([[1.23378944 - 1.03450204e+00j, 0.52344952 - 1.79651918e+00j,
+                            0.12704520 - 9.59589442e-01j, -0.10264952 + 5.93058891e-01j,
+                            0.71979565 + 1.11022302e-16j],
+                           [0.52344952 - 1.79651918e+00j, 0.13901056 - 1.67385947e+00j,
+                            0.20575480 - 1.11030991e+00j, 0.26900435 + 0.00000000e+00j,
+                            -0.10264952 - 5.93058891e-01j],
+                           [0.12704520 - 9.59589442e-01j, 0.20575480 - 1.11030991e+00j,
+                            -0.53760000 + 0.00000000e+00j, 0.20575480 + 1.11030991e+00j,
+                            0.12704520 + 9.59589442e-01j],
+                           [-0.10264952 + 5.93058891e-01j, 0.26900435 + 0.00000000e+00j,
+                            0.20575480 + 1.11030991e+00j, 0.13901056 + 1.67385947e+00j,
+                            0.52344952 + 1.79651918e+00j],
+                           [0.71979565 - 1.11022302e-16j, -0.10264952 - 5.93058891e-01j,
+                            0.12704520 + 9.59589442e-01j, 0.52344952 + 1.79651918e+00j,
+                            1.23378944 + 1.03450204e+00j]])
+
+        window = np.array([[0, 0, 0, 1, 1], 
+                  [0, 0, 1, 1, 1],
+                  [0, 1, 1, 1, 0],
+                  [1, 1, 1, 0, 0],
+                  [1, 1, 0, 0, 0]])
+
+        assert bs.lc == self.lc
+        assert np.isclose(bs.fs, 1)
+        assert bs.maxlag == 2
+        assert bs.n == 5
+        assert np.allclose(bs.lags, lags)
+        assert np.allclose(bs.freq, freq)
+        assert np.allclose(bs.bispec, bispec)
+        assert bs.window_name == 'uniform'
+        assert np.allclose(bs.window, window)
+
+    def test_bispectrum_window_parzen(self):
+        bs = Bispectrum(self.lc, maxlag=2, window='parzen')
+        bispec = np.array([[0.32973576 - 7.99387943e-02j, 0.30089706 - 1.04641240e-01j,
+                            0.27257082 - 3.99693972e-02j, 0.28390294 + 2.47024460e-02j,
+                            0.31923282 + 6.93889390e-18j],
+                           [0.30089706 - 1.04641240e-01j, 0.29306424 - 1.29343686e-01j,
+                            0.27122918 - 6.46718431e-02j, 0.26556718 + 0.00000000e+00j,
+                            0.28390294 - 2.47024460e-02j],
+                           [0.27257082 - 3.99693972e-02j, 0.27122918 - 6.46718431e-02j,
+                            0.27040000 + 0.00000000e+00j, 0.27122918 + 6.46718431e-02j,
+                            0.27257082 + 3.99693972e-02j],
+                           [0.28390294 + 2.47024460e-02j, 0.26556718 + 0.00000000e+00j,
+                            0.27122918 + 6.46718431e-02j, 0.29306424 + 1.29343686e-01j,
+                            0.30089706 + 1.04641240e-01j],
+                           [0.31923282 - 6.93889390e-18j, 0.28390294 - 2.47024460e-02j,
+                            0.27257082 + 3.99693972e-02j, 0.30089706 + 1.04641240e-01j,
+                            0.32973576 + 7.99387943e-02j]])
+
+        window = np.array([[0, 0, 0, 0, 0, ],
+                  [0., 0., 0.0625, 0.0625, 0.],
+                  [0., 0.0625, 1., 0.0625, 0.],
+                  [0., 0.0625, 0.0625, 0., 0.],
+                  [0., 0., 0., 0., 0.]])
+
+        assert bs.lc == self.lc
+        assert np.isclose(bs.fs, 1)
+        assert bs.maxlag == 2
+        assert bs.n == 5
+        assert np.allclose(bs.bispec, bispec)
+        assert bs.window_name == 'parzen'
+        assert np.allclose(bs.window, window)
+
+    def test_bispectrum_window_hamming(self):
+        bs = Bispectrum(self.lc, maxlag=2, window='hamming')
+        bispec = np.array([[0.49072469 - 3.67258307e-01j, 0.34762420 - 4.91066236e-01j,
+                            0.21848648 - 1.93948024e-01j, 0.26176961 + 1.19866471e-01j,
+                            0.43090809 + 5.55111512e-17j],
+                           [0.34762420 - 4.91066236e-01j, 0.30777863 - 5.94236424e-01j,
+                            0.21286804 - 3.03495625e-01j, 0.19173603 + 5.55111512e-17j,
+                            0.26176961 - 1.19866471e-01j],
+                           [0.21848648 - 1.93948024e-01j, 0.21286804 - 3.03495625e-01j,
+                            0.19471176 + 0.00000000e+00j, 0.21286804 + 3.03495625e-01j,
+                            0.21848648 + 1.93948024e-01j],
+                           [0.26176961 + 1.19866471e-01j, 0.19173603 - 5.55111512e-17j,
+                            0.21286804 + 3.03495625e-01j, 0.30777863 + 5.94236424e-01j,
+                            0.34762420 + 4.91066236e-01j],
+                           [0.43090809 - 5.55111512e-17j, 0.26176961 - 1.19866471e-01j,
+                            0.21848648 + 1.93948024e-01j, 0.34762420 + 4.91066236e-01j,
+                            0.49072469 + 3.67258307e-01j]])
+
+        window = np.array([[0, 0, 0, 0.023328, 0.0064],
+                           [0., 0., 0.2916, 0.2916, 0.023328],
+                           [0., 0.2916, 1., 0.2916, 0.],
+                           [0.023328, 0.2916, 0.2916, 0., 0.],
+                           [0.0064, 0.023328, 0., 0., 0.]])
+
+        assert bs.lc == self.lc
+        assert np.isclose(bs.fs, 1)
+        assert bs.maxlag == 2
+        assert bs.n == 5
+        assert np.allclose(bs.bispec, bispec)
+        assert bs.window_name == 'hamming'
+        assert np.allclose(bs.window, window)
+
+    def test_bispectrum_window_hanning(self):
+        bs = Bispectrum(self.lc, maxlag=2, window='hanning')
+        bispec = np.array([[0.45494303 - 3.19755177e-01j, 0.33958823 - 4.18564961e-01j,
+                            0.22628328 - 1.59877589e-01j, 0.27161177 + 9.88097838e-02j,
+                            0.41293126 + 5.55111512e-17j],
+                           [0.33958823 - 4.18564961e-01j, 0.30825697 - 5.17374745e-01j,
+                            0.22091672 - 2.58687372e-01j, 0.19826874 + 5.55111512e-17j,
+                            0.27161177 - 9.88097838e-02j],
+                           [0.22628328 - 1.59877589e-01j, 0.22091672 - 2.58687372e-01j,
+                            0.21760000 + 0.00000000e+00j, 0.22091672 + 2.58687372e-01j,
+                            0.22628328 + 1.59877589e-01j],
+                           [0.27161177 + 9.88097838e-02j, 0.19826874 - 5.55111512e-17j,
+                            0.22091672 + 2.58687372e-01j, 0.30825697 + 5.17374745e-01j,
+                            0.33958823 + 4.18564961e-01j],
+                           [0.41293126 - 5.55111512e-17j, 0.27161177 - 9.88097838e-02j,
+                            0.22628328 + 1.59877589e-01j, 0.33958823 + 4.18564961e-01j,
+                            0.45494303 + 3.19755177e-01j]])
+
+        window = np.array([[0., 0., 0., 0., 0.],
+                           [0., 0., 0.25, 0.25, 0.],
+                           [0., 0.25, 1., 0.25, 0.],
+                           [0., 0.25, 0.25, 0., 0.],
+                           [0., 0., 0., 0., 0.]]
+                          )
+
+        assert bs.lc == self.lc
+        assert np.isclose(bs.fs, 1)
+        assert bs.maxlag == 2
+        assert bs.n == 5
+        assert np.allclose(bs.bispec, bispec)
+        assert bs.window_name == 'hanning'
+        assert np.allclose(bs.window, window)
+
+    def test_bispectrum_window_triangular(self):
+        bs = Bispectrum(self.lc, maxlag=2, window='triangular')
+        bispec = np.array([[0.82428321 - 7.24678086e-01j, 0.42949926 - 1.11847388e+00j,
+                            0.16365995 - 5.32196997e-01j, 0.12187354 + 3.28915833e-01j,
+                            0.57999943 + 5.55111512e-17j],
+                           [0.42949926 - 1.11847388e+00j, 0.25368159 - 1.17255377e+00j,
+                            0.18598485 - 6.91254876e-01j, 0.18948537 + 1.11022302e-16j,
+                            0.12187354 - 3.28915833e-01j],
+                           [0.16365995 - 5.32196997e-01j, 0.18598485 - 6.91254876e-01j,
+                            -0.09896960 + 0.00000000e+00j, 0.18598485 + 6.91254876e-01j,
+                            0.16365995 + 5.32196997e-01j],
+                           [0.12187354 + 3.28915833e-01j, 0.18948537 - 1.11022302e-16j,
+                            0.18598485 + 6.91254876e-01j, 0.25368159 + 1.17255377e+00j,
+                            0.42949926 + 1.11847388e+00j],
+                           [0.57999943 - 5.55111512e-17j, 0.12187354 - 3.28915833e-01j,
+                            0.16365995 + 5.32196997e-01j, 0.42949926 + 1.11847388e+00j,
+                            0.82428321 + 7.24678086e-01j]])
+
+        window = np.array([[0., 0., 0., 0.384, 0.36],
+                           [0., 0., 0.64, 0.64, 0.384],
+                           [0., 0.64, 1., 0.64, 0.],
+                           [0.384, 0.64, 0.64, 0., 0.],
+                           [0.36, 0.384, 0., 0., 0.]])
+
+        assert bs.lc == self.lc
+        assert np.isclose(bs.fs, 1)
+        assert bs.maxlag == 2
+        assert bs.n == 5
+        assert np.allclose(bs.bispec, bispec)
+        assert bs.window_name == 'triangular'
+        assert np.allclose(bs.window, window)
+
+    def test_bispectrum_window_welch(self):
+        bs = Bispectrum(self.lc, maxlag=2, window='welch')
+        bispec = np.array([[0.66362182 - 7.19449149e-01j, 0.40407352 - 9.41771162e-01j,
+                            0.14913738 - 3.59724574e-01j, 0.25112648 + 2.22322014e-01j,
+                            0.56909534 + 5.55111512e-17j],
+                           [0.40407352 - 9.41771162e-01j, 0.33357818 - 1.16409318e+00j,
+                            0.13706262 - 5.82046588e-01j, 0.08610466 + 0.00000000e+00j,
+                            0.25112648 - 2.22322014e-01j],
+                           [0.14913738 - 3.59724574e-01j, 0.13706262 - 5.82046588e-01j,
+                            0.12960000 + 0.00000000e+00j, 0.13706262 + 5.82046588e-01j,
+                            0.14913738 + 3.59724574e-01j],
+                           [0.25112648 + 2.22322014e-01j, 0.08610466 + 0.00000000e+00j,
+                            0.13706262 + 5.82046588e-01j, 0.33357818 + 1.16409318e+00j,
+                            0.40407352 + 9.41771162e-01j],
+                           [0.56909534 - 5.55111512e-17j, 0.25112648 - 2.22322014e-01j,
+                            0.14913738 + 3.59724574e-01j, 0.40407352 + 9.41771162e-01j,
+                            0.66362182 + 7.19449149e-01j]])
+
+        window = np.array([[0., 0., 0., 0., 0.],
+                           [0., 0., 0.5625, 0.5625, 0.],
+                           [0., 0.5625, 1., 0.5625, 0.],
+                           [0., 0.5625, 0.5625, 0., 0.],
+                           [0., 0., 0., 0., 0.]])
+
+        assert bs.lc == self.lc
+        assert np.isclose(bs.fs, 1)
+        assert bs.maxlag == 2
+        assert bs.n == 5
+        assert np.allclose(bs.bispec, bispec)
+        assert bs.window_name == 'welch'
+        assert np.allclose(bs.window, window)
+
+    def test_bispectrum_window_blackmann(self):
+        bs = Bispectrum(self.lc, maxlag=2, window='blackmann')
+        bispec = np.array([[0.36998520 - 1.56242334e-01j, 0.31320687 - 2.04896068e-01j,
+                            0.25789699 - 7.84933643e-02j, 0.27972923 + 4.85115670e-02j,
+                            0.34901011 + 2.77555756e-17j],
+                           [0.31320687 - 2.04896068e-01j, 0.29778797 - 2.52805407e-01j,
+                            0.25527601 - 1.26632734e-01j, 0.24440392 + 2.77555756e-17j,
+                            0.27972923 - 4.85115670e-02j],
+                           [0.25789699 - 7.84933643e-02j, 0.25527601 - 1.26632734e-01j,
+                            0.25316762 + 0.00000000e+00j, 0.25527601 + 1.26632734e-01j,
+                            0.25789699 + 7.84933643e-02j],
+                           [0.27972923 + 4.85115670e-02j, 0.24440392 - 2.77555756e-17j,
+                            0.25527601 + 1.26632734e-01j, 0.29778797 + 2.52805407e-01j,
+                            0.31320687 + 2.04896068e-01j],
+                           [0.34901011 - 2.77555756e-17j, 0.27972923 - 4.85115670e-02j,
+                            0.25789699 + 7.84933643e-02j, 0.31320687 + 2.04896068e-01j,
+                            0.36998520 + 1.56242334e-01j]])
+
+        window = np.array([[0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 8.41430799e-04,
+                            4.73205937e-05],
+                           [0.00000000e+00, 0.00000000e+00, 1.22318645e-01, 1.22318645e-01,
+                            8.41430799e-04],
+                           [0.00000000e+00, 1.22318645e-01, 9.99997000e-01, 1.22318645e-01,
+                            0.00000000e+00],
+                           [8.41430799e-04, 1.22318645e-01, 1.22318645e-01, 0.00000000e+00,
+                            0.00000000e+00],
+                           [4.73205937e-05, 8.41430799e-04, 0.00000000e+00, 0.00000000e+00,
+                            0.00000000e+00]])
+
+        assert bs.lc == self.lc
+        assert np.isclose(bs.fs, 1)
+        assert bs.maxlag == 2
+        assert bs.n == 5
+        assert np.allclose(bs.bispec, bispec)
+        assert bs.window_name == 'blackmann'
+        assert np.allclose(bs.window, window)
+
+    def test_bispectrum_window_flattop(self):
+        bs = Bispectrum(self.lc, maxlag=2, window='flat-top')
+        bispec = np.array([[28.90860041 - 0.40702734j, 28.76176133 - 0.53280571j,
+                            28.61753157 - 0.20351367j, 28.67523175 + 0.12577837j, 28.85512219 + 0.j, ],
+                           [28.76176133 - 0.53280571j, 28.72187869 - 0.65858408j,
+                            28.61070029 - 0.32929204j, 28.58187089 + 0.j, 28.67523175 - 0.12577837j],
+                           [28.61753157 - 0.20351367j, 28.61070029 - 0.32929204j, 28.60647832 + 0.j,
+                            28.61070029 + 0.32929204j, 28.61753157 + 0.20351367j],
+                           [28.67523175 + 0.12577837j, 28.58187089 + 0.j, 28.61070029 + 0.32929204j,
+                            28.72187869 + 0.65858408j, 28.76176133 + 0.53280571j],
+                           [28.85512219 + 0.j, 28.67523175 - 0.12577837j,
+                            28.61753157 + 0.20351367j, 28.76176133 + 0.53280571j,
+                            28.90860041 + 0.40702734j]])
+
+        window = np.array([[0.00000000e+00, -0.00000000e+00, 0.00000000e+00, 5.95391791e-18,
+                            3.48773876e-32],
+                           [-0.00000000e+00, 0.00000000e+00, 3.18233584e-01, 3.18233584e-01,
+                            5.95391791e-18],
+                           [0.00000000e+00, 3.18233584e-01, 9.96392115e+01, 3.18233584e-01,
+                            0.00000000e+00],
+                           [5.95391791e-18, 3.18233584e-01, 3.18233584e-01, 0.00000000e+00,
+                            -0.00000000e+00],
+                           [3.48773876e-32, 5.95391791e-18, 0.00000000e+00, -0.00000000e+00,
+                            0.00000000e+00]]
+                          )
+
+        assert bs.lc == self.lc
+        assert np.isclose(bs.fs, 1)
+        assert bs.maxlag == 2
+        assert bs.n == 5
+        assert np.allclose(bs.bispec, bispec)
+        assert bs.window_name == 'flat-top'
+        assert np.allclose(bs.window, window)
+
+    def test_bad_window(self):
+        window_bad = 123
+        with pytest.raises(TypeError):
+            bs = Bispectrum(self.lc, maxlag=2, window=window_bad)
+
+    def test_not_available_window(self):
+        window_not = 'kaiser'
+        with pytest.raises(ValueError):
+            bs = Bispectrum(self.lc, maxlag=2, window=window_not)
+    
     @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)
diff --git a/stingray/tests/test_utils.py b/stingray/tests/test_utils.py
index 25b5ed93..2c666973 100644
--- a/stingray/tests/test_utils.py
+++ b/stingray/tests/test_utils.py
@@ -114,3 +114,77 @@ class TestUtils(object):
 
         with pytest.raises(ValueError):
             utils.get_random_state('foobar')
+
+class TestCreateWindow(object):
+    @classmethod
+    def setup_class(cls):
+        cls.N = 5
+        cls.uniform_window = 'uniform'
+        cls.parzen_window = 'parzen'
+        cls.hamming_window = 'hamming'
+        cls.hanning_window = 'hanning'
+        cls.triangular_window = 'triangular'
+        cls.welch_window = 'welch'
+        cls.blackmann_window = 'blackmann'
+        cls.flattop_window = 'flat-top'
+
+    def test_bad_N(self):
+        N_bad = 'abc'
+        with pytest.raises(TypeError):
+            window = utils.create_window(N_bad, self.uniform_window)
+
+    def test_bad_window_type(self):
+        window_bad = 123
+        with pytest.raises(TypeError):
+            window = utils.create_window(self.N, window_bad)
+
+    def test_not_available_window(self):
+        window_not = 'kaiser'
+        with pytest.raises(ValueError):
+            window = utils.create_window(self.N, window_not)
+
+    def test_N_equals_zero(self):
+        N = 0
+        window = utils.create_window(N)
+        assert len(window) == 0
+
+    def test_uniform_window(self):
+        result = np.ones(self.N)
+        window = utils.create_window(self.N)
+        assert np.allclose(window, result)
+
+    def test_parzen_window(self):
+        result = np.array([0, 0.25, 1, 0.25, 0])
+        window = utils.create_window(self.N, self.parzen_window)
+        assert np.allclose(window, result)
+
+    def test_hamming_window(self):
+        result = np.array([0.08, 0.54, 1, 0.54, 0.08])
+        window = utils.create_window(self.N, self.hamming_window)
+        assert np.allclose(window, result)
+
+    def test_hanning_window(self):
+        result = np.array([0, 0.5, 1, 0.5, 0])
+        window = utils.create_window(self.N, self.hanning_window)
+        assert np.allclose(window, result)
+
+    def test_triangular_window(self):
+        result = np.array([0.6, 0.8, 1, 0.8, 0.6])
+        window = utils.create_window(self.N, self.triangular_window)
+        assert np.allclose(window, result)
+
+    def test_welch_window(self):
+        result = np.array([0, 0.75, 1, 0.75, 0])
+        window = utils.create_window(self.N, self.welch_window)
+        assert np.allclose(window, result)
+
+    def test_blackmann_window(self):
+        result = np.array([0.006879, 0.349741, 0.999999, 0.349741, 0.006879])
+        window = utils.create_window(self.N, self.blackmann_window)
+        assert np.allclose(window, result)
+
+    def test_flat_top_window(self):
+        result = np.array([8.67361738e-17, -2.62000000e-01, 4.63600000e+00, -2.62000000e-01,
+                           8.67361738e-17])
+        window = utils.create_window(self.N, self.flattop_window)
+        assert np.allclose(window, result)