Borß–Martin windows in Python

borss_martin.py

# -*- coding: utf-8 -*-
"""
Created on Tue Aug 16 2016

Author of paper is anglicized as "borss" in URLs and email.
"""
from __future__ import division, print_function, absolute_import
from scipy.special import comb
from numpy import zeros_like, linspace, pi, sin, concatenate, array
import numpy as np


def _len_guards(M):
    """Handle small or incorrect window lengths"""
    if int(M) != M or M < 0:
        raise ValueError('Window length M must be a non-negative integer')
    return M <= 1


# Odd-length windows can be periodic, too.
# Same as https://github.com/scipy/scipy/pull/6483
def _extend(M, sym):
    """Extend window by 1 sample if needed for DFT-even symmetry"""
    if not sym:
        return M + 1, True
    else:
        return M, False


def _truncate(w, needed):
    """Truncate window by 1 sample if needed for DFT-even symmetry"""
    if needed:
        return w[:-1]
    else:
        return w


def borss_martin(M, noverlap, order=1, sym=True):
    u"""
    Return a window with constant overlap-add for arbitrary overlap.

    Generate a window function of length `M` that has the constant
    overlap-add ("COLA") property for `noverlap` samples of overlap.  The
    window is optimized for maximally-steep sidelobe fall-off for a given
    `order`.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero, an empty
        array is returned.
    noverlap : int
        Number of samples that each window overlaps the next.
    order : int
        Number of cosine terms in the expansion of the __TODO___
        Order of the ___TODO___.  Defaults to 1.  See Notes for details.
    sym : bool, optional
        When True, generates a symmetric window, for use in filter design.
        When False, generates a periodic window, for use in spectral analysis.

    Returns
    -------
    w : ndarray
        The window, with the maximum value normalized to 1 (though the value 1
        might not appear if the number of samples is even and sym is True).

    See Also
    --------
    bartlett, tukey, hann
    stft
    check_COLA

    Notes
    -----
    Borß and Martin point out that a window with COLA condition can be created
    for any arbitrary overlap by convolving 2 windows: The first a
    rectangular window of length equal to the hop size, and the second a
    window with an area of 1 and length equal to the overlap size.

    They then develop a family of such windows which optimizes the steepness
    of sidelobe fall-off for a given order `N`:

    - With ``order=0``, a trapezoidal window is produced, becoming a Bartlett
      window at 50% overlap.  This has side-lobe fall-off of 12 dB/octave.
    - With ``order=1``, a more smooth window is produced, equivalent to a
      Tukey window for overlap < 50%, or a Hann window at 50% overlap.
      This has side-lobe fall-off of 18 dB/octave.
    - Higher orders exhibit wider main lobes and faster fall-off, with slope of
      ``(N+2)*6`` dB/octave.

    Another way of stating this is that the boundaries of the window
    have `N` continuous derivatives.

    References
    ----------
    .. [1] C. Borß and R. Martin, "On the construction of window functions
           with constant-overlap-add constraint for arbitrary window shifts",
           2012 IEEE International Conference on Acoustics, Speech and Signal
           Processing (ICASSP), Kyoto, 2012, pp. 337-340.
           http://dx.doi.org/10.1109/ICASSP.2012.6287885

    Examples
    --------
    >>> scipy.signal.borss_martin(12, 14)  make a real example
    array([ 0.00000773,  0.00346009,  0.04652002,  0.22973712,  0.59988532,
            0.9456749 ,  0.9456749 ,  0.59988532,  0.22973712,  0.04652002,
            0.00346009,  0.00000773])

    Plot the window and the frequency response:

    >>> from numpy.fft import fft, fftshift
    >>> pass
    >>> # show a line for the fall-off rate

    show the side-lobe fall-off for different orders

    """
    if _len_guards(M):
        return np.ones(M)
    M, is_periodic = _extend(M, sym)

    L = noverlap
    if is_periodic:
        L = L + 1

    if L >= M:
        raise ValueError('Overlap length must be less than total length')
    if L < 0 or int(L) != L:
        raise ValueError('Overlap must be a non-negative integer')

    N = order

    if N < 0 or int(N) != N:
        raise ValueError('Order must be a non-negative integer')

    # includes both 0 and 1, so overlap works [0, 1/2, 1] + [1, 1/2, 0]
    t = linspace(-pi/2, pi/2, L)

    # Tack on a time value that will be the peak of the final function if
    # ramps are overlapping
    if L > M//2:
        t_val = (M - 1) / (L - 1) * pi/2 - pi/2
        t = concatenate((t, [t_val]))

    r = zeros_like(t)

    """
    From the paper, calculated the time-domain cos terms. Rect is equivalent
    to unit steps up and down, and convolution with step is equivalent to
    integration, so each end is the integral of the cos terms.  Numerical
    convolution/summation won't work, so calculated the antiderivatives of
    cos terms in the interval that made them simplest (symmetrical about the
    origin).

    1. Generate delta functions in frequency domain from binomial coefficients
    2. Transform to time domain (creates g(t) hump function)
    3. Normalize to area of 1

    Convolution with boxcar = (convolution with step) - (convolution with
    delayed step), and convolution with step = integration, so:

    4. Integrate over time, producing a curved ramp.  Formulas are simpler if
    the ramp is shifted so it's centered on the origin:
    """

    for k in range(N//2 + 1):
        # "where the individual weights are the binomial coefficients"
        wgt = comb(N, N//2 - k)

        # Odd N
        if N % 2:
            r += wgt * sin((2*k + 1)*t) / (2*k + 1)
        # Even N
        elif k == 0:
            r += wgt * t
        else:
            r += wgt * sin(2*k*t) / k

    if N % 2:  # Odd N
        r /= 2**N / comb(N/2., N//2)
    else:  # Even N
        r /= comb(N, N//2) * pi

    if L > M//2:
        # Split out normalization constant
        norm = r[-1]
        r = r[:-1]
    else:
        norm = 1/2.

    # Pad with flat line to complete ramped step
    head = concatenate((r, np.full(M - L, 0.5)))

    # Combine step up and step down
    win = head + head[::-1]

    # Normalize amplitude to 1
    win /= 2 * norm

    return _truncate(win, is_periodic)


##########################
# TESTS
##########################

from numpy.testing import assert_allclose, assert_raises
import windows as signal
from spectral import check_COLA
import pytest


def compare_win(w1, w2, L):
    assert check_COLA(w1, len(w1), L)
    assert check_COLA(w2, len(w2), L)
    assert_allclose(w1, w2)


def test_0_order_bartlett():
    """
    0th-order:
    A special case is produced for TS = TL = TT /2 (50% overlap) where the
    tapered window becomes the well known Bartlett (triangular) window.
    """
    bm = borss_martin(6, 3, 0, False)
    ba = signal.bartlett(6, False)
    w = [0, 1/3, 2/3, 3/3, 2/3, 1/3]
    compare_win(w, bm, 3)
    compare_win(ba, bm, 3)

    bm = borss_martin(7, 4, 0, True)
    ba = signal.bartlett(7, True)
    w = array([0, 1, 2, 3, 2, 1, 0])/3
    compare_win(w, bm, 4)
    compare_win(ba, bm, 4)

    # only difference is normalization:
    # borss_martin is a trapezoid with top = 1
    # bartlett is a triangle with an intersample peak = 1
    bm = borss_martin(6, 3, 0, True)
    ba = signal.bartlett(6, True)
    w = array([0, 2, 4, 4, 2, 0])/4
    compare_win(w, bm, 3)
    compare_win(ba * 5/4, bm, 3)

    bm = borss_martin(7, 4, 0, False)
    ba = signal.bartlett(7, False)
    w = array([0, 1, 2, 3, 3, 2, 1])/3
    compare_win(w, bm, 4)
    compare_win(ba * 3.5/3, bm, 4)


def test_0_overlap():
    # All ones is the only thing that makes sense for 0 overlap, but breaks
    # convention that otherwise starts on 0.
    # TODO: Special case or raise Exception?

    # TODO: Why do these work correctly?
    w = [1, 1, 1, 1, 1, 1]
    bm = borss_martin(6, 0, 0, True)
    compare_win(w, bm, 0)

    w = [1, 1, 1, 1, 1, 1, 1]
    bm = borss_martin(7, 0, 0, True)
    compare_win(w, bm, 0)

    w = [1, 1, 1, 1, 1, 1, 1]
    bm = borss_martin(7, 0, 1, True)
    compare_win(w, bm, 0)

#    w = [1, 1, 1, 1, 1, 1]
#    bm = borss_martin(6, 0, 0, False)
#    compare_win(w, bm, 0)
#
#    w = [1, 1, 1, 1, 1, 1, 1]
#    bm = borss_martin(7, 0, 0, False)
#    compare_win(w, bm, 0)


def test_1_overlap():
    w = [0, 1, 1, 1, 1, 1]
    bm = borss_martin(6, 1, 0, False)
    compare_win(w, bm, 1)

    w = [0, 1, 1, 1, 1, 1, 1]
    bm = borss_martin(7, 1, 0, False)
    compare_win(w, bm, 1)

    # Only way for sym to work is to be 1/2 on the ends, which breaks
    # convention that otherwise starts on 0.
    # TODO: Special case or raise Exception?

#    bm = borss_martin(6, 1, 0, True)
#    w = [1/2, 1, 1, 1, 1, 1/2]
#    compare_win(w, bm, 1)

#    bm = borss_martin(7, 1, 0, True)
#    w = [1/2, 1, 1, 1, 1, 1, 1/2]
#    compare_win(w, bm, 1)


def test_0_order_minority():
    # Where overlap is less than hop size

    # Small overlaps

    # Even length, periodic
    w = [0, 1, 1, 1, 1, 1]
    bm = borss_martin(6, 1, 0, False)
    compare_win(w, bm, 1)

    w = [0, 1/2, 1, 1, 1, 1/2]
    bm = borss_martin(6, 2, 0, False)
    compare_win(w, bm, 2)

    w = [0, 1/3, 2/3, 1, 2/3, 1/3]
    bm = borss_martin(6, 3, 0, False)
    compare_win(w, bm, 3)

    # Odd length, periodic
    bm = borss_martin(7, 1, 0, False)
    w = [0, 1, 1, 1, 1, 1, 1]
    compare_win(w, bm, 1)

    bm = borss_martin(7, 2, 0, False)
    w = [0, 1/2, 1, 1, 1, 1, 1/2]
    compare_win(w, bm, 2)

    bm = borss_martin(7, 3, 0, False)
    w = [0, 1/3, 2/3, 1, 1, 2/3, 1/3]
    compare_win(w, bm, 3)

    # Even length, symmetrical
    bm = borss_martin(6, 2, 0, True)
    w = [0, 1, 1, 1, 1, 0]
    compare_win(w, bm, 2)

    bm = borss_martin(6, 3, 0, True)
    w = [0, 1/2, 1, 1, 1/2, 0]
    compare_win(w, bm, 3)

    bm = borss_martin(7, 2, 0, True)
    w = [0, 1, 1, 1, 1, 1, 0]
    compare_win(w, bm, 2)

    bm = borss_martin(7, 3, 0, True)
    w = [0, 1/2, 1, 1, 1, 1/2, 0]
    compare_win(w, bm, 3)

    # Periodic is a truncated version of symmetrical:
    # Even periodic from odd symmetrical
    bm = borss_martin(8, 3, 0, False)
    w = array([0, 1, 2, 3, 3, 3, 2, 1]) / 3
    compare_win(w, bm, 3)

    bm = borss_martin(9, 4, 0, True)
    w = array([0, 1, 2, 3, 3, 3, 2, 1, 0]) / 3
    compare_win(w, bm, 4)

    # Odd periodic from even symmetrical
    bm = borss_martin(9, 3, 0, False)
    w = array([0, 1, 2, 3, 3, 3, 3, 2, 1]) / 3
    compare_win(w, bm, 3)

    bm = borss_martin(10, 4, 0, True)
    w = array([0, 1, 2, 3, 3, 3, 3, 2, 1, 0]) / 3
    compare_win(w, bm, 4)


def test_0_order_majority():
    # Where overlap is greater than hop size
    w = [0, 1/2, 1, 1, 1, 1/2]
    bm = borss_martin(6, 4, 0, False)
    compare_win(bm, w, 4)  # CHAHNGED


def test_1_order_tukey():
    # "the 1st-order window w1(t) is a Tukey window for TS > TL"
    bm = borss_martin(9, 3, 1, sym=True)
    tu = signal.tukey(9, alpha=0.5, sym=True)
    w = array([0, 1, 2, 2, 2, 2, 2, 1, 0]) / 2
    compare_win(bm, w, 3)
    compare_win(bm, tu, 3)

    # alpha = 0.75 -> 0.75/2 = 0.375 overlap = 3 samples overlap
    bm = borss_martin(8, 3, 1, sym=False)
    tu = signal.tukey(8, alpha=.75, sym=False)
    w = array([0, 1, 3, 4, 4, 4, 3, 1]) / 4
    compare_win(bm, w, 3)
    compare_win(bm, tu, 3)

    bm = borss_martin(9, 4, 1, sym=True)
    tu = signal.tukey(9, alpha=.75, sym=True)
    w = array([0, 1, 3, 4, 4, 4, 3, 1, 0]) / 4
    compare_win(bm, w, 4)
    compare_win(bm, tu, 4)


def test_1_order_hann():
    # "the 1st-order window w1(t) is ... a Hann window for TS = TL"
    bm = borss_martin(6, 3, 1, False)
    ha = signal.hann(6, False)
    compare_win(ha, bm, 3)

    bm = borss_martin(7, 4, 1, True)
    ha = signal.hann(7, True)
    compare_win(ha, bm, 4)

    bm = borss_martin(71, 36, 1, True)
    ha = signal.hann(71, True)
    compare_win(ha, bm, 36)

#    # signal.hann(6, True) is not COLA for 3 overlap, so shouldn't be equal
#    bm = borss_martin(6, 3, 1, True)
#    ha = signal.hann(6, True)
#    compare_win(ha, bm, 3)


def test_2_order():
    pass
    # TODO: add tests


def test_invalid():
    # Invalid lengths
    assert_raises(ValueError, borss_martin, -6, 3)
    assert_raises(ValueError, borss_martin, 6.3, 3)
    assert_raises(ValueError, borss_martin, 'spam', 3)

    # Invalid orders
    assert_raises(ValueError, borss_martin, 6, 3, -2)
    assert_raises(ValueError, borss_martin, 6, 3, 2.5)
    assert_raises(TypeError, borss_martin, 6, 3, 'egg')

    # Invalid overlaps
    assert_raises(ValueError, borss_martin, 6, 6)
    assert_raises(ValueError, borss_martin, 6, 8)
    assert_raises(ValueError, borss_martin, 6, -3)
    assert_raises(ValueError, borss_martin, 6, 3.2)
    assert_raises(TypeError, borss_martin, 6, 'bacon')


# TODO: add tests for slope of fall-off rate


if __name__ == "__main__":
    pytest.main(['--tb=short', __file__])

    # reproduce n=3 50% overlap

    M = 1000
    L = int(0.5 * M)
    N = 3
    import matplotlib.pyplot as plt
    from scipy.fftpack import fft, fftshift

    bm = borss_martin(M, L, N, False)

    assert check_COLA(bm, M, L)

    t = linspace(-1/2, 1/2, M)
    plt.figure('time', figsize=(10.1, 7.7245))
    plt.plot(t, bm, lw=2)
    plt.xlim(-5/8, 5/8)
    plt.ylim(-0.5, 1.5)
    plt.xticks((-1/2, -1/4, 0, 1/4, 1/2))
    plt.grid(True)

    plt.figure('freq', figsize=(9.675, 7.375))
    window = bm
    A = fft(window, 100001) / (len(window)/2.0)
    freq = np.linspace(-0.5, 0.5, len(A))
    freq = np.linspace(-M/2, M/2, len(A))
    response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
    plt.plot(freq, response, lw=2)
    plt.axis([-10, 10, -105, 5])
    plt.xticks((-10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10))
    plt.grid(True)

    # Confirmed that all graphs match the paper, both time and frequency

Fig. 3 n=3 50 pct overlap freq.png

Fig. 3 n=3 50 pct overlap time.png