20+ functions for generating different waveforms

Question

I added a guard clause to the functions in scipy.signal.windows, but the way they are currently written means the same 11 lines are now repeated in every function.

Each function has a unique core and a unique docstring, then some common guards and preconditioning around it:

if int(M) != M or M < 0:
    raise ValueError('Window length M must be a non-negative integer')
if M == 0:
    return np.array([])
if M == 1:
    return np.ones(1, 'd')
odd = M % 2
if not sym and not odd:
    M = M + 1

...

if not sym and not odd:
    w = w[:-1]
return w   

So I thought I'd factor out the common code to make it more DRY. However I can't find a really elegant way to do it:

• Moving the redundant code into helper functions doesn't work because the early exits don't work nested inside another function. Adding more code to catch the early exits kind of defeats the purpose.
• Wrapping them in a decorator changes the function signature from blackman(M, sym=True) to blackman(M, *args), since some window functions have different numbers of arguments than others.
• Splitting the unique guts of each function out into their own _core functions would work, but seems a little mangled, with the actual code being separated from the docstrings.

The complete file is here. Removed parts to make it less lengthy, while still showing some docstrings and functions with different numbers of arguments:

"""The suite of window functions."""
from __future__ import division, print_function, absolute_import

import warnings

import numpy as np
from scipy import fftpack, linalg, special
from scipy._lib.six import string_types

__all__ = ['boxcar', 'triang', 'parzen', 'bohman', 'blackman', 'nuttall',
           'blackmanharris', 'flattop', 'bartlett', 'hanning', 'barthann',
           'hamming', 'kaiser', 'gaussian', 'general_gaussian', 'chebwin',
           'slepian', 'cosine', 'hann', 'exponential', 'tukey', 'get_window']
...    

def triang(M, sym=True):
    """Return a triangular window.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero or less, an empty
        array is returned.
    sym : bool, optional
        When True (default), 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
        does not appear if `M` is even and `sym` is True).

    Examples
    --------
    Plot the window and its frequency response:

    >>> from scipy import signal
    >>> from scipy.fftpack import fft, fftshift
    >>> import matplotlib.pyplot as plt

    >>> window = signal.triang(51)
    >>> plt.plot(window)
    >>> plt.title("Triangular window")
    >>> plt.ylabel("Amplitude")
    >>> plt.xlabel("Sample")

    >>> plt.figure()
    >>> A = fft(window, 2048) / (len(window)/2.0)
    >>> freq = np.linspace(-0.5, 0.5, len(A))
    >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
    >>> plt.plot(freq, response)
    >>> plt.axis([-0.5, 0.5, -120, 0])
    >>> plt.title("Frequency response of the triangular window")
    >>> plt.ylabel("Normalized magnitude [dB]")
    >>> plt.xlabel("Normalized frequency [cycles per sample]")

    """
    if int(M) != M or M < 0:
        raise ValueError('Window length M must be a non-negative integer')
    if M == 0:
        return np.array([])
    if M == 1:
        return np.ones(1, 'd')
    odd = M % 2
    if not sym and not odd:
        M = M + 1
    n = np.arange(1, (M + 1) // 2 + 1)
    if M % 2 == 0:
        w = (2 * n - 1.0) / M
        w = np.r_[w, w[::-1]]
    else:
        w = 2 * n / (M + 1.0)
        w = np.r_[w, w[-2::-1]]

    if not sym and not odd:
        w = w[:-1]
    return w


def parzen(M, sym=True):
    """Return a Parzen window.
...    
    """
    if int(M) != M or M < 0:
        raise ValueError('Window length M must be a non-negative integer')
    if M == 0:
        return np.array([])
    if M == 1:
        return np.ones(1, 'd')
    odd = M % 2
    if not sym and not odd:
        M = M + 1
    n = np.arange(-(M - 1) / 2.0, (M - 1) / 2.0 + 0.5, 1.0)
    na = np.extract(n < -(M - 1) / 4.0, n)
    nb = np.extract(abs(n) <= (M - 1) / 4.0, n)
    wa = 2 * (1 - np.abs(na) / (M / 2.0)) ** 3.0
    wb = (1 - 6 * (np.abs(nb) / (M / 2.0)) ** 2.0 +
          6 * (np.abs(nb) / (M / 2.0)) ** 3.0)
    w = np.r_[wa, wb, wa[::-1]]
    if not sym and not odd:
        w = w[:-1]
    return w

...


def tukey(M, alpha=0.5, sym=True):
    r"""Return a Tukey window, also known as a tapered cosine window.
... 
    """
    if int(M) != M or M < 0:
        raise ValueError('Window length M must be a non-negative integer')
    if M == 0:
        return np.array([])
    if M == 1:
        return np.ones(1, 'd')

    if alpha <= 0:
        return np.ones(M, 'd')
    elif alpha >= 1.0:
        return hann(M, sym=sym)

    odd = M % 2
    if not sym and not odd:
        M = M + 1

    n = np.arange(0, M)
    width = int(np.floor(alpha*(M-1)/2.0))
    n1 = n[0:width+1]
    n2 = n[width+1:M-width-1]
    n3 = n[M-width-1:]

    w1 = 0.5 * (1 + np.cos(np.pi * (-1 + 2.0*n1/alpha/(M-1))))
    w2 = np.ones(n2.shape)
    w3 = 0.5 * (1 + np.cos(np.pi * (-2.0/alpha + 1 + 2.0*n3/alpha/(M-1))))

    w = np.concatenate((w1, w2, w3))

    if not sym and not odd:
        w = w[:-1]
    return w

...    

def slepian(M, width, sym=True):
    """Return a digital Slepian (DPSS) window.
...
    """
    if int(M) != M or M < 0:
        raise ValueError('Window length M must be a non-negative integer')
    if M == 0:
        return np.array([])
    if M == 1:
        return np.ones(1, 'd')
    odd = M % 2
    if not sym and not odd:
        M = M + 1

    # our width is the full bandwidth
    width = width / 2
    # to match the old version
    width = width / 2
    m = np.arange(M, dtype='d')
    H = np.zeros((2, M))
    H[0, 1:] = m[1:] * (M - m[1:]) / 2
    H[1, :] = ((M - 1 - 2 * m) / 2)**2 * np.cos(2 * np.pi * width)

    _, win = linalg.eig_banded(H, select='i', select_range=(M-1, M-1))
    win = win.ravel() / win.max()

    if not sym and not odd:
        win = win[:-1]
    return win
...    

def exponential(M, center=None, tau=1., sym=True):
    r"""Return an exponential (or Poisson) window.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero or less, an empty
        array is returned.
    center : float, optional
        Parameter defining the center location of the window function.
        The default value if not given is ``center = (M-1) / 2``.  This
        parameter must take its default value for symmetric windows.
    tau : float, optional
        Parameter defining the decay.  For ``center = 0`` use
        ``tau = -(M-1) / ln(x)`` if ``x`` is the fraction of the window
        remaining at the end.
    sym : bool, optional
        When True (default), 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
        does not appear if `M` is even and `sym` is True).

    Notes
    -----
    The Exponential window is defined as

    .. math::  w(n) = e^{-|n-center| / \tau}

    References
    ----------
    S. Gade and H. Herlufsen, "Windows to FFT analysis (Part I)",
    Technical Review 3, Bruel & Kjaer, 1987.

    Examples
    --------
    Plot the symmetric window and its frequency response:

    >>> from scipy import signal
    >>> from scipy.fftpack import fft, fftshift
    >>> import matplotlib.pyplot as plt

    >>> M = 51
    >>> tau = 3.0
    >>> window = signal.exponential(M, tau=tau)
    >>> plt.plot(window)
    >>> plt.title("Exponential Window (tau=3.0)")
    >>> plt.ylabel("Amplitude")
    >>> plt.xlabel("Sample")

    >>> plt.figure()
    >>> A = fft(window, 2048) / (len(window)/2.0)
    >>> freq = np.linspace(-0.5, 0.5, len(A))
    >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
    >>> plt.plot(freq, response)
    >>> plt.axis([-0.5, 0.5, -35, 0])
    >>> plt.title("Frequency response of the Exponential window (tau=3.0)")
    >>> plt.ylabel("Normalized magnitude [dB]")
    >>> plt.xlabel("Normalized frequency [cycles per sample]")

    This function can also generate non-symmetric windows:

    >>> tau2 = -(M-1) / np.log(0.01)
    >>> window2 = signal.exponential(M, 0, tau2, False)
    >>> plt.figure()
    >>> plt.plot(window2)
    >>> plt.ylabel("Amplitude")
    >>> plt.xlabel("Sample")
    """
    if sym and center is not None:
        raise ValueError("If sym==True, center must be None.")
    if int(M) != M or M < 0:
        raise ValueError('Window length M must be a non-negative integer')
    if M == 0:
        return np.array([])
    if M == 1:
        return np.ones(1, 'd')
    odd = M % 2
    if not sym and not odd:
        M = M + 1
    if center is None:
        center = (M-1) / 2

    n = np.arange(0, M)
    w = np.exp(-np.abs(n-center) / tau)
    if not sym and not odd:
        w = w[:-1]

    return w

...

Answer 1

I do believe that the decorator way is the proper one. There are some tradeoff using it but you can easily overcome your main concern using functools.wraps: it will reuse the name, docstring and signature of the decorated function for the wrapper.

The wrapper within the decorator should be aware of both M and sym; this is where things can get tricky depending on how you want things to behave.

A first approach could be:

from functools import wraps


def argument_checker(func):
    @wraps(func)
    def wrapper(M, *args, sym=True, **kwargs):
        if int(M) != M or M < 0:
            raise ValueError('Window length M must be a non-negative integer')
        if M == 0:
            return np.array([])
        if M == 1:
            return np.ones(1, 'd')

        if not sym and not M % 2:
            return func(M + 1, *args, **kwargs)[:-1]
        else:
            return func(M, *args, **kwargs)
    return wrapper

I just modified the boilerplate code to group the two if not sym and not odd into a single if as there is func that abstract the window building. However this solution is Python 3 only as the signature of wrapper will raise a SyntaxError in Python 2, and as regard to your imports from __future__, I'm assuming Python 2 here. So you might want:

def argument_checker(func):
    @wraps(func)
    def wrapper(M, *args, **kwargs):
        sym = kwargs.pop('sym', True)
        ...

But this poses the issue of sym needing to be a keyword argument, and not a positional one. Again, this can be made more explicit using Python 3 syntax for keyword-only arguments:

def triang(M, *, sym=True):
    ...

or

def tukey(M, alpha=0.5, *, sym=True):
    ...

Which, when called naturally using triang(150, False), for instance, will raise TypeError: triang() takes 1 positional argument but 2 were given.

So I'd go for this second version and keep the current signature of the functions for now, but if you can use Python 3 instead, it will be made much more explicit. I would, however, remove the *args parameter from the wrapper to make it accept only keyword arguments.

But then, there is the case of tukey that add an other layer of check into the checks. To handle that properly, it is easier to split the decorator into two, more meaningful, tasks:

from functools import wraps


def argument_checker(func):
    @wraps(func)
    def wrapper(M, **kwargs):
        if int(M) != M or M < 0:
            raise ValueError('Window length M must be a non-negative integer')
        if M == 0:
            return np.array([])
        if M == 1:
            return np.ones(1, 'd')

        return func(M, **kwargs)
    return wrapper


def parity_handler(func):
    @wraps(func)
    def wrapper(M, **kwargs):
        sym = kwargs.get('sym', True)
        if not sym and not M % 2:
            return func(M + 1, **kwargs)[:-1]
        else:
            return func(M, **kwargs)
    return wrapper

So you can

@argument_checker
@parity_handler
def triang(M, sym=True):
    ...

But most importantly:

def alpha_handler(func):
    @wraps(func)
    def wrapper(M, **kwargs):
        alpha = kwargs.get('alpha', 0.5)
        if alpha <= 0:
            return np.ones(M, 'd')
        elif alpha >= 1.0:
            return hann(M, sym=kwargs.get('sym', True))
        return func(M, **kwargs)
    return wrapper

@argument_checker
@alpha_handler
@parity_handler
def tukey(...

And the whole module could become:

"""The suite of window functions."""
from __future__ import division, print_function, absolute_import

import warnings
from functools import wraps

import numpy as np
from scipy import fftpack, linalg, special
from scipy._lib.six import string_types

__all__ = ['boxcar', 'triang', 'parzen', 'bohman', 'blackman', 'nuttall',
           'blackmanharris', 'flattop', 'bartlett', 'hanning', 'barthann',
           'hamming', 'kaiser', 'gaussian', 'general_gaussian', 'chebwin',
           'slepian', 'cosine', 'hann', 'exponential', 'tukey', 'get_window']
...    


def argument_checker(func):
    @wraps(func)
    def wrapper(M, **kwargs):
        if int(M) != M or M < 0:
            raise ValueError('Window length M must be a non-negative integer')
        if M == 0:
            return np.array([])
        if M == 1:
            return np.ones(1, 'd')

        return func(M, **kwargs)
    return wrapper


def parity_handler(func):
    @wraps(func)
    def wrapper(M, **kwargs):
        sym = kwargs.get('sym', True)
        if not sym and not M % 2:
            return func(M + 1, **kwargs)[:-1]
        else:
            return func(M, **kwargs)
    return wrapper


@argument_checker
@parity_handler
def triang(M, sym=True):
    """Return a triangular window.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero or less, an empty
        array is returned.
    sym : bool, optional
        When True (default), 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
        does not appear if `M` is even and `sym` is True).

    Examples
    --------
    Plot the window and its frequency response:

    >>> from scipy import signal
    >>> from scipy.fftpack import fft, fftshift
    >>> import matplotlib.pyplot as plt

    >>> window = signal.triang(51)
    >>> plt.plot(window)
    >>> plt.title("Triangular window")
    >>> plt.ylabel("Amplitude")
    >>> plt.xlabel("Sample")

    >>> plt.figure()
    >>> A = fft(window, 2048) / (len(window)/2.0)
    >>> freq = np.linspace(-0.5, 0.5, len(A))
    >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
    >>> plt.plot(freq, response)
    >>> plt.axis([-0.5, 0.5, -120, 0])
    >>> plt.title("Frequency response of the triangular window")
    >>> plt.ylabel("Normalized magnitude [dB]")
    >>> plt.xlabel("Normalized frequency [cycles per sample]")

    """
    n = np.arange(1, (M + 1) // 2 + 1)
    if M % 2 == 0:
        w = (2 * n - 1.0) / M
        w = np.r_[w, w[::-1]]
    else:
        w = 2 * n / (M + 1.0)
        w = np.r_[w, w[-2::-1]]
    return w


@argument_checker
@parity_handler
def parzen(M, sym=True):
    """Return a Parzen window.
...    
    """
    n = np.arange(-(M - 1) / 2.0, (M - 1) / 2.0 + 0.5, 1.0)
    na = np.extract(n < -(M - 1) / 4.0, n)
    nb = np.extract(abs(n) <= (M - 1) / 4.0, n)
    wa = 2 * (1 - np.abs(na) / (M / 2.0)) ** 3.0
    wb = (1 - 6 * (np.abs(nb) / (M / 2.0)) ** 2.0 +
          6 * (np.abs(nb) / (M / 2.0)) ** 3.0)
    w = np.r_[wa, wb, wa[::-1]]
    return w

...


def alpha_handler(func):
    @wraps(func)
    def wrapper(M, **kwargs):
        alpha = kwargs.get('alpha', 0.5)
        if alpha <= 0:
            return np.ones(M, 'd')
        elif alpha >= 1.0:
            return hann(M, sym=kwargs.get('sym', True))
        return func(M, **kwargs)
    return wrapper

@argument_checker
@alpha_handler
@parity_handler
def tukey(M, alpha=0.5, sym=True):
    r"""Return a Tukey window, also known as a tapered cosine window.
... 
    """    
    n = np.arange(0, M)
    width = int(np.floor(alpha*(M-1)/2.0))
    n1 = n[0:width+1]
    n2 = n[width+1:M-width-1]
    n3 = n[M-width-1:]

    w1 = 0.5 * (1 + np.cos(np.pi * (-1 + 2.0*n1/alpha/(M-1))))
    w2 = np.ones(n2.shape)
    w3 = 0.5 * (1 + np.cos(np.pi * (-2.0/alpha + 1 + 2.0*n3/alpha/(M-1))))

    w = np.concatenate((w1, w2, w3))

    return w

...    

@argument_checker
@parity_handler
def slepian(M, width, sym=True):
    """Return a digital Slepian (DPSS) window.
...
    """    
    # our width is the full bandwidth
    width = width / 2
    # to match the old version
    width = width / 2
    m = np.arange(M, dtype='d')
    H = np.zeros((2, M))
    H[0, 1:] = m[1:] * (M - m[1:]) / 2
    H[1, :] = ((M - 1 - 2 * m) / 2)**2 * np.cos(2 * np.pi * width)

    _, win = linalg.eig_banded(H, select='i', select_range=(M-1, M-1))
    win = win.ravel() / win.max()

    return win
...    

@argument_checker
@parity_handler
def exponential(M, center=None, tau=1., sym=True):
    r"""Return an exponential (or Poisson) window.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero or less, an empty
        array is returned.
    center : float, optional
        Parameter defining the center location of the window function.
        The default value if not given is ``center = (M-1) / 2``.  This
        parameter must take its default value for symmetric windows.
    tau : float, optional
        Parameter defining the decay.  For ``center = 0`` use
        ``tau = -(M-1) / ln(x)`` if ``x`` is the fraction of the window
        remaining at the end.
    sym : bool, optional
        When True (default), 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
        does not appear if `M` is even and `sym` is True).

    Notes
    -----
    The Exponential window is defined as

    .. math::  w(n) = e^{-|n-center| / \tau}

    References
    ----------
    S. Gade and H. Herlufsen, "Windows to FFT analysis (Part I)",
    Technical Review 3, Bruel & Kjaer, 1987.

    Examples
    --------
    Plot the symmetric window and its frequency response:

    >>> from scipy import signal
    >>> from scipy.fftpack import fft, fftshift
    >>> import matplotlib.pyplot as plt

    >>> M = 51
    >>> tau = 3.0
    >>> window = signal.exponential(M, tau=tau)
    >>> plt.plot(window)
    >>> plt.title("Exponential Window (tau=3.0)")
    >>> plt.ylabel("Amplitude")
    >>> plt.xlabel("Sample")

    >>> plt.figure()
    >>> A = fft(window, 2048) / (len(window)/2.0)
    >>> freq = np.linspace(-0.5, 0.5, len(A))
    >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
    >>> plt.plot(freq, response)
    >>> plt.axis([-0.5, 0.5, -35, 0])
    >>> plt.title("Frequency response of the Exponential window (tau=3.0)")
    >>> plt.ylabel("Normalized magnitude [dB]")
    >>> plt.xlabel("Normalized frequency [cycles per sample]")

    This function can also generate non-symmetric windows:

    >>> tau2 = -(M-1) / np.log(0.01)
    >>> window2 = signal.exponential(M, 0, tau2, False)
    >>> plt.figure()
    >>> plt.plot(window2)
    >>> plt.ylabel("Amplitude")
    >>> plt.xlabel("Sample")
    """
    if sym and center is not None:
        raise ValueError("If sym==True, center must be None.")
    if center is None:
        center = (M-1) / 2

    n = np.arange(0, M)
    w = np.exp(-np.abs(n-center) / tau)

    return w

...

You can also blend the same approach than for tukey into exponential if you want to check for if sym and center is not None: earlier.

Answer 2

I'm very new to thinking functional-programming-style, but I'm going to give it a shot.

The way I see it, each of these functions looks like

postprocess( core( preprocess( M ) ) )

We don't actually want to make new _core functions, however, so perhaps a little bit of sufficiently generic and adaptable boilerplate is acceptable. Unfortunately, this will come at the cost of code complexity, but hopefully, comments will illustrate intent, and make code maintenance easier in future.

To start off with, we could have

def _preprocess(M, sym):
    if int(M) != M or M < 0:
        raise ValueError('Window length M must be a non-negative integer')
    if M == 0:
        raise ReturnException(np.array([]))
    if M == 1:
        return ReturnException(np.ones(1, 'd'))
    odd = M % 2
    if not sym and not odd:
        M = M + 1
    return M

so that the starting boiler plate is reduced to something like:

def triang(M, sym=True):
    try:
        M = _preprocess(M, sym)
    except ReturnException as r:
        return r.value

ReturnException is something I'm defining a la "MyError" defined in the Python documentation.

But of course, all of this will not work for core-function-specific boilerplate. But coming to that later...

Then, we have a post-processing function:

def _postprocess(w, sym):
    if not sym and not odd:
        w = w[:-1]
    return w

and in each of the core functions, simply:

return _postprocess(w, sym)

If you really want to explicitly have function-specific steps happening at some stage during the pre-processing process, then you might want to consider adding a hook for that too (for example, in the tukey function in your question)

def _preprocess(M, sym, core_specific_fn, *args, **kwargs):
        if int(M) != M or M < 0:
    raise ValueError('Window length M must be a non-negative integer')
if M == 0:
    return np.array([])
if M == 1:
    return np.ones(1, 'd')

M, sym = core_specific_fn(M, sym, *args, **kwargs)

odd = M % 2
if not sym and not odd:
    M = M + 1

with the core_specific_fn being defined in tukey itself, and raising a return exception too. Preprocess will be called after this "sub-function" has been defined.

But in the tukey case, this is really not necessary. Just in case it is necessary somewhere else though...

But I should add that all this might not be worth the extra complexity it adds to the code. Especially if things like the additional hook from tukey is a common occurrence across the twenty functions.