1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993199419951996199719981999200020012002200320042005200620072008200920102011201220132014201520162017201820192020202120222023202420252026202720282029203020312032203320342035203620372038203920402041204220432044204520462047204820492050205120522053205420552056205720582059206020612062206320642065206620672068206920702071207220732074207520762077207820792080208120822083208420852086208720882089209020912092209320942095209620972098209921002101210221032104210521062107 |
- """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 <https://ccrma.stanford.edu/~jos/sasp/Kaiser_DPSS_Windows_Compared.html>`_):
- >>> 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)
|