Source code for spectrum.window

"""
.. topic:: This module contains tapering windows utilities.

    .. autosummary::

        Window
        window_visu
        create_window
        window_hann
        window_hamming
        window_bartlett
        window_bartlett_hann
        window_blackman
        window_blackman_harris
        window_blackman_nuttall
        window_bohman
        window_chebwin
        window_cosine
        window_flattop
        window_gaussian
        window_hamming
        window_hann
        window_kaiser
        window_lanczos
        window_nuttall
        window_parzen
        window_taylor
        window_tukey


    .. codeauthor:: Thomas Cokelaer 2011

    :References: See [Nuttall]_, [Marple]_, [Harris]_
"""
import numpy as np
from numpy import pi, cos, arange, array, sin, exp, sinc, linspace, \
    sqrt, ones, log, array, prod

from spectrum import tools as stools

#__all__ = ["window_names","Window","create_window"]


window_names = {
                 'bartlett_hann': 'window_bartlett_hann',
                 'blackman_harris': 'window_blackman_harris',
                 'blackman_nuttall': 'window_blackman_nuttall',
                 'bohman':'window_bohman',
                 'blackman': 'window_blackman',
                 'chebwin':'window_chebwin',
                 'gaussian': 'window_gaussian',
                 'hamming':'window_hamming',
                 'kaiser': 'window_kaiser',
                 'lanczos': 'window_lanczos',
                 'sinc': 'window_lanczos',
                 'poisson': 'window_poisson',
                 'tukey':'window_tukey',
                 'nuttall':'window_nuttall',
                 'parzen':'window_parzen',
                 'flattop': 'window_flattop',
                 'riesz':'window_riesz',
                 'riemann':'window_riemann',
                 'hann':'window_hann',
                 'hanning':'window_hann',
                 'poisson_hanning': 'window_poisson_hanning',
                 'rectangular':'window_rectangle',
                 'rectangle': 'window_rectangle',
                 'bartlett':'window_bartlett',
                 'triangular':     'window_bartlett',
                 'cosine': 'window_cosine',
                 'sine':   'window_cosine',
                 'cauchy':'window_cauchy',
                 'taylor':'window_taylor',
                 }


[docs]class Window(object): r"""Window tapering object This class provides utilities to manipulate tapering windows. Plotting functions allows to visualise the time and frequency response. It is also possible to retrieve relevant quantities such as the equivalent noise band width. The following examples illustrates the usage. First, we create the window by providing a name and a size:: from spectrum import Window w = Window(64, 'hamming') The window has been computed and the data is stored in:: w.data This object contains plotting methods so that you can see the time or frequency response. .. plot:: :include-source: :width: 80% from spectrum.window import Window w = Window(64, 'hamming') w.plot_frequencies() Some windows may accept optional arguments. For instance, :func:`window_blackman` accepts an optional argument called :math:`\alpha` as well as :class:`Window`. Indeed, we use the factory :func:`~spectrum.window.create_window`, which manage all the optional arguments. So you can write:: w = Window(64, 'blackman', alpha=1) .. seealso:: :func:`create_window`. .. rubric:: Constructor: """ def __init__(self, N, name=None, norm=True, **kargs): """Create a tapering window object :param N: the window length :param name: the type of window, e.g., 'Hann' :param norm: normalise the window in frequency domain (for plotting) :param kargs: any of :func:`create_window` valid optional arguments. .. rubric:: Attributes: * data: time series data * frequencies: getter to the frequency series * response: getter to the PSD * enbw: getter to the Equivalent noise band width. """ #input attributse assert N>0 , "First argument must be positive" if name is None or name not in list(window_names.keys()): raise ValueError("second argument must be a valid window name %s" % list(window_names.keys())) self.__N = N self.__name = name self.__norm = norm #compute the window data and Fourier domain data self.__data = create_window(N, name, **kargs) self.__frequencies = None self.__response = None # equivalent noise BW and other quantities. self.__enbw = enbw(self.data) def _getNorm(self): return self.__norm norm = property(fget=_getNorm, doc="Getter of the normalisation flag (True by default)") def _getN(self): return self.__N N = property(fget=_getN, doc="Getter for the window length") def _getENBW(self): return self.__enbw enbw = property(fget=_getENBW, doc="getter for the equivalent noise band\ width. See :func:`enbw` function") def _getMeanSquare(self): return np.sum(self.data**2)/self.N mean_square = property(fget=_getMeanSquare, doc="returns :math:`\frac{w^2}{N}`") def _getName(self): return self.__name name = property(fget=_getName, doc="Getter for the window name") def _getData(self): return self.__data data = property(fget=_getData, doc="Getter for the window values (in time)") def _getF(self): if self.__response is None: self.compute_response() self.__frequencies = linspace(-0.5, 0.5, len(self.__response)) return self.__frequencies frequencies = property(fget=_getF, doc="Getter for the frequency array") def _getResponse(self): if self.__response is None: self.compute_response() return self.__response response = property(fget=_getResponse, doc="Getter for the frequency \ response. See :meth:`compute_response`")
[docs] def compute_response(self, **kargs): """Compute the window data frequency response :param norm: True by default. normalised the frequency data. :param int NFFT: total length of the final data sets( 2048 by default. if less than data length, then NFFT is set to the data length*2). The response is stored in :attr:`response`. .. note:: Units are dB (20 log10) since we plot the frequency response) """ from numpy.fft import fft, fftshift norm = kargs.get('norm', self.norm) # do some padding. Default is max(2048, data.len*2) NFFT = kargs.get('NFFT', 2048) if NFFT < len(self.data): NFFT = self.data.size * 2 # compute the fft modulus A = fft(self.data, NFFT) mag = abs(fftshift(A)) # do we want to normalise the data if norm is True: mag = mag / max(mag) response = 20. * stools.log10(mag) # factor 20 we are looking at the response # not the powe #response = clip(response,mindB,100) self.__response = response
[docs] def plot_frequencies(self, mindB=None, maxdB=None, norm=True): """Plot the window in the frequency domain :param mindB: change the default lower y bound :param maxdB: change the default upper lower bound :param bool norm: if True, normalise the frequency response. .. plot:: :width: 80% :include-source: from spectrum.window import Window w = Window(64, name='hamming') w.plot_frequencies() """ from pylab import plot, title, xlim, grid, ylim, xlabel, ylabel # recompute the response self.compute_response(norm=norm) plot(self.frequencies, self.response) title("ENBW=%2.1f" % (self.enbw)) ylabel('Frequency response (dB)') xlabel('Fraction of sampling frequency') # define the plot limits xlim(-0.5, 0.5) y0, y1 = ylim() if mindB: y0 = mindB if maxdB is not None: y1 = maxdB else: y1 = max(self.response) ylim(y0, y1) grid(True)
[docs] def plot_window(self): """Plot the window in the time domain .. plot:: :width: 80% :include-source: from spectrum.window import Window w = Window(64, name='hamming') w.plot_window() """ from pylab import plot, xlim, grid, title, ylabel, axis x = linspace(0, 1, self.N) xlim(0, 1) plot(x, self.data) grid(True) title('%s Window (%s points)' % (self.name.capitalize(), self.N)) ylabel('Amplitude') axis([0, 1, 0, 1.1])
[docs] def plot_time_freq(self, mindB=-100, maxdB=None, norm=True, yaxis_label_position="right"): """Plotting method to plot both time and frequency domain results. See :meth:`plot_frequencies` for the optional arguments. .. plot:: :width: 80% :include-source: from spectrum.window import Window w = Window(64, name='hamming') w.plot_time_freq() """ from pylab import subplot, gca subplot(1, 2, 1) self.plot_window() subplot(1, 2, 2) self.plot_frequencies(mindB=mindB, maxdB=maxdB, norm=norm) if yaxis_label_position=="left": try: tight_layout() except: pass else: ax = gca() ax.yaxis.set_label_position("right")
[docs] def info(self): """Print object information such as length and name""" print(self)
def __str__(self): msg = 'Window object:\n' msg += 'Name: %s\n' % self.name msg += 'Length: %s\n' % self.N msg += 'NFFT (for the frequency response): %s\n' % self.response.size msg += 'ENBW=%s\n' % self.enbw return msg
[docs]def create_window(N, name=None, **kargs): r"""Returns the N-point window given a valid name :param int N: window size :param str name: window name (default is *rectangular*). Valid names are stored in :func:`~spectrum.window.window_names`. :param kargs: optional arguments are: * *beta*: argument of the :func:`window_kaiser` function (default is 8.6) * *attenuation*: argument of the :func:`window_chebwin` function (default is 50dB) * *alpha*: argument of the 1. :func:`window_gaussian` function (default is 2.5) 2. :func:`window_blackman` function (default is 0.16) 3. :func:`window_poisson` function (default is 2) 4. :func:`window_cauchy` function (default is 3) * *mode*: argument :func:`window_flattop` function (default is *symmetric*, can be *periodic*) * *r*: argument of the :func:`window_tukey` function (default is 0.5). The following windows have been simply wrapped from existing librairies like NumPy: * **Rectangular**: :func:`window_rectangle`, * **Bartlett** or Triangular: see :func:`window_bartlett`, * **Hanning** or Hann: see :func:`window_hann`, * **Hamming**: see :func:`window_hamming`, * **Kaiser**: see :func:`window_kaiser`, * **chebwin**: see :func:`window_chebwin`. The following windows have been implemented from scratch: * **Blackman**: See :func:`window_blackman` * **Bartlett-Hann** : see :func:`window_bartlett_hann` * **cosine or sine**: see :func:`window_cosine` * **gaussian**: see :func:`window_gaussian` * **Bohman**: see :func:`window_bohman` * **Lanczos or sinc**: see :func:`window_lanczos` * **Blackman Harris**: see :func:`window_blackman_harris` * **Blackman Nuttall**: see :func:`window_blackman_nuttall` * **Nuttall**: see :func:`window_nuttall` * **Tukey**: see :func:`window_tukey` * **Parzen**: see :func:`window_parzen` * **Flattop**: see :func:`window_flattop` * **Riesz**: see :func:`window_riesz` * **Riemann**: see :func:`window_riemann` * **Poisson**: see :func:`window_poisson` * **Poisson-Hanning**: see :func:`window_poisson_hanning` .. todo:: on request taylor, potter, Bessel, expo, rife-vincent, Kaiser-Bessel derived (KBD) .. plot:: :width: 80% :include-source: from pylab import plot, legend from spectrum import create_window data = create_window(51, 'hamming') plot(data, label='hamming') data = create_window(51, 'kaiser') plot(data, label='kaiser') legend() .. plot:: :width: 80% :include-source: from pylab import plot, log10, linspace, fft, clip from spectrum import create_window, fftshift A = fft(create_window(51, 'hamming'), 2048) / 25.5 mag = abs(fftshift(A)) freq = linspace(-0.5,0.5,len(A)) response = 20 * log10(mag) mindB = -60 response = clip(response,mindB,100) plot(freq, response) .. seealso:: :func:`window_visu`, :func:`Window`, :mod:`spectrum.dpss` """ if name is None: name = 'rectangle' name = name.lower() assert name in list(window_names.keys()), \ """window name %s not implemented or incorrect. Try to use one of %s"""\ % (name, window_names) # create the function name f = eval(window_names[name]) windows_with_parameters = \ {'kaiser': {'beta': eval(window_names['kaiser']).__defaults__[0]}, 'blackman': {'alpha': eval(window_names['blackman']).__defaults__[0]}, 'cauchy': {'alpha': eval(window_names['cauchy']).__defaults__[0]}, 'flattop': {'mode': eval(window_names['flattop']).__defaults__[0]}, 'gaussian': {'alpha': eval(window_names['gaussian']).__defaults__[0]}, 'chebwin': {'attenuation':eval(window_names['chebwin']).__defaults__[0]}, 'tukey': {'r':eval(window_names['tukey']).__defaults__[0]}, 'poisson': {'alpha': eval(window_names['poisson']).__defaults__[0]}, 'poisson_hanning': {'alpha': eval(window_names['poisson_hanning']).__defaults__[0]}, 'taylor': {'nbar': eval(window_names['taylor']).__defaults__[0], 'sll': eval(window_names['taylor']).__defaults__[0]}, } if name not in list(windows_with_parameters.keys()): if len(kargs) == 0: # no parameters, so we directly call the function w = f(N) else: raise ValueError(""" Parameters do not match any of the window. The window provided do not expect any parameters. Try to remove the parameters""") elif name in list(windows_with_parameters.keys()): # user optional parameters are provided, scan them: dargs = {} for arg in list(kargs.keys()): # check that the parameters are valid, and obtain the default value try: default = windows_with_parameters[name][arg] except: raise ValueError(""" Invalid optional argument (%s) for %s window. Valid optional arguments are (%s)""" % \ (arg, name, list(windows_with_parameters[name].keys()))) # add the user parameter to the list of parameters dargs[arg] = kargs.get(arg, default) # call the proper function with the optional arguments w = f(N, **dargs) return w
[docs]def enbw(data): r"""Computes the equivalent noise bandwidth .. math:: ENBW = N \frac{\sum_{n=1}^{N} w_n^2}{\left(\sum_{n=1}^{N} w_n \right)^2} .. doctest:: >>> from spectrum import create_window, enbw >>> w = create_window(64, 'rectangular') >>> enbw(w) 1.0 The following table contains the ENBW values for some of the implemented windows in this module (with N=16384). They have been double checked against litterature (Source: [Harris]_, [Marple]_). If not present, it means that it has not been checked. =================== ============ ============= name ENBW litterature =================== ============ ============= rectangular 1. 1. triangle 1.3334 1.33 Hann 1.5001 1.5 Hamming 1.3629 1.36 blackman 1.7268 1.73 kaiser 1.7 blackmanharris,4 2.004 2. riesz 1.2000 1.2 riemann 1.32 1.3 parzen 1.917 1.92 tukey 0.25 1.102 1.1 bohman 1.7858 1.79 poisson 2 1.3130 1.3 hanningpoisson 0.5 1.609 1.61 cauchy 1.489 1.48 lanczos 1.3 =================== ============ ============= """ N = len(data) return N * np.sum(data**2) / np.sum(data)**2
def _kaiser(n, beta): """Independant Kaiser window For the definition of the Kaiser window, see A. V. Oppenheim & R. W. Schafer, "Discrete-Time Signal Processing". The continuous version of width n centered about x=0 is: .. note:: 2 times slower than scipy.kaiser """ from scipy.special import iv as besselI m = n - 1 k = arange(0, m) k = 2. * beta / m * sqrt (k * (m - k)) w = besselI (0, k) / besselI (0, beta) return w
[docs]def window_visu(N=51, name='hamming', **kargs): """A Window visualisation tool :param N: length of the window :param name: name of the window :param NFFT: padding used by the FFT :param mindB: the minimum frequency power in dB :param maxdB: the maximum frequency power in dB :param kargs: optional arguments passed to :func:`create_window` This function plot the window shape and its equivalent in the Fourier domain. .. plot:: :width: 80% :include-source: from spectrum import window_visu window_visu(64, 'kaiser', beta=8.) """ # get the default parameters mindB = kargs.pop('mindB', -100) maxdB = kargs.pop('maxdB', None) norm = kargs.pop('norm', True) # create a window object w = Window(N, name, **kargs) # plot the time and frequency windows w.plot_time_freq(mindB=mindB, maxdB=maxdB, norm=norm)
[docs]def window_rectangle(N): r"""Kaiser window :param N: window length .. plot:: :width: 80% :include-source: from spectrum import window_visu window_visu(64, 'rectangle') """ return ones(N)
[docs]def window_kaiser(N, beta=8.6, method='numpy'): r"""Kaiser window :param N: window length :param beta: kaiser parameter (default is 8.6) To obtain a Kaiser window that designs an FIR filter with sidelobe attenuation of :math:`\alpha` dB, use the following :math:`\beta` where :math:`\beta = \pi \alpha`. .. math:: w_n = \frac{I_0\left(\pi\alpha\sqrt{1-\left(\frac{2n}{M}-1\right)^2}\right)} {I_0(\pi \alpha)} where * :math:`I_0` is the zeroth order Modified Bessel function of the first kind. * :math:`\alpha` is a real number that determines the shape of the window. It determines the trade-off between main-lobe width and side lobe level. * the length of the sequence is N=M+1. The Kaiser window can approximate many other windows by varying the :math:`\beta` parameter: ===== ======================== beta Window shape ===== ======================== 0 Rectangular 5 Similar to a Hamming 6 Similar to a Hanning 8.6 Similar to a Blackman ===== ======================== .. plot:: :width: 80% :include-source: from pylab import plot, legend, xlim from spectrum import window_kaiser N = 64 for beta in [1,2,4,8,16]: plot(window_kaiser(N, beta), label='beta='+str(beta)) xlim(0,N) legend() .. plot:: :width: 80% :include-source: from spectrum import window_visu window_visu(64, 'kaiser', beta=8.) .. seealso:: numpy.kaiser, :func:`spectrum.window.create_window` """ if N == 1: return ones(1) if method == 'numpy': from numpy import kaiser return kaiser(N, beta) else: return _kaiser(N, beta)
[docs]def window_blackman(N, alpha=0.16): r"""Blackman window :param N: window length .. math:: a_0 - a_1 \cos(\frac{2\pi n}{N-1}) +a_2 \cos(\frac{4\pi n }{N-1}) with .. math:: a_0 = (1-\alpha)/2, a_1=0.5, a_2=\alpha/2 \rm{\;and\; \alpha}=0.16 When :math:`\alpha=0.16`, this is the unqualified Blackman window with :math:`a_0=0.48` and :math:`a_2=0.08`. .. plot:: :width: 80% :include-source: from spectrum import window_visu window_visu(64, 'blackman') .. note:: Although Numpy implements a blackman window for :math:`\alpha=0.16`, this implementation is valid for any :math:`\alpha`. .. seealso:: numpy.blackman, :func:`create_window`, :class:`Window` """ a0 = (1. - alpha)/2. a1 = 0.5 a2 = alpha/2. if (N == 1): win = array([1.]) else: k = arange(0, N)/float(N-1.) win = a0 - a1 * cos (2 * pi * k) + a2 * cos (4 * pi * k) return win
#from numpy import blackman #return blackman(N)
[docs]def window_bartlett(N): r"""Bartlett window (wrapping of numpy.bartlett) also known as Fejer :param int N: window length The Bartlett window is defined as .. math:: w(n) = \frac{2}{N-1} \left( \frac{N-1}{2} - \left|n - \frac{N-1}{2}\right| \right) .. plot:: :width: 80% :include-source: from spectrum import window_visu window_visu(64, 'bartlett') .. seealso:: numpy.bartlett, :func:`create_window`, :class:`Window`. """ from numpy import bartlett return bartlett(N)
[docs]def window_hamming(N): r"""Hamming window :param N: window length The Hamming window is defined as .. math:: 0.54 -0.46 \cos\left(\frac{2\pi n}{N-1}\right) \qquad 0 \leq n \leq M-1 .. plot:: :width: 80% :include-source: from spectrum import window_visu window_visu(64, 'hamming') .. seealso:: numpy.hamming, :func:`create_window`, :class:`Window`. """ from numpy import hamming return hamming(N)
[docs]def window_hann(N): r"""Hann window (or Hanning). (wrapping of numpy.bartlett) :param int N: window length The Hanning window is also known as the Cosine Bell. Usually, it is called Hann window, to avoid confusion with the Hamming window. .. math:: w(n) = 0.5\left(1- \cos\left(\frac{2\pi n}{N-1}\right)\right) \qquad 0 \leq n \leq M-1 .. plot:: :width: 80% :include-source: from spectrum import window_visu window_visu(64, 'hanning') .. seealso:: numpy.hanning, :func:`create_window`, :class:`Window`. """ from numpy import hanning return hanning(N)
[docs]def window_gaussian(N, alpha=2.5): r"""Gaussian window :param N: window length .. math:: \exp^{-0.5 \left( \sigma\frac{n}{N/2} \right)^2} with :math:`\frac{N-1}{2}\leq n \leq \frac{N-1}{2}`. .. note:: N-1 is used to be in agreement with octave convention. The ENBW of 1.4 is also in agreement with [Harris]_ .. plot:: :width: 80% :include-source: from spectrum import window_visu window_visu(64, 'gaussian', alpha=2.5) .. seealso:: scipy.signal.gaussian, :func:`create_window` """ t = linspace(-(N-1)/2., (N-1)/2., N) #t = linspace(-(N)/2., (N)/2., N) w = exp(-0.5*(alpha * t/(N/2.))**2.) return w
[docs]def window_chebwin(N, attenuation=50): """Cheb window :param N: window length .. plot:: :width: 80% :include-source: from spectrum import window_visu window_visu(64, 'chebwin', attenuation=50) .. seealso:: scipy.signal.chebwin, :func:`create_window`, :class:`Window` """ import scipy.signal return scipy.signal.chebwin(N, attenuation)
[docs]def window_cosine(N): r"""Cosine tapering window also known as sine window. :param N: window length .. math:: w(n) = \cos\left(\frac{\pi n}{N-1} - \frac{\pi}{2}\right) = \sin \left(\frac{\pi n}{N-1}\right) .. plot:: :width: 80% :include-source: from spectrum import window_visu window_visu(64, 'cosine') .. seealso:: :func:`create_window`, :class:`Window` """ if N ==1: return ones(1) n = arange(0, N) win = sin(pi*n/(N-1.)) return win
[docs]def window_lanczos(N): r"""Lanczos window also known as sinc window. :param N: window length .. math:: w(n) = sinc \left( \frac{2n}{N-1} - 1 \right) .. plot:: :width: 80% :include-source: from spectrum import window_visu window_visu(64, 'lanczos') .. seealso:: :func:`create_window`, :class:`Window` """ if N ==1: return ones(1) n = linspace(-N/2., N/2., N) win = sinc(2*n/(N-1.)) return win
[docs]def window_bartlett_hann(N): r"""Bartlett-Hann window :param N: window length .. math:: w(n) = a_0 + a_1 \left| \frac{n}{N-1} -\frac{1}{2}\right| - a_2 \cos \left( \frac{2\pi n}{N-1} \right) with :math:`a_0 = 0.62`, :math:`a_1 = 0.48` and :math:`a_2=0.38` .. plot:: :width: 80% :include-source: from spectrum import window_visu window_visu(64, 'bartlett_hann') .. seealso:: :func:`create_window`, :class:`Window` """ if N == 1: return ones(1) n = arange(0, N) a0 = 0.62 a1 = 0.48 a2 = 0.38 win = a0 - a1 *abs(n/(N-1.)-0.5) -a2 * cos(2*pi*n/(N-1.)) return win
def _coeff4(N, a0, a1, a2, a3): """a common internal function to some window functions with 4 coeffs For the blackmna harris for instance, the results are identical to octave if N is odd but not for even values...if n =0 whatever N is, the w(0) must be equal to a0-a1+a2-a3, which is the case here, but not in octave...""" if N == 1: return ones(1) n = arange(0, N) N1 = N - 1. w = a0 -a1*cos(2.*pi*n / N1) + a2*cos(4.*pi*n / N1) - a3*cos(6.*pi*n / N1) return w
[docs]def window_nuttall(N): r"""Nuttall tapering window :param N: window length .. math:: w(n) = a_0 - a_1 \cos\left(\frac{2\pi n}{N-1}\right)+ a_2 \cos\left(\frac{4\pi n}{N-1}\right)- a_3 \cos\left(\frac{6\pi n}{N-1}\right) with :math:`a_0 = 0.355768`, :math:`a_1 = 0.487396`, :math:`a_2=0.144232` and :math:`a_3=0.012604` .. plot:: :width: 80% :include-source: from spectrum import window_visu window_visu(64, 'nuttall', mindB=-80) .. seealso:: :func:`create_window`, :class:`Window` """ a0 = 0.355768 a1 = 0.487396 a2 = 0.144232 a3 = 0.012604 return _coeff4(N, a0, a1, a2, a3)
[docs]def window_blackman_nuttall(N): r"""Blackman Nuttall window returns a minimum, 4-term Blackman-Harris window. The window is minimum in the sense that its maximum sidelobes are minimized. The coefficients for this window differ from the Blackman-Harris window coefficients and produce slightly lower sidelobes. :param N: window length .. math:: w(n) = a_0 - a_1 \cos\left(\frac{2\pi n}{N-1}\right)+ a_2 \cos\left(\frac{4\pi n}{N-1}\right)- a_3 \cos\left(\frac{6\pi n}{N-1}\right) with :math:`a_0 = 0.3635819`, :math:`a_1 = 0.4891775`, :math:`a_2=0.1365995` and :math:`0_3=.0106411` .. plot:: :width: 80% :include-source: from spectrum import window_visu window_visu(64, 'blackman_nuttall', mindB=-80) .. seealso:: :func:`spectrum.window.create_window` .. seealso:: :func:`create_window`, :class:`Window` """ a0 = 0.3635819 a1 = 0.4891775 a2 = 0.1365995 a3 = 0.0106411 return _coeff4(N, a0, a1, a2, a3)
[docs]def window_blackman_harris(N): r"""Blackman Harris window :param N: window length .. math:: w(n) = a_0 - a_1 \cos\left(\frac{2\pi n}{N-1}\right)+ a_2 \cos\left(\frac{4\pi n}{N-1}\right)- a_3 \cos\left(\frac{6\pi n}{N-1}\right) =============== ========= coeff value =============== ========= :math:`a_0` 0.35875 :math:`a_1` 0.48829 :math:`a_2` 0.14128 :math:`a_3` 0.01168 =============== ========= .. plot:: :width: 80% :include-source: from spectrum import window_visu window_visu(64, 'blackman_harris', mindB=-80) .. seealso:: :func:`spectrum.window.create_window` .. seealso:: :func:`create_window`, :class:`Window` """ a0 = 0.35875 a1 = 0.48829 a2 = 0.14128 a3 = 0.01168 return _coeff4(N, a0, a1, a2, a3)
[docs]def window_bohman(N): r"""Bohman tapering window :param N: window length .. math:: w(n) = (1-|x|) \cos (\pi |x|) + \frac{1}{\pi} \sin(\pi |x|) where x is a length N vector of linearly spaced values between -1 and 1. .. plot:: :width: 80% :include-source: from spectrum import window_visu window_visu(64, 'bohman') .. seealso:: :func:`create_window`, :class:`Window` """ x = linspace(-1, 1, N) w = (1.-abs(x)) * cos(pi*abs(x)) + 1./pi * sin(pi*abs(x)) return w
[docs]def window_tukey(N, r=0.5): """Tukey tapering window (or cosine-tapered window) :param N: window length :param r: defines the ratio between the constant section and the cosine section. It has to be between 0 and 1. The function returns a Hanning window for `r=0` and a full box for `r=1`. .. plot:: :width: 80% :include-source: from spectrum import window_visu window_visu(64, 'tukey') window_visu(64, 'tukey', r=1) .. math:: 0.5 (1+cos(2pi/r (x-r/2))) for 0<=x<r/2 .. math:: 0.5 (1+cos(2pi/r (x-1+r/2))) for x>=1-r/2 .. seealso:: :func:`create_window`, :class:`Window` """ assert r>=0 and r<=1 , "r must be in [0,1]" if N==1: return ones(1) if r == 0: return ones(N) elif r == 1: return window_hann(N) else: from numpy import flipud, concatenate, where ## cosine-tapered window x = linspace(0, 1, N) x1 = where(x<r/2.) w = 0.5*(1+cos(2*pi/r*(x[x1[0]]-r/2))) w = concatenate((w, ones(N-len(w)*2), flipud(w))) return w
[docs]def window_parzen(N): r"""Parsen tapering window (also known as de la Valle-Poussin) :param N: window length Parzen windows are piecewise cubic approximations of Gaussian windows. Parzen window sidelobes fall off as :math:`1/\omega^4`. if :math:`0\leq|x|\leq (N-1)/4`: .. math:: w(n) = 1-6 \left( \frac{|n|}{N/2} \right)^2 +6 \left( \frac{|n|}{N/2}\right)^3 if :math:`(N-1)/4\leq|x|\leq (N-1)/2` .. math:: w(n) = 2 \left(1- \frac{|n|}{N/2}\right)^3 .. plot:: :width: 80% :include-source: from spectrum import window_visu window_visu(64, 'parzen') .. seealso:: :func:`create_window`, :class:`Window` """ from numpy import where, concatenate n = linspace(-(N-1)/2., (N-1)/2., N) n1 = n[where(abs(n)<=(N-1)/4.)[0]] n2 = n[where(n>(N-1)/4.)[0]] n3 = n[where(n<-(N-1)/4.)[0]] w1 = 1. -6.*(abs(n1)/(N/2.))**2 + 6*(abs(n1)/(N/2.))**3 w2 = 2.*(1-abs(n2)/(N/2.))**3 w3 = 2.*(1-abs(n3)/(N/2.))**3 w = concatenate((w3, w1, w2)) return w
[docs]def window_flattop(N, mode='symmetric',precision=None): r"""Flat-top tapering window Returns symmetric or periodic flat top window. :param N: window length :param mode: way the data are normalised. If mode is *symmetric*, then divide n by N-1. IF mode is *periodic*, divide by N, to be consistent with octave code. When using windows for filter design, the *symmetric* mode should be used (default). When using windows for spectral analysis, the *periodic* mode should be used. The mathematical form of the flat-top window in the symmetric case is: .. math:: w(n) = a_0 - a_1 \cos\left(\frac{2\pi n}{N-1}\right) + a_2 \cos\left(\frac{4\pi n}{N-1}\right) - a_3 \cos\left(\frac{6\pi n}{N-1}\right) + a_4 \cos\left(\frac{8\pi n}{N-1}\right) ===== ============= coeff value ===== ============= a0 0.21557895 a1 0.41663158 a2 0.277263158 a3 0.083578947 a4 0.006947368 ===== ============= .. plot:: :width: 80% :include-source: from spectrum import window_visu window_visu(64, 'bohman') .. seealso:: :func:`create_window`, :class:`Window` """ assert mode in ['periodic', 'symmetric'] t = arange(0, N) # FIXME: N=1 for mode = periodic ? if mode == 'periodic': x = 2*pi*t/float(N) else: if N ==1: return ones(1) x = 2*pi*t/float(N-1) a0 = 0.21557895 a1 = 0.41663158 a2 = 0.277263158 a3 = 0.083578947 a4 = 0.006947368 if precision == 'octave': #to compare with octave, same as above but less precise d = 4.6402 a0 = 1./d a1 = 1.93/d a2 = 1.29/d a3 = 0.388/d a4 = 0.0322/d w = a0-a1*cos(x)+a2*cos(2*x)-a3*cos(3*x)+a4*cos(4*x) return w
[docs]def window_taylor(N, nbar=4, sll=-30): """Taylor tapering window Taylor windows allows you to make tradeoffs between the mainlobe width and sidelobe level (sll). Implemented as described by Carrara, Goodman, and Majewski in 'Spotlight Synthetic Aperture Radar: Signal Processing Algorithms' Pages 512-513 :param N: window length :param float nbar: :param float sll: The default values gives equal height sidelobes (nbar) and maximum sidelobe level (sll). .. warning:: not implemented .. seealso:: :func:`create_window`, :class:`Window` """ B = 10**(-sll/20) A = log(B + sqrt(B**2 - 1))/pi s2 = nbar**2 / (A**2 + (nbar - 0.5)**2) ma = arange(1,nbar) def calc_Fm(m): numer = (-1)**(m+1) * prod(1-m**2/s2/(A**2 + (ma - 0.5)**2)) denom = 2* prod([ 1-m**2/j**2 for j in ma if j != m]) return numer/denom Fm = array([calc_Fm(m) for m in ma]) def W(n): return 2 * np.sum(Fm * cos(2*pi*ma*(n-N/2 + 1/2)/N)) + 1 w = array([W(n) for n in range(N)]) # normalize (Note that this is not described in the original text) scale = W((N-1)/2) w /= scale return w
[docs]def window_riesz(N): r"""Riesz tapering window :param N: window length .. math:: w(n) = 1 - \left| \frac{n}{N/2} \right|^2 with :math:`-N/2 \leq n \leq N/2`. .. plot:: :width: 80% :include-source: from spectrum import window_visu window_visu(64, 'riesz') .. seealso:: :func:`create_window`, :class:`Window` """ n = linspace(-N/2., (N)/2., N) w = 1 - abs(n/(N/2.))**2. return w
[docs]def window_riemann(N): r"""Riemann tapering window :param int N: window length .. math:: w(n) = 1 - \left| \frac{n}{N/2} \right|^2 with :math:`-N/2 \leq n \leq N/2`. .. plot:: :width: 80% :include-source: from spectrum import window_visu window_visu(64, 'riesz') .. seealso:: :func:`create_window`, :class:`Window` """ n = linspace(-N/2., (N)/2., N) w = sin(n/float(N)*2.*pi) / (n / float(N)*2.*pi) return w
[docs]def window_poisson(N, alpha=2): r"""Poisson tapering window :param int N: window length .. math:: w(n) = \exp^{-\alpha \frac{|n|}{N/2} } with :math:`-N/2 \leq n \leq N/2`. .. plot:: :width: 80% :include-source: from spectrum import window_visu window_visu(64, 'poisson') window_visu(64, 'poisson', alpha=3) window_visu(64, 'poisson', alpha=4) .. seealso:: :func:`create_window`, :class:`Window` """ n = linspace(-N/2., (N)/2., N) w = exp(-alpha * abs(n)/(N/2.)) return w
[docs]def window_poisson_hanning(N, alpha=2): r"""Hann-Poisson tapering window This window is constructed as the product of the Hanning and Poisson windows. The parameter **alpha** is the Poisson parameter. :param int N: window length :param float alpha: parameter of the poisson window .. plot:: :width: 80% :include-source: from spectrum import window_visu window_visu(64, 'poisson_hanning', alpha=0.5) window_visu(64, 'poisson_hanning', alpha=1) window_visu(64, 'poisson_hanning') .. seealso:: :func:`window_poisson`, :func:`window_hann` """ w1 = window_hann(N) w2 = window_poisson(N, alpha=alpha) return w1*w2
[docs]def window_cauchy(N, alpha=3): r"""Cauchy tapering window :param int N: window length :param float alpha: parameter of the poisson window .. math:: w(n) = \frac{1}{1+\left(\frac{\alpha*n}{N/2}\right)**2} .. plot:: :width: 80% :include-source: from spectrum import window_visu window_visu(64, 'cauchy', alpha=3) window_visu(64, 'cauchy', alpha=4) window_visu(64, 'cauchy', alpha=5) .. seealso:: :func:`window_poisson`, :func:`window_hann` """ n = linspace(-N/2., (N)/2., N) w = 1./(1.+ (alpha*n/(N/2.))**2) return w