From ae633e87e91660c1feba5879076172ddfaced7f0 Mon Sep 17 00:00:00 2001 From: nyanpasu64 Date: Sat, 22 Dec 2018 01:30:44 -0800 Subject: [PATCH] Extract Scipy windows module --- corrscope/utils/scipy_windows.py | 2107 ++++++++++++++++++++++++++++++ 1 file changed, 2107 insertions(+) create mode 100644 corrscope/utils/scipy_windows.py diff --git a/corrscope/utils/scipy_windows.py b/corrscope/utils/scipy_windows.py new file mode 100644 index 0000000..7a15c30 --- /dev/null +++ b/corrscope/utils/scipy_windows.py @@ -0,0 +1,2107 @@ +"""The suite of window functions.""" + +from __future__ import division, print_function, absolute_import + +import operator +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_cosine','general_gaussian', + 'general_hamming', 'chebwin', 'slepian', 'cosine', 'hann', + 'exponential', 'tukey', 'dpss', 'get_window'] + + +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 + + +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 general_cosine(M, a, sym=True): + r""" + Generic weighted sum of cosine terms window + + Parameters + ---------- + M : int + Number of points in the output window + a : array_like + Sequence of weighting coefficients. This uses the convention of being + centered on the origin, so these will typically all be positive + numbers, not alternating sign. + 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. + + References + ---------- + .. [1] A. Nuttall, "Some windows with very good sidelobe behavior," IEEE + Transactions on Acoustics, Speech, and Signal Processing, vol. 29, + no. 1, pp. 84-91, Feb 1981. :doi:`10.1109/TASSP.1981.1163506`. + .. [2] Heinzel G. et al., "Spectrum and spectral density estimation by the + Discrete Fourier transform (DFT), including a comprehensive list of + window functions and some new flat-top windows", February 15, 2002 + https://holometer.fnal.gov/GH_FFT.pdf + + Examples + -------- + Heinzel describes a flat-top window named "HFT90D" with formula: [2]_ + + .. math:: w_j = 1 - 1.942604 \cos(z) + 1.340318 \cos(2z) + - 0.440811 \cos(3z) + 0.043097 \cos(4z) + + where + + .. math:: z = \frac{2 \pi j}{N}, j = 0...N - 1 + + Since this uses the convention of starting at the origin, to reproduce the + window, we need to convert every other coefficient to a positive number: + + >>> HFT90D = [1, 1.942604, 1.340318, 0.440811, 0.043097] + + The paper states that the highest sidelobe is at -90.2 dB. Reproduce + Figure 42 by plotting the window and its frequency response, and confirm + the sidelobe level in red: + + >>> from scipy.signal.windows import general_cosine + >>> from scipy.fftpack import fft, fftshift + >>> import matplotlib.pyplot as plt + + >>> window = general_cosine(1000, HFT90D, sym=False) + >>> plt.plot(window) + >>> plt.title("HFT90D window") + >>> plt.ylabel("Amplitude") + >>> plt.xlabel("Sample") + + >>> plt.figure() + >>> A = fft(window, 10000) / (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([-50/1000, 50/1000, -140, 0]) + >>> plt.title("Frequency response of the HFT90D window") + >>> plt.ylabel("Normalized magnitude [dB]") + >>> plt.xlabel("Normalized frequency [cycles per sample]") + >>> plt.axhline(-90.2, color='red') + >>> plt.show() + """ + if _len_guards(M): + return np.ones(M) + M, needs_trunc = _extend(M, sym) + + fac = np.linspace(-np.pi, np.pi, M) + w = np.zeros(M) + for k in range(len(a)): + w += a[k] * np.cos(k * fac) + + return _truncate(w, needs_trunc) + + +def boxcar(M, sym=True): + """Return a boxcar or rectangular window. + + Also known as a rectangular window or Dirichlet window, this is equivalent + to no window at all. + + Parameters + ---------- + M : int + Number of points in the output window. If zero or less, an empty + array is returned. + sym : bool, optional + Whether the window is symmetric. (Has no effect for boxcar.) + + Returns + ------- + w : ndarray + The window, with the maximum value normalized to 1. + + 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.boxcar(51) + >>> plt.plot(window) + >>> plt.title("Boxcar 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 boxcar window") + >>> plt.ylabel("Normalized magnitude [dB]") + >>> plt.xlabel("Normalized frequency [cycles per sample]") + + """ + if _len_guards(M): + return np.ones(M) + M, needs_trunc = _extend(M, sym) + + w = np.ones(M, float) + + return _truncate(w, needs_trunc) + + +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). + + See Also + -------- + bartlett : A triangular window that touches zero + + 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 _len_guards(M): + return np.ones(M) + M, needs_trunc = _extend(M, sym) + + 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 _truncate(w, needs_trunc) + + +def parzen(M, sym=True): + """Return a Parzen 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). + + References + ---------- + .. [1] E. Parzen, "Mathematical Considerations in the Estimation of + Spectra", Technometrics, Vol. 3, No. 2 (May, 1961), pp. 167-190 + + 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.parzen(51) + >>> plt.plot(window) + >>> plt.title("Parzen 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 Parzen window") + >>> plt.ylabel("Normalized magnitude [dB]") + >>> plt.xlabel("Normalized frequency [cycles per sample]") + + """ + if _len_guards(M): + return np.ones(M) + M, needs_trunc = _extend(M, sym) + + 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 _truncate(w, needs_trunc) + + +def bohman(M, sym=True): + """Return a Bohman 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.bohman(51) + >>> plt.plot(window) + >>> plt.title("Bohman 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 Bohman window") + >>> plt.ylabel("Normalized magnitude [dB]") + >>> plt.xlabel("Normalized frequency [cycles per sample]") + + """ + if _len_guards(M): + return np.ones(M) + M, needs_trunc = _extend(M, sym) + + fac = np.abs(np.linspace(-1, 1, M)[1:-1]) + w = (1 - fac) * np.cos(np.pi * fac) + 1.0 / np.pi * np.sin(np.pi * fac) + w = np.r_[0, w, 0] + + return _truncate(w, needs_trunc) + + +def blackman(M, sym=True): + r""" + Return a Blackman window. + + The Blackman window is a taper formed by using the first three terms of + a summation of cosines. It was designed to have close to the minimal + leakage possible. It is close to optimal, only slightly worse than a + Kaiser 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). + + Notes + ----- + The Blackman window is defined as + + .. math:: w(n) = 0.42 - 0.5 \cos(2\pi n/M) + 0.08 \cos(4\pi n/M) + + The "exact Blackman" window was designed to null out the third and fourth + sidelobes, but has discontinuities at the boundaries, resulting in a + 6 dB/oct fall-off. This window is an approximation of the "exact" window, + which does not null the sidelobes as well, but is smooth at the edges, + improving the fall-off rate to 18 dB/oct. [3]_ + + Most references to the Blackman window come from the signal processing + literature, where it is used as one of many windowing functions for + smoothing values. It is also known as an apodization (which means + "removing the foot", i.e. smoothing discontinuities at the beginning + and end of the sampled signal) or tapering function. It is known as a + "near optimal" tapering function, almost as good (by some measures) + as the Kaiser window. + + References + ---------- + .. [1] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power + spectra, Dover Publications, New York. + .. [2] Oppenheim, A.V., and R.W. Schafer. Discrete-Time Signal Processing. + Upper Saddle River, NJ: Prentice-Hall, 1999, pp. 468-471. + .. [3] Harris, Fredric J. (Jan 1978). "On the use of Windows for Harmonic + Analysis with the Discrete Fourier Transform". Proceedings of the + IEEE 66 (1): 51-83. :doi:`10.1109/PROC.1978.10837`. + + 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.blackman(51) + >>> plt.plot(window) + >>> plt.title("Blackman 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 Blackman window") + >>> plt.ylabel("Normalized magnitude [dB]") + >>> plt.xlabel("Normalized frequency [cycles per sample]") + + """ + # Docstring adapted from NumPy's blackman function + return general_cosine(M, [0.42, 0.50, 0.08], sym) + + +def nuttall(M, sym=True): + """Return a minimum 4-term Blackman-Harris window according to Nuttall. + + This variation is called "Nuttall4c" by Heinzel. [2]_ + + 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). + + References + ---------- + .. [1] A. Nuttall, "Some windows with very good sidelobe behavior," IEEE + Transactions on Acoustics, Speech, and Signal Processing, vol. 29, + no. 1, pp. 84-91, Feb 1981. :doi:`10.1109/TASSP.1981.1163506`. + .. [2] Heinzel G. et al., "Spectrum and spectral density estimation by the + Discrete Fourier transform (DFT), including a comprehensive list of + window functions and some new flat-top windows", February 15, 2002 + https://holometer.fnal.gov/GH_FFT.pdf + + 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.nuttall(51) + >>> plt.plot(window) + >>> plt.title("Nuttall 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 Nuttall window") + >>> plt.ylabel("Normalized magnitude [dB]") + >>> plt.xlabel("Normalized frequency [cycles per sample]") + + """ + return general_cosine(M, [0.3635819, 0.4891775, 0.1365995, 0.0106411], sym) + + +def blackmanharris(M, sym=True): + """Return a minimum 4-term Blackman-Harris 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.blackmanharris(51) + >>> plt.plot(window) + >>> plt.title("Blackman-Harris 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 Blackman-Harris window") + >>> plt.ylabel("Normalized magnitude [dB]") + >>> plt.xlabel("Normalized frequency [cycles per sample]") + + """ + return general_cosine(M, [0.35875, 0.48829, 0.14128, 0.01168], sym) + + +def flattop(M, sym=True): + """Return a flat top 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). + + Notes + ----- + Flat top windows are used for taking accurate measurements of signal + amplitude in the frequency domain, with minimal scalloping error from the + center of a frequency bin to its edges, compared to others. This is a + 5th-order cosine window, with the 5 terms optimized to make the main lobe + maximally flat. [1]_ + + References + ---------- + .. [1] D'Antona, Gabriele, and A. Ferrero, "Digital Signal Processing for + Measurement Systems", Springer Media, 2006, p. 70 + :doi:`10.1007/0-387-28666-7`. + + 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.flattop(51) + >>> plt.plot(window) + >>> plt.title("Flat top 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 flat top window") + >>> plt.ylabel("Normalized magnitude [dB]") + >>> plt.xlabel("Normalized frequency [cycles per sample]") + + """ + a = [0.21557895, 0.41663158, 0.277263158, 0.083578947, 0.006947368] + return general_cosine(M, a, sym) + + +def bartlett(M, sym=True): + r""" + Return a Bartlett window. + + The Bartlett window is very similar to a triangular window, except + that the end points are at zero. It is often used in signal + processing for tapering a signal, without generating too much + ripple in the frequency domain. + + 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 triangular window, with the first and last samples equal to zero + and the maximum value normalized to 1 (though the value 1 does not + appear if `M` is even and `sym` is True). + + See Also + -------- + triang : A triangular window that does not touch zero at the ends + + Notes + ----- + The Bartlett window is defined as + + .. math:: w(n) = \frac{2}{M-1} \left( + \frac{M-1}{2} - \left|n - \frac{M-1}{2}\right| + \right) + + Most references to the Bartlett window come from the signal + processing literature, where it is used as one of many windowing + functions for smoothing values. Note that convolution with this + window produces linear interpolation. It is also known as an + apodization (which means"removing the foot", i.e. smoothing + discontinuities at the beginning and end of the sampled signal) or + tapering function. The Fourier transform of the Bartlett is the product + of two sinc functions. + Note the excellent discussion in Kanasewich. [2]_ + + References + ---------- + .. [1] M.S. Bartlett, "Periodogram Analysis and Continuous Spectra", + Biometrika 37, 1-16, 1950. + .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", + The University of Alberta Press, 1975, pp. 109-110. + .. [3] A.V. Oppenheim and R.W. Schafer, "Discrete-Time Signal + Processing", Prentice-Hall, 1999, pp. 468-471. + .. [4] Wikipedia, "Window function", + https://en.wikipedia.org/wiki/Window_function + .. [5] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling, + "Numerical Recipes", Cambridge University Press, 1986, page 429. + + 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.bartlett(51) + >>> plt.plot(window) + >>> plt.title("Bartlett 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 Bartlett window") + >>> plt.ylabel("Normalized magnitude [dB]") + >>> plt.xlabel("Normalized frequency [cycles per sample]") + + """ + # Docstring adapted from NumPy's bartlett function + if _len_guards(M): + return np.ones(M) + M, needs_trunc = _extend(M, sym) + + n = np.arange(0, M) + w = np.where(np.less_equal(n, (M - 1) / 2.0), + 2.0 * n / (M - 1), 2.0 - 2.0 * n / (M - 1)) + + return _truncate(w, needs_trunc) + + +def hann(M, sym=True): + r""" + Return a Hann window. + + The Hann window is a taper formed by using a raised cosine or sine-squared + with ends that touch zero. + + 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). + + Notes + ----- + The Hann window is defined as + + .. math:: w(n) = 0.5 - 0.5 \cos\left(\frac{2\pi{n}}{M-1}\right) + \qquad 0 \leq n \leq M-1 + + The window was named for Julius von Hann, an Austrian meteorologist. It is + also known as the Cosine Bell. It is sometimes erroneously referred to as + the "Hanning" window, from the use of "hann" as a verb in the original + paper and confusion with the very similar Hamming window. + + Most references to the Hann window come from the signal processing + literature, where it is used as one of many windowing functions for + smoothing values. It is also known as an apodization (which means + "removing the foot", i.e. smoothing discontinuities at the beginning + and end of the sampled signal) or tapering function. + + References + ---------- + .. [1] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power + spectra, Dover Publications, New York. + .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", + The University of Alberta Press, 1975, pp. 106-108. + .. [3] Wikipedia, "Window function", + https://en.wikipedia.org/wiki/Window_function + .. [4] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling, + "Numerical Recipes", Cambridge University Press, 1986, page 425. + + 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.hann(51) + >>> plt.plot(window) + >>> plt.title("Hann 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 Hann window") + >>> plt.ylabel("Normalized magnitude [dB]") + >>> plt.xlabel("Normalized frequency [cycles per sample]") + + """ + # Docstring adapted from NumPy's hanning function + return general_hamming(M, 0.5, sym) + + +@np.deprecate(new_name='scipy.signal.windows.hann') +def hanning(*args, **kwargs): + return hann(*args, **kwargs) + + +def tukey(M, alpha=0.5, sym=True): + r"""Return a Tukey window, also known as a tapered cosine window. + + Parameters + ---------- + M : int + Number of points in the output window. If zero or less, an empty + array is returned. + alpha : float, optional + Shape parameter of the Tukey window, representing the fraction of the + window inside the cosine tapered region. + If zero, the Tukey window is equivalent to a rectangular window. + If one, the Tukey window is equivalent to a Hann window. + 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). + + References + ---------- + .. [1] Harris, Fredric J. (Jan 1978). "On the use of Windows for Harmonic + Analysis with the Discrete Fourier Transform". Proceedings of the + IEEE 66 (1): 51-83. :doi:`10.1109/PROC.1978.10837` + .. [2] Wikipedia, "Window function", + https://en.wikipedia.org/wiki/Window_function#Tukey_window + + 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.tukey(51) + >>> plt.plot(window) + >>> plt.title("Tukey window") + >>> plt.ylabel("Amplitude") + >>> plt.xlabel("Sample") + >>> plt.ylim([0, 1.1]) + + >>> 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 Tukey window") + >>> plt.ylabel("Normalized magnitude [dB]") + >>> plt.xlabel("Normalized frequency [cycles per sample]") + + """ + if _len_guards(M): + return np.ones(M) + + if alpha <= 0: + return np.ones(M, 'd') + elif alpha >= 1.0: + return hann(M, sym=sym) + + M, needs_trunc = _extend(M, sym) + + 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 _truncate(w, needs_trunc) + + +def barthann(M, sym=True): + """Return a modified Bartlett-Hann 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.barthann(51) + >>> plt.plot(window) + >>> plt.title("Bartlett-Hann 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 Bartlett-Hann window") + >>> plt.ylabel("Normalized magnitude [dB]") + >>> plt.xlabel("Normalized frequency [cycles per sample]") + + """ + if _len_guards(M): + return np.ones(M) + M, needs_trunc = _extend(M, sym) + + n = np.arange(0, M) + fac = np.abs(n / (M - 1.0) - 0.5) + w = 0.62 - 0.48 * fac + 0.38 * np.cos(2 * np.pi * fac) + + return _truncate(w, needs_trunc) + + +def general_hamming(M, alpha, sym=True): + r"""Return a generalized Hamming window. + + The generalized Hamming window is constructed by multiplying a rectangular + window by one period of a cosine function [1]_. + + Parameters + ---------- + M : int + Number of points in the output window. If zero or less, an empty + array is returned. + alpha : float + The window coefficient, :math:`\alpha` + 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 generalized Hamming window is defined as + + .. math:: w(n) = \alpha - \left(1 - \alpha\right) \cos\left(\frac{2\pi{n}}{M-1}\right) + \qquad 0 \leq n \leq M-1 + + Both the common Hamming window and Hann window are special cases of the + generalized Hamming window with :math:`\alpha` = 0.54 and :math:`\alpha` = + 0.5, respectively [2]_. + + See Also + -------- + hamming, hann + + Examples + -------- + The Sentinel-1A/B Instrument Processing Facility uses generalized Hamming + windows in the processing of spaceborne Synthetic Aperture Radar (SAR) + data [3]_. The facility uses various values for the :math:`\alpha` + parameter based on operating mode of the SAR instrument. Some common + :math:`\alpha` values include 0.75, 0.7 and 0.52 [4]_. As an example, we + plot these different windows. + + >>> from scipy.signal.windows import general_hamming + >>> from scipy.fftpack import fft, fftshift + >>> import matplotlib.pyplot as plt + + >>> fig1, spatial_plot = plt.subplots() + >>> spatial_plot.set_title("Generalized Hamming Windows") + >>> spatial_plot.set_ylabel("Amplitude") + >>> spatial_plot.set_xlabel("Sample") + + >>> fig2, freq_plot = plt.subplots() + >>> freq_plot.set_title("Frequency Responses") + >>> freq_plot.set_ylabel("Normalized magnitude [dB]") + >>> freq_plot.set_xlabel("Normalized frequency [cycles per sample]") + + >>> for alpha in [0.75, 0.7, 0.52]: + ... window = general_hamming(41, alpha) + ... spatial_plot.plot(window, label="{:.2f}".format(alpha)) + ... 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()))) + ... freq_plot.plot(freq, response, label="{:.2f}".format(alpha)) + >>> freq_plot.legend(loc="upper right") + >>> spatial_plot.legend(loc="upper right") + + References + ---------- + .. [1] DSPRelated, "Generalized Hamming Window Family", + https://www.dsprelated.com/freebooks/sasp/Generalized_Hamming_Window_Family.html + .. [2] Wikipedia, "Window function", + https://en.wikipedia.org/wiki/Window_function + .. [3] Riccardo Piantanida ESA, "Sentinel-1 Level 1 Detailed Algorithm + Definition", + https://sentinel.esa.int/documents/247904/1877131/Sentinel-1-Level-1-Detailed-Algorithm-Definition + .. [4] Matthieu Bourbigot ESA, "Sentinel-1 Product Definition", + https://sentinel.esa.int/documents/247904/1877131/Sentinel-1-Product-Definition + """ + return general_cosine(M, [alpha, 1. - alpha], sym) + + +def hamming(M, sym=True): + r"""Return a Hamming window. + + The Hamming window is a taper formed by using a raised cosine with + non-zero endpoints, optimized to minimize the nearest side lobe. + + 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). + + Notes + ----- + The Hamming window is defined as + + .. math:: w(n) = 0.54 - 0.46 \cos\left(\frac{2\pi{n}}{M-1}\right) + \qquad 0 \leq n \leq M-1 + + The Hamming was named for R. W. Hamming, an associate of J. W. Tukey and + is described in Blackman and Tukey. It was recommended for smoothing the + truncated autocovariance function in the time domain. + Most references to the Hamming window come from the signal processing + literature, where it is used as one of many windowing functions for + smoothing values. It is also known as an apodization (which means + "removing the foot", i.e. smoothing discontinuities at the beginning + and end of the sampled signal) or tapering function. + + References + ---------- + .. [1] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power + spectra, Dover Publications, New York. + .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", The + University of Alberta Press, 1975, pp. 109-110. + .. [3] Wikipedia, "Window function", + https://en.wikipedia.org/wiki/Window_function + .. [4] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling, + "Numerical Recipes", Cambridge University Press, 1986, page 425. + + 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.hamming(51) + >>> plt.plot(window) + >>> plt.title("Hamming 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 Hamming window") + >>> plt.ylabel("Normalized magnitude [dB]") + >>> plt.xlabel("Normalized frequency [cycles per sample]") + + """ + # Docstring adapted from NumPy's hamming function + return general_hamming(M, 0.54, sym) + + +def kaiser(M, beta, sym=True): + r"""Return a Kaiser window. + + The Kaiser window is a taper formed by using a Bessel function. + + Parameters + ---------- + M : int + Number of points in the output window. If zero or less, an empty + array is returned. + beta : float + Shape parameter, determines trade-off between main-lobe width and + side lobe level. As beta gets large, the window narrows. + 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 Kaiser window is defined as + + .. math:: w(n) = I_0\left( \beta \sqrt{1-\frac{4n^2}{(M-1)^2}} + \right)/I_0(\beta) + + with + + .. math:: \quad -\frac{M-1}{2} \leq n \leq \frac{M-1}{2}, + + where :math:`I_0` is the modified zeroth-order Bessel function. + + The Kaiser was named for Jim Kaiser, who discovered a simple approximation + to the DPSS window based on Bessel functions. + The Kaiser window is a very good approximation to the Digital Prolate + Spheroidal Sequence, or Slepian window, which is the transform which + maximizes the energy in the main lobe of the window relative to total + energy. + + The Kaiser can approximate other windows by varying the beta parameter. + (Some literature uses alpha = beta/pi.) [4]_ + + ==== ======================= + beta Window shape + ==== ======================= + 0 Rectangular + 5 Similar to a Hamming + 6 Similar to a Hann + 8.6 Similar to a Blackman + ==== ======================= + + A beta value of 14 is probably a good starting point. Note that as beta + gets large, the window narrows, and so the number of samples needs to be + large enough to sample the increasingly narrow spike, otherwise NaNs will + be returned. + + Most references to the Kaiser window come from the signal processing + literature, where it is used as one of many windowing functions for + smoothing values. It is also known as an apodization (which means + "removing the foot", i.e. smoothing discontinuities at the beginning + and end of the sampled signal) or tapering function. + + References + ---------- + .. [1] J. F. Kaiser, "Digital Filters" - Ch 7 in "Systems analysis by + digital computer", Editors: F.F. Kuo and J.F. Kaiser, p 218-285. + John Wiley and Sons, New York, (1966). + .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", The + University of Alberta Press, 1975, pp. 177-178. + .. [3] Wikipedia, "Window function", + https://en.wikipedia.org/wiki/Window_function + .. [4] F. J. Harris, "On the use of windows for harmonic analysis with the + discrete Fourier transform," Proceedings of the IEEE, vol. 66, + no. 1, pp. 51-83, Jan. 1978. :doi:`10.1109/PROC.1978.10837`. + + 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.kaiser(51, beta=14) + >>> plt.plot(window) + >>> plt.title(r"Kaiser window ($\beta$=14)") + >>> 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(r"Frequency response of the Kaiser window ($\beta$=14)") + >>> plt.ylabel("Normalized magnitude [dB]") + >>> plt.xlabel("Normalized frequency [cycles per sample]") + + """ + # Docstring adapted from NumPy's kaiser function + if _len_guards(M): + return np.ones(M) + M, needs_trunc = _extend(M, sym) + + n = np.arange(0, M) + alpha = (M - 1) / 2.0 + w = (special.i0(beta * np.sqrt(1 - ((n - alpha) / alpha) ** 2.0)) / + special.i0(beta)) + + return _truncate(w, needs_trunc) + + +def gaussian(M, std, sym=True): + r"""Return a Gaussian window. + + Parameters + ---------- + M : int + Number of points in the output window. If zero or less, an empty + array is returned. + std : float + The standard deviation, sigma. + 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 Gaussian window is defined as + + .. math:: w(n) = e^{ -\frac{1}{2}\left(\frac{n}{\sigma}\right)^2 } + + 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.gaussian(51, std=7) + >>> plt.plot(window) + >>> plt.title(r"Gaussian window ($\sigma$=7)") + >>> 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(r"Frequency response of the Gaussian window ($\sigma$=7)") + >>> plt.ylabel("Normalized magnitude [dB]") + >>> plt.xlabel("Normalized frequency [cycles per sample]") + + """ + if _len_guards(M): + return np.ones(M) + M, needs_trunc = _extend(M, sym) + + n = np.arange(0, M) - (M - 1.0) / 2.0 + sig2 = 2 * std * std + w = np.exp(-n ** 2 / sig2) + + return _truncate(w, needs_trunc) + + +def general_gaussian(M, p, sig, sym=True): + r"""Return a window with a generalized Gaussian shape. + + Parameters + ---------- + M : int + Number of points in the output window. If zero or less, an empty + array is returned. + p : float + Shape parameter. p = 1 is identical to `gaussian`, p = 0.5 is + the same shape as the Laplace distribution. + sig : float + The standard deviation, sigma. + 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 generalized Gaussian window is defined as + + .. math:: w(n) = e^{ -\frac{1}{2}\left|\frac{n}{\sigma}\right|^{2p} } + + the half-power point is at + + .. math:: (2 \log(2))^{1/(2 p)} \sigma + + 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.general_gaussian(51, p=1.5, sig=7) + >>> plt.plot(window) + >>> plt.title(r"Generalized Gaussian window (p=1.5, $\sigma$=7)") + >>> 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(r"Freq. resp. of the gen. Gaussian " + ... "window (p=1.5, $\sigma$=7)") + >>> plt.ylabel("Normalized magnitude [dB]") + >>> plt.xlabel("Normalized frequency [cycles per sample]") + + """ + if _len_guards(M): + return np.ones(M) + M, needs_trunc = _extend(M, sym) + + n = np.arange(0, M) - (M - 1.0) / 2.0 + w = np.exp(-0.5 * np.abs(n / sig) ** (2 * p)) + + return _truncate(w, needs_trunc) + + +# `chebwin` contributed by Kumar Appaiah. +def chebwin(M, at, sym=True): + r"""Return a Dolph-Chebyshev window. + + Parameters + ---------- + M : int + Number of points in the output window. If zero or less, an empty + array is returned. + at : float + Attenuation (in dB). + 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 always normalized to 1 + + Notes + ----- + This window optimizes for the narrowest main lobe width for a given order + `M` and sidelobe equiripple attenuation `at`, using Chebyshev + polynomials. It was originally developed by Dolph to optimize the + directionality of radio antenna arrays. + + Unlike most windows, the Dolph-Chebyshev is defined in terms of its + frequency response: + + .. math:: W(k) = \frac + {\cos\{M \cos^{-1}[\beta \cos(\frac{\pi k}{M})]\}} + {\cosh[M \cosh^{-1}(\beta)]} + + where + + .. math:: \beta = \cosh \left [\frac{1}{M} + \cosh^{-1}(10^\frac{A}{20}) \right ] + + and 0 <= abs(k) <= M-1. A is the attenuation in decibels (`at`). + + The time domain window is then generated using the IFFT, so + power-of-two `M` are the fastest to generate, and prime number `M` are + the slowest. + + The equiripple condition in the frequency domain creates impulses in the + time domain, which appear at the ends of the window. + + References + ---------- + .. [1] C. Dolph, "A current distribution for broadside arrays which + optimizes the relationship between beam width and side-lobe level", + Proceedings of the IEEE, Vol. 34, Issue 6 + .. [2] Peter Lynch, "The Dolph-Chebyshev Window: A Simple Optimal Filter", + American Meteorological Society (April 1997) + http://mathsci.ucd.ie/~plynch/Publications/Dolph.pdf + .. [3] F. J. Harris, "On the use of windows for harmonic analysis with the + discrete Fourier transforms", Proceedings of the IEEE, Vol. 66, + No. 1, January 1978 + + 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.chebwin(51, at=100) + >>> plt.plot(window) + >>> plt.title("Dolph-Chebyshev window (100 dB)") + >>> 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 Dolph-Chebyshev window (100 dB)") + >>> plt.ylabel("Normalized magnitude [dB]") + >>> plt.xlabel("Normalized frequency [cycles per sample]") + + """ + if np.abs(at) < 45: + warnings.warn("This window is not suitable for spectral analysis " + "for attenuation values lower than about 45dB because " + "the equivalent noise bandwidth of a Chebyshev window " + "does not grow monotonically with increasing sidelobe " + "attenuation when the attenuation is smaller than " + "about 45 dB.") + if _len_guards(M): + return np.ones(M) + M, needs_trunc = _extend(M, sym) + + # compute the parameter beta + order = M - 1.0 + beta = np.cosh(1.0 / order * np.arccosh(10 ** (np.abs(at) / 20.))) + k = np.r_[0:M] * 1.0 + x = beta * np.cos(np.pi * k / M) + # Find the window's DFT coefficients + # Use analytic definition of Chebyshev polynomial instead of expansion + # from scipy.special. Using the expansion in scipy.special leads to errors. + p = np.zeros(x.shape) + p[x > 1] = np.cosh(order * np.arccosh(x[x > 1])) + p[x < -1] = (2 * (M % 2) - 1) * np.cosh(order * np.arccosh(-x[x < -1])) + p[np.abs(x) <= 1] = np.cos(order * np.arccos(x[np.abs(x) <= 1])) + + # Appropriate IDFT and filling up + # depending on even/odd M + if M % 2: + w = np.real(fftpack.fft(p)) + n = (M + 1) // 2 + w = w[:n] + w = np.concatenate((w[n - 1:0:-1], w)) + else: + p = p * np.exp(1.j * np.pi / M * np.r_[0:M]) + w = np.real(fftpack.fft(p)) + n = M // 2 + 1 + w = np.concatenate((w[n - 1:0:-1], w[1:n])) + w = w / max(w) + + return _truncate(w, needs_trunc) + + +def slepian(M, width, sym=True): + """Return a digital Slepian (DPSS) window. + + Used to maximize the energy concentration in the main lobe. Also called + the digital prolate spheroidal sequence (DPSS). + + .. note:: Deprecated in SciPy 1.1. + `slepian` will be removed in a future version of SciPy, it is + replaced by `dpss`, which uses the standard definition of a + digital Slepian window. + + Parameters + ---------- + M : int + Number of points in the output window. If zero or less, an empty + array is returned. + width : float + Bandwidth + 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 always normalized to 1 + + See Also + -------- + dpss + + References + ---------- + .. [1] D. Slepian & H. O. Pollak: "Prolate spheroidal wave functions, + Fourier analysis and uncertainty-I," Bell Syst. Tech. J., vol.40, + pp.43-63, 1961. https://archive.org/details/bstj40-1-43 + .. [2] H. J. Landau & H. O. Pollak: "Prolate spheroidal wave functions, + Fourier analysis and uncertainty-II," Bell Syst. Tech. J. , vol.40, + pp.65-83, 1961. https://archive.org/details/bstj40-1-65 + + 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.slepian(51, width=0.3) + >>> plt.plot(window) + >>> plt.title("Slepian (DPSS) window (BW=0.3)") + >>> 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 Slepian window (BW=0.3)") + >>> plt.ylabel("Normalized magnitude [dB]") + >>> plt.xlabel("Normalized frequency [cycles per sample]") + + """ + warnings.warn('slepian is deprecated and will be removed in a future ' + 'version, use dpss instead', DeprecationWarning) + if _len_guards(M): + return np.ones(M) + M, needs_trunc = _extend(M, sym) + + # 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 _truncate(win, needs_trunc) + + +def cosine(M, sym=True): + """Return a window with a simple cosine shape. + + 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). + + Notes + ----- + + .. versionadded:: 0.13.0 + + 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.cosine(51) + >>> plt.plot(window) + >>> plt.title("Cosine 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 cosine window") + >>> plt.ylabel("Normalized magnitude [dB]") + >>> plt.xlabel("Normalized frequency [cycles per sample]") + >>> plt.show() + + """ + if _len_guards(M): + return np.ones(M) + M, needs_trunc = _extend(M, sym) + + w = np.sin(np.pi / M * (np.arange(0, M) + .5)) + + return _truncate(w, needs_trunc) + + +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 _len_guards(M): + return np.ones(M) + M, needs_trunc = _extend(M, sym) + + if center is None: + center = (M-1) / 2 + + n = np.arange(0, M) + w = np.exp(-np.abs(n-center) / tau) + + return _truncate(w, needs_trunc) + + +def dpss(M, NW, Kmax=None, sym=True, norm=None, return_ratios=False): + """ + Compute the Discrete Prolate Spheroidal Sequences (DPSS). + + DPSS (or Slepian sequences) are often used in multitaper power spectral + density estimation (see [1]_). The first window in the sequence can be + used to maximize the energy concentration in the main lobe, and is also + called the Slepian window. + + Parameters + ---------- + M : int + Window length. + NW : float + Standardized half bandwidth corresponding to ``2*NW = BW/f0 = BW*N*dt`` + where ``dt`` is taken as 1. + Kmax : int | None, optional + Number of DPSS windows to return (orders ``0`` through ``Kmax-1``). + If None (default), return only a single window of shape ``(M,)`` + instead of an array of windows of shape ``(Kmax, M)``. + 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. + norm : {2, 'approximate', 'subsample'} | None, optional + If 'approximate' or 'subsample', then the windows are normalized by the + maximum, and a correction scale-factor for even-length windows + is applied either using ``M**2/(M**2+NW)`` ("approximate") or + a FFT-based subsample shift ("subsample"), see Notes for details. + If None, then "approximate" is used when ``Kmax=None`` and 2 otherwise + (which uses the l2 norm). + return_ratios : bool, optional + If True, also return the concentration ratios in addition to the + windows. + + Returns + ------- + v : ndarray, shape (Kmax, N) or (N,) + The DPSS windows. Will be 1D if `Kmax` is None. + r : ndarray, shape (Kmax,) or float, optional + The concentration ratios for the windows. Only returned if + `return_ratios` evaluates to True. Will be 0D if `Kmax` is None. + + Notes + ----- + This computation uses the tridiagonal eigenvector formulation given + in [2]_. + + The default normalization for ``Kmax=None``, i.e. window-generation mode, + simply using the l-infinity norm would create a window with two unity + values, which creates slight normalization differences between even and odd + orders. The approximate correction of ``M**2/float(M**2+NW)`` for even + sample numbers is used to counteract this effect (see Examples below). + + For very long signals (e.g., 1e6 elements), it can be useful to compute + windows orders of magnitude shorter and use interpolation (e.g., + `scipy.interpolate.interp1d`) to obtain tapers of length `M`, + but this in general will not preserve orthogonality between the tapers. + + .. versionadded:: 1.1 + + References + ---------- + .. [1] Percival DB, Walden WT. Spectral Analysis for Physical Applications: + Multitaper and Conventional Univariate Techniques. + Cambridge University Press; 1993. + .. [2] Slepian, D. Prolate spheroidal wave functions, Fourier analysis, and + uncertainty V: The discrete case. Bell System Technical Journal, + Volume 57 (1978), 1371430. + .. [3] Kaiser, JF, Schafer RW. On the Use of the I0-Sinh Window for + Spectrum Analysis. IEEE Transactions on Acoustics, Speech and + Signal Processing. ASSP-28 (1): 105-107; 1980. + + Examples + -------- + We can compare the window to `kaiser`, which was invented as an alternative + that was easier to calculate [3]_ (example adapted from + `here `_): + + >>> import numpy as np + >>> import matplotlib.pyplot as plt + >>> from scipy.signal import windows, freqz + >>> N = 51 + >>> fig, axes = plt.subplots(3, 2, figsize=(5, 7)) + >>> for ai, alpha in enumerate((1, 3, 5)): + ... win_dpss = windows.dpss(N, alpha) + ... beta = alpha*np.pi + ... win_kaiser = windows.kaiser(N, beta) + ... for win, c in ((win_dpss, 'k'), (win_kaiser, 'r')): + ... win /= win.sum() + ... axes[ai, 0].plot(win, color=c, lw=1.) + ... axes[ai, 0].set(xlim=[0, N-1], title=r'$\\alpha$ = %s' % alpha, + ... ylabel='Amplitude') + ... w, h = freqz(win) + ... axes[ai, 1].plot(w, 20 * np.log10(np.abs(h)), color=c, lw=1.) + ... axes[ai, 1].set(xlim=[0, np.pi], + ... title=r'$\\beta$ = %0.2f' % beta, + ... ylabel='Magnitude (dB)') + >>> for ax in axes.ravel(): + ... ax.grid(True) + >>> axes[2, 1].legend(['DPSS', 'Kaiser']) + >>> fig.tight_layout() + >>> plt.show() + + And here are examples of the first four windows, along with their + concentration ratios: + + >>> M = 512 + >>> NW = 2.5 + >>> win, eigvals = windows.dpss(M, NW, 4, return_ratios=True) + >>> fig, ax = plt.subplots(1) + >>> ax.plot(win.T, linewidth=1.) + >>> ax.set(xlim=[0, M-1], ylim=[-0.1, 0.1], xlabel='Samples', + ... title='DPSS, M=%d, NW=%0.1f' % (M, NW)) + >>> ax.legend(['win[%d] (%0.4f)' % (ii, ratio) + ... for ii, ratio in enumerate(eigvals)]) + >>> fig.tight_layout() + >>> plt.show() + + Using a standard :math:`l_{\\infty}` norm would produce two unity values + for even `M`, but only one unity value for odd `M`. This produces uneven + window power that can be counteracted by the approximate correction + ``M**2/float(M**2+NW)``, which can be selected by using + ``norm='approximate'`` (which is the same as ``norm=None`` when + ``Kmax=None``, as is the case here). Alternatively, the slower + ``norm='subsample'`` can be used, which uses subsample shifting in the + frequency domain (FFT) to compute the correction: + + >>> Ms = np.arange(1, 41) + >>> factors = (50, 20, 10, 5, 2.0001) + >>> energy = np.empty((3, len(Ms), len(factors))) + >>> for mi, M in enumerate(Ms): + ... for fi, factor in enumerate(factors): + ... NW = M / float(factor) + ... # Corrected using empirical approximation (default) + ... win = windows.dpss(M, NW) + ... energy[0, mi, fi] = np.sum(win ** 2) / np.sqrt(M) + ... # Corrected using subsample shifting + ... win = windows.dpss(M, NW, norm='subsample') + ... energy[1, mi, fi] = np.sum(win ** 2) / np.sqrt(M) + ... # Uncorrected (using l-infinity norm) + ... win /= win.max() + ... energy[2, mi, fi] = np.sum(win ** 2) / np.sqrt(M) + >>> fig, ax = plt.subplots(1) + >>> hs = ax.plot(Ms, energy[2], '-o', markersize=4, + ... markeredgecolor='none') + >>> leg = [hs[-1]] + >>> for hi, hh in enumerate(hs): + ... h1 = ax.plot(Ms, energy[0, :, hi], '-o', markersize=4, + ... color=hh.get_color(), markeredgecolor='none', + ... alpha=0.66) + ... h2 = ax.plot(Ms, energy[1, :, hi], '-o', markersize=4, + ... color=hh.get_color(), markeredgecolor='none', + ... alpha=0.33) + ... if hi == len(hs) - 1: + ... leg.insert(0, h1[0]) + ... leg.insert(0, h2[0]) + >>> ax.set(xlabel='M (samples)', ylabel=r'Power / $\\sqrt{M}$') + >>> ax.legend(leg, ['Uncorrected', r'Corrected: $\\frac{M^2}{M^2+NW}$', + ... 'Corrected (subsample)']) + >>> fig.tight_layout() + + """ # noqa: E501 + if _len_guards(M): + return np.ones(M) + if norm is None: + norm = 'approximate' if Kmax is None else 2 + known_norms = (2, 'approximate', 'subsample') + if norm not in known_norms: + raise ValueError('norm must be one of %s, got %s' + % (known_norms, norm)) + if Kmax is None: + singleton = True + Kmax = 1 + else: + singleton = False + Kmax = operator.index(Kmax) + if not 0 < Kmax <= M: + raise ValueError('Kmax must be greater than 0 and less than M') + if NW >= M/2.: + raise ValueError('NW must be less than M/2.') + if NW <= 0: + raise ValueError('NW must be positive') + M, needs_trunc = _extend(M, sym) + W = float(NW) / M + nidx = np.arange(M) + + # Here we want to set up an optimization problem to find a sequence + # whose energy is maximally concentrated within band [-W,W]. + # Thus, the measure lambda(T,W) is the ratio between the energy within + # that band, and the total energy. This leads to the eigen-system + # (A - (l1)I)v = 0, where the eigenvector corresponding to the largest + # eigenvalue is the sequence with maximally concentrated energy. The + # collection of eigenvectors of this system are called Slepian + # sequences, or discrete prolate spheroidal sequences (DPSS). Only the + # first K, K = 2NW/dt orders of DPSS will exhibit good spectral + # concentration + # [see https://en.wikipedia.org/wiki/Spectral_concentration_problem] + + # Here we set up an alternative symmetric tri-diagonal eigenvalue + # problem such that + # (B - (l2)I)v = 0, and v are our DPSS (but eigenvalues l2 != l1) + # the main diagonal = ([N-1-2*t]/2)**2 cos(2PIW), t=[0,1,2,...,N-1] + # and the first off-diagonal = t(N-t)/2, t=[1,2,...,N-1] + # [see Percival and Walden, 1993] + d = ((M - 1 - 2 * nidx) / 2.) ** 2 * np.cos(2 * np.pi * W) + e = nidx[1:] * (M - nidx[1:]) / 2. + + # only calculate the highest Kmax eigenvalues + w, windows = linalg.eigh_tridiagonal( + d, e, select='i', select_range=(M - Kmax, M - 1)) + w = w[::-1] + windows = windows[:, ::-1].T + + # By convention (Percival and Walden, 1993 pg 379) + # * symmetric tapers (k=0,2,4,...) should have a positive average. + fix_even = (windows[::2].sum(axis=1) < 0) + for i, f in enumerate(fix_even): + if f: + windows[2 * i] *= -1 + # * antisymmetric tapers should begin with a positive lobe + # (this depends on the definition of "lobe", here we'll take the first + # point above the numerical noise, which should be good enough for + # sufficiently smooth functions, and more robust than relying on an + # algorithm that uses max(abs(w)), which is susceptible to numerical + # noise problems) + thresh = max(1e-7, 1. / M) + for i, w in enumerate(windows[1::2]): + if w[w * w > thresh][0] < 0: + windows[2 * i + 1] *= -1 + + # Now find the eigenvalues of the original spectral concentration problem + # Use the autocorr sequence technique from Percival and Walden, 1993 pg 390 + if return_ratios: + dpss_rxx = _fftautocorr(windows) + r = 4 * W * np.sinc(2 * W * nidx) + r[0] = 2 * W + ratios = np.dot(dpss_rxx, r) + if singleton: + ratios = ratios[0] + # Deal with sym and Kmax=None + if norm != 2: + windows /= windows.max() + if M % 2 == 0: + if norm == 'approximate': + correction = M**2 / float(M**2 + NW) + else: + s = np.fft.rfft(windows[0]) + shift = -(1 - 1./M) * np.arange(1, M//2 + 1) + s[1:] *= 2 * np.exp(-1j * np.pi * shift) + correction = M / s.real.sum() + windows *= correction + # else we're already l2 normed, so do nothing + if needs_trunc: + windows = windows[:, :-1] + if singleton: + windows = windows[0] + return (windows, ratios) if return_ratios else windows + + +def _fftautocorr(x): + """Compute the autocorrelation of a real array and crop the result.""" + N = x.shape[-1] + use_N = fftpack.next_fast_len(2*N-1) + x_fft = np.fft.rfft(x, use_N, axis=-1) + cxy = np.fft.irfft(x_fft * x_fft.conj(), n=use_N)[:, :N] + # Or equivalently (but in most cases slower): + # cxy = np.array([np.convolve(xx, yy[::-1], mode='full') + # for xx, yy in zip(x, x)])[:, N-1:2*N-1] + return cxy + + +_win_equiv_raw = { + ('barthann', 'brthan', 'bth'): (barthann, False), + ('bartlett', 'bart', 'brt'): (bartlett, False), + ('blackman', 'black', 'blk'): (blackman, False), + ('blackmanharris', 'blackharr', 'bkh'): (blackmanharris, False), + ('bohman', 'bman', 'bmn'): (bohman, False), + ('boxcar', 'box', 'ones', + 'rect', 'rectangular'): (boxcar, False), + ('chebwin', 'cheb'): (chebwin, True), + ('cosine', 'halfcosine'): (cosine, False), + ('exponential', 'poisson'): (exponential, True), + ('flattop', 'flat', 'flt'): (flattop, False), + ('gaussian', 'gauss', 'gss'): (gaussian, True), + ('general gaussian', 'general_gaussian', + 'general gauss', 'general_gauss', 'ggs'): (general_gaussian, True), + ('hamming', 'hamm', 'ham'): (hamming, False), + ('hanning', 'hann', 'han'): (hann, False), + ('kaiser', 'ksr'): (kaiser, True), + ('nuttall', 'nutl', 'nut'): (nuttall, False), + ('parzen', 'parz', 'par'): (parzen, False), + ('slepian', 'slep', 'optimal', 'dpss', 'dss'): (slepian, True), + ('triangle', 'triang', 'tri'): (triang, False), + ('tukey', 'tuk'): (tukey, True), +} + +# Fill dict with all valid window name strings +_win_equiv = {} +for k, v in _win_equiv_raw.items(): + for key in k: + _win_equiv[key] = v[0] + +# Keep track of which windows need additional parameters +_needs_param = set() +for k, v in _win_equiv_raw.items(): + if v[1]: + _needs_param.update(k) + + +def get_window(window, Nx, fftbins=True): + """ + Return a window. + + Parameters + ---------- + window : string, float, or tuple + The type of window to create. See below for more details. + Nx : int + The number of samples in the window. + fftbins : bool, optional + If True (default), create a "periodic" window, ready to use with + `ifftshift` and be multiplied by the result of an FFT (see also + `fftpack.fftfreq`). + If False, create a "symmetric" window, for use in filter design. + + Returns + ------- + get_window : ndarray + Returns a window of length `Nx` and type `window` + + Notes + ----- + Window types: + + `boxcar`, `triang`, `blackman`, `hamming`, `hann`, `bartlett`, + `flattop`, `parzen`, `bohman`, `blackmanharris`, `nuttall`, + `barthann`, `kaiser` (needs beta), `gaussian` (needs standard + deviation), `general_gaussian` (needs power, width), `slepian` + (needs width), `dpss` (needs normalized half-bandwidth), + `chebwin` (needs attenuation), `exponential` (needs decay scale), + `tukey` (needs taper fraction) + + If the window requires no parameters, then `window` can be a string. + + If the window requires parameters, then `window` must be a tuple + with the first argument the string name of the window, and the next + arguments the needed parameters. + + If `window` is a floating point number, it is interpreted as the beta + parameter of the `kaiser` window. + + Each of the window types listed above is also the name of + a function that can be called directly to create a window of + that type. + + Examples + -------- + >>> from scipy import signal + >>> signal.get_window('triang', 7) + array([ 0.125, 0.375, 0.625, 0.875, 0.875, 0.625, 0.375]) + >>> signal.get_window(('kaiser', 4.0), 9) + array([ 0.08848053, 0.29425961, 0.56437221, 0.82160913, 0.97885093, + 0.97885093, 0.82160913, 0.56437221, 0.29425961]) + >>> signal.get_window(4.0, 9) + array([ 0.08848053, 0.29425961, 0.56437221, 0.82160913, 0.97885093, + 0.97885093, 0.82160913, 0.56437221, 0.29425961]) + + """ + sym = not fftbins + try: + beta = float(window) + except (TypeError, ValueError): + args = () + if isinstance(window, tuple): + winstr = window[0] + if len(window) > 1: + args = window[1:] + elif isinstance(window, string_types): + if window in _needs_param: + raise ValueError("The '" + window + "' window needs one or " + "more parameters -- pass a tuple.") + else: + winstr = window + else: + raise ValueError("%s as window type is not supported." % + str(type(window))) + + try: + winfunc = _win_equiv[winstr] + except KeyError: + raise ValueError("Unknown window type.") + + params = (Nx,) + args + (sym,) + else: + winfunc = kaiser + params = (Nx, beta, sym) + + return winfunc(*params)