windows.py 72 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993199419951996199719981999200020012002200320042005200620072008200920102011201220132014201520162017201820192020202120222023202420252026202720282029203020312032203320342035203620372038203920402041204220432044204520462047204820492050205120522053205420552056205720582059206020612062206320642065206620672068206920702071207220732074207520762077207820792080208120822083208420852086208720882089209020912092209320942095209620972098209921002101210221032104210521062107
  1. """The suite of window functions."""
  2. from __future__ import division, print_function, absolute_import
  3. import operator
  4. import warnings
  5. import numpy as np
  6. from scipy import fftpack, linalg, special
  7. from scipy._lib.six import string_types
  8. __all__ = ['boxcar', 'triang', 'parzen', 'bohman', 'blackman', 'nuttall',
  9. 'blackmanharris', 'flattop', 'bartlett', 'hanning', 'barthann',
  10. 'hamming', 'kaiser', 'gaussian', 'general_cosine','general_gaussian',
  11. 'general_hamming', 'chebwin', 'slepian', 'cosine', 'hann',
  12. 'exponential', 'tukey', 'dpss', 'get_window']
  13. def _len_guards(M):
  14. """Handle small or incorrect window lengths"""
  15. if int(M) != M or M < 0:
  16. raise ValueError('Window length M must be a non-negative integer')
  17. return M <= 1
  18. def _extend(M, sym):
  19. """Extend window by 1 sample if needed for DFT-even symmetry"""
  20. if not sym:
  21. return M + 1, True
  22. else:
  23. return M, False
  24. def _truncate(w, needed):
  25. """Truncate window by 1 sample if needed for DFT-even symmetry"""
  26. if needed:
  27. return w[:-1]
  28. else:
  29. return w
  30. def general_cosine(M, a, sym=True):
  31. r"""
  32. Generic weighted sum of cosine terms window
  33. Parameters
  34. ----------
  35. M : int
  36. Number of points in the output window
  37. a : array_like
  38. Sequence of weighting coefficients. This uses the convention of being
  39. centered on the origin, so these will typically all be positive
  40. numbers, not alternating sign.
  41. sym : bool, optional
  42. When True (default), generates a symmetric window, for use in filter
  43. design.
  44. When False, generates a periodic window, for use in spectral analysis.
  45. References
  46. ----------
  47. .. [1] A. Nuttall, "Some windows with very good sidelobe behavior," IEEE
  48. Transactions on Acoustics, Speech, and Signal Processing, vol. 29,
  49. no. 1, pp. 84-91, Feb 1981. :doi:`10.1109/TASSP.1981.1163506`.
  50. .. [2] Heinzel G. et al., "Spectrum and spectral density estimation by the
  51. Discrete Fourier transform (DFT), including a comprehensive list of
  52. window functions and some new flat-top windows", February 15, 2002
  53. https://holometer.fnal.gov/GH_FFT.pdf
  54. Examples
  55. --------
  56. Heinzel describes a flat-top window named "HFT90D" with formula: [2]_
  57. .. math:: w_j = 1 - 1.942604 \cos(z) + 1.340318 \cos(2z)
  58. - 0.440811 \cos(3z) + 0.043097 \cos(4z)
  59. where
  60. .. math:: z = \frac{2 \pi j}{N}, j = 0...N - 1
  61. Since this uses the convention of starting at the origin, to reproduce the
  62. window, we need to convert every other coefficient to a positive number:
  63. >>> HFT90D = [1, 1.942604, 1.340318, 0.440811, 0.043097]
  64. The paper states that the highest sidelobe is at -90.2 dB. Reproduce
  65. Figure 42 by plotting the window and its frequency response, and confirm
  66. the sidelobe level in red:
  67. >>> from scipy.signal.windows import general_cosine
  68. >>> from scipy.fftpack import fft, fftshift
  69. >>> import matplotlib.pyplot as plt
  70. >>> window = general_cosine(1000, HFT90D, sym=False)
  71. >>> plt.plot(window)
  72. >>> plt.title("HFT90D window")
  73. >>> plt.ylabel("Amplitude")
  74. >>> plt.xlabel("Sample")
  75. >>> plt.figure()
  76. >>> A = fft(window, 10000) / (len(window)/2.0)
  77. >>> freq = np.linspace(-0.5, 0.5, len(A))
  78. >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
  79. >>> plt.plot(freq, response)
  80. >>> plt.axis([-50/1000, 50/1000, -140, 0])
  81. >>> plt.title("Frequency response of the HFT90D window")
  82. >>> plt.ylabel("Normalized magnitude [dB]")
  83. >>> plt.xlabel("Normalized frequency [cycles per sample]")
  84. >>> plt.axhline(-90.2, color='red')
  85. >>> plt.show()
  86. """
  87. if _len_guards(M):
  88. return np.ones(M)
  89. M, needs_trunc = _extend(M, sym)
  90. fac = np.linspace(-np.pi, np.pi, M)
  91. w = np.zeros(M)
  92. for k in range(len(a)):
  93. w += a[k] * np.cos(k * fac)
  94. return _truncate(w, needs_trunc)
  95. def boxcar(M, sym=True):
  96. """Return a boxcar or rectangular window.
  97. Also known as a rectangular window or Dirichlet window, this is equivalent
  98. to no window at all.
  99. Parameters
  100. ----------
  101. M : int
  102. Number of points in the output window. If zero or less, an empty
  103. array is returned.
  104. sym : bool, optional
  105. Whether the window is symmetric. (Has no effect for boxcar.)
  106. Returns
  107. -------
  108. w : ndarray
  109. The window, with the maximum value normalized to 1.
  110. Examples
  111. --------
  112. Plot the window and its frequency response:
  113. >>> from scipy import signal
  114. >>> from scipy.fftpack import fft, fftshift
  115. >>> import matplotlib.pyplot as plt
  116. >>> window = signal.boxcar(51)
  117. >>> plt.plot(window)
  118. >>> plt.title("Boxcar window")
  119. >>> plt.ylabel("Amplitude")
  120. >>> plt.xlabel("Sample")
  121. >>> plt.figure()
  122. >>> A = fft(window, 2048) / (len(window)/2.0)
  123. >>> freq = np.linspace(-0.5, 0.5, len(A))
  124. >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
  125. >>> plt.plot(freq, response)
  126. >>> plt.axis([-0.5, 0.5, -120, 0])
  127. >>> plt.title("Frequency response of the boxcar window")
  128. >>> plt.ylabel("Normalized magnitude [dB]")
  129. >>> plt.xlabel("Normalized frequency [cycles per sample]")
  130. """
  131. if _len_guards(M):
  132. return np.ones(M)
  133. M, needs_trunc = _extend(M, sym)
  134. w = np.ones(M, float)
  135. return _truncate(w, needs_trunc)
  136. def triang(M, sym=True):
  137. """Return a triangular window.
  138. Parameters
  139. ----------
  140. M : int
  141. Number of points in the output window. If zero or less, an empty
  142. array is returned.
  143. sym : bool, optional
  144. When True (default), generates a symmetric window, for use in filter
  145. design.
  146. When False, generates a periodic window, for use in spectral analysis.
  147. Returns
  148. -------
  149. w : ndarray
  150. The window, with the maximum value normalized to 1 (though the value 1
  151. does not appear if `M` is even and `sym` is True).
  152. See Also
  153. --------
  154. bartlett : A triangular window that touches zero
  155. Examples
  156. --------
  157. Plot the window and its frequency response:
  158. >>> from scipy import signal
  159. >>> from scipy.fftpack import fft, fftshift
  160. >>> import matplotlib.pyplot as plt
  161. >>> window = signal.triang(51)
  162. >>> plt.plot(window)
  163. >>> plt.title("Triangular window")
  164. >>> plt.ylabel("Amplitude")
  165. >>> plt.xlabel("Sample")
  166. >>> plt.figure()
  167. >>> A = fft(window, 2048) / (len(window)/2.0)
  168. >>> freq = np.linspace(-0.5, 0.5, len(A))
  169. >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
  170. >>> plt.plot(freq, response)
  171. >>> plt.axis([-0.5, 0.5, -120, 0])
  172. >>> plt.title("Frequency response of the triangular window")
  173. >>> plt.ylabel("Normalized magnitude [dB]")
  174. >>> plt.xlabel("Normalized frequency [cycles per sample]")
  175. """
  176. if _len_guards(M):
  177. return np.ones(M)
  178. M, needs_trunc = _extend(M, sym)
  179. n = np.arange(1, (M + 1) // 2 + 1)
  180. if M % 2 == 0:
  181. w = (2 * n - 1.0) / M
  182. w = np.r_[w, w[::-1]]
  183. else:
  184. w = 2 * n / (M + 1.0)
  185. w = np.r_[w, w[-2::-1]]
  186. return _truncate(w, needs_trunc)
  187. def parzen(M, sym=True):
  188. """Return a Parzen window.
  189. Parameters
  190. ----------
  191. M : int
  192. Number of points in the output window. If zero or less, an empty
  193. array is returned.
  194. sym : bool, optional
  195. When True (default), generates a symmetric window, for use in filter
  196. design.
  197. When False, generates a periodic window, for use in spectral analysis.
  198. Returns
  199. -------
  200. w : ndarray
  201. The window, with the maximum value normalized to 1 (though the value 1
  202. does not appear if `M` is even and `sym` is True).
  203. References
  204. ----------
  205. .. [1] E. Parzen, "Mathematical Considerations in the Estimation of
  206. Spectra", Technometrics, Vol. 3, No. 2 (May, 1961), pp. 167-190
  207. Examples
  208. --------
  209. Plot the window and its frequency response:
  210. >>> from scipy import signal
  211. >>> from scipy.fftpack import fft, fftshift
  212. >>> import matplotlib.pyplot as plt
  213. >>> window = signal.parzen(51)
  214. >>> plt.plot(window)
  215. >>> plt.title("Parzen window")
  216. >>> plt.ylabel("Amplitude")
  217. >>> plt.xlabel("Sample")
  218. >>> plt.figure()
  219. >>> A = fft(window, 2048) / (len(window)/2.0)
  220. >>> freq = np.linspace(-0.5, 0.5, len(A))
  221. >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
  222. >>> plt.plot(freq, response)
  223. >>> plt.axis([-0.5, 0.5, -120, 0])
  224. >>> plt.title("Frequency response of the Parzen window")
  225. >>> plt.ylabel("Normalized magnitude [dB]")
  226. >>> plt.xlabel("Normalized frequency [cycles per sample]")
  227. """
  228. if _len_guards(M):
  229. return np.ones(M)
  230. M, needs_trunc = _extend(M, sym)
  231. n = np.arange(-(M - 1) / 2.0, (M - 1) / 2.0 + 0.5, 1.0)
  232. na = np.extract(n < -(M - 1) / 4.0, n)
  233. nb = np.extract(abs(n) <= (M - 1) / 4.0, n)
  234. wa = 2 * (1 - np.abs(na) / (M / 2.0)) ** 3.0
  235. wb = (1 - 6 * (np.abs(nb) / (M / 2.0)) ** 2.0 +
  236. 6 * (np.abs(nb) / (M / 2.0)) ** 3.0)
  237. w = np.r_[wa, wb, wa[::-1]]
  238. return _truncate(w, needs_trunc)
  239. def bohman(M, sym=True):
  240. """Return a Bohman window.
  241. Parameters
  242. ----------
  243. M : int
  244. Number of points in the output window. If zero or less, an empty
  245. array is returned.
  246. sym : bool, optional
  247. When True (default), generates a symmetric window, for use in filter
  248. design.
  249. When False, generates a periodic window, for use in spectral analysis.
  250. Returns
  251. -------
  252. w : ndarray
  253. The window, with the maximum value normalized to 1 (though the value 1
  254. does not appear if `M` is even and `sym` is True).
  255. Examples
  256. --------
  257. Plot the window and its frequency response:
  258. >>> from scipy import signal
  259. >>> from scipy.fftpack import fft, fftshift
  260. >>> import matplotlib.pyplot as plt
  261. >>> window = signal.bohman(51)
  262. >>> plt.plot(window)
  263. >>> plt.title("Bohman window")
  264. >>> plt.ylabel("Amplitude")
  265. >>> plt.xlabel("Sample")
  266. >>> plt.figure()
  267. >>> A = fft(window, 2048) / (len(window)/2.0)
  268. >>> freq = np.linspace(-0.5, 0.5, len(A))
  269. >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
  270. >>> plt.plot(freq, response)
  271. >>> plt.axis([-0.5, 0.5, -120, 0])
  272. >>> plt.title("Frequency response of the Bohman window")
  273. >>> plt.ylabel("Normalized magnitude [dB]")
  274. >>> plt.xlabel("Normalized frequency [cycles per sample]")
  275. """
  276. if _len_guards(M):
  277. return np.ones(M)
  278. M, needs_trunc = _extend(M, sym)
  279. fac = np.abs(np.linspace(-1, 1, M)[1:-1])
  280. w = (1 - fac) * np.cos(np.pi * fac) + 1.0 / np.pi * np.sin(np.pi * fac)
  281. w = np.r_[0, w, 0]
  282. return _truncate(w, needs_trunc)
  283. def blackman(M, sym=True):
  284. r"""
  285. Return a Blackman window.
  286. The Blackman window is a taper formed by using the first three terms of
  287. a summation of cosines. It was designed to have close to the minimal
  288. leakage possible. It is close to optimal, only slightly worse than a
  289. Kaiser window.
  290. Parameters
  291. ----------
  292. M : int
  293. Number of points in the output window. If zero or less, an empty
  294. array is returned.
  295. sym : bool, optional
  296. When True (default), generates a symmetric window, for use in filter
  297. design.
  298. When False, generates a periodic window, for use in spectral analysis.
  299. Returns
  300. -------
  301. w : ndarray
  302. The window, with the maximum value normalized to 1 (though the value 1
  303. does not appear if `M` is even and `sym` is True).
  304. Notes
  305. -----
  306. The Blackman window is defined as
  307. .. math:: w(n) = 0.42 - 0.5 \cos(2\pi n/M) + 0.08 \cos(4\pi n/M)
  308. The "exact Blackman" window was designed to null out the third and fourth
  309. sidelobes, but has discontinuities at the boundaries, resulting in a
  310. 6 dB/oct fall-off. This window is an approximation of the "exact" window,
  311. which does not null the sidelobes as well, but is smooth at the edges,
  312. improving the fall-off rate to 18 dB/oct. [3]_
  313. Most references to the Blackman window come from the signal processing
  314. literature, where it is used as one of many windowing functions for
  315. smoothing values. It is also known as an apodization (which means
  316. "removing the foot", i.e. smoothing discontinuities at the beginning
  317. and end of the sampled signal) or tapering function. It is known as a
  318. "near optimal" tapering function, almost as good (by some measures)
  319. as the Kaiser window.
  320. References
  321. ----------
  322. .. [1] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power
  323. spectra, Dover Publications, New York.
  324. .. [2] Oppenheim, A.V., and R.W. Schafer. Discrete-Time Signal Processing.
  325. Upper Saddle River, NJ: Prentice-Hall, 1999, pp. 468-471.
  326. .. [3] Harris, Fredric J. (Jan 1978). "On the use of Windows for Harmonic
  327. Analysis with the Discrete Fourier Transform". Proceedings of the
  328. IEEE 66 (1): 51-83. :doi:`10.1109/PROC.1978.10837`.
  329. Examples
  330. --------
  331. Plot the window and its frequency response:
  332. >>> from scipy import signal
  333. >>> from scipy.fftpack import fft, fftshift
  334. >>> import matplotlib.pyplot as plt
  335. >>> window = signal.blackman(51)
  336. >>> plt.plot(window)
  337. >>> plt.title("Blackman window")
  338. >>> plt.ylabel("Amplitude")
  339. >>> plt.xlabel("Sample")
  340. >>> plt.figure()
  341. >>> A = fft(window, 2048) / (len(window)/2.0)
  342. >>> freq = np.linspace(-0.5, 0.5, len(A))
  343. >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
  344. >>> plt.plot(freq, response)
  345. >>> plt.axis([-0.5, 0.5, -120, 0])
  346. >>> plt.title("Frequency response of the Blackman window")
  347. >>> plt.ylabel("Normalized magnitude [dB]")
  348. >>> plt.xlabel("Normalized frequency [cycles per sample]")
  349. """
  350. # Docstring adapted from NumPy's blackman function
  351. return general_cosine(M, [0.42, 0.50, 0.08], sym)
  352. def nuttall(M, sym=True):
  353. """Return a minimum 4-term Blackman-Harris window according to Nuttall.
  354. This variation is called "Nuttall4c" by Heinzel. [2]_
  355. Parameters
  356. ----------
  357. M : int
  358. Number of points in the output window. If zero or less, an empty
  359. array is returned.
  360. sym : bool, optional
  361. When True (default), generates a symmetric window, for use in filter
  362. design.
  363. When False, generates a periodic window, for use in spectral analysis.
  364. Returns
  365. -------
  366. w : ndarray
  367. The window, with the maximum value normalized to 1 (though the value 1
  368. does not appear if `M` is even and `sym` is True).
  369. References
  370. ----------
  371. .. [1] A. Nuttall, "Some windows with very good sidelobe behavior," IEEE
  372. Transactions on Acoustics, Speech, and Signal Processing, vol. 29,
  373. no. 1, pp. 84-91, Feb 1981. :doi:`10.1109/TASSP.1981.1163506`.
  374. .. [2] Heinzel G. et al., "Spectrum and spectral density estimation by the
  375. Discrete Fourier transform (DFT), including a comprehensive list of
  376. window functions and some new flat-top windows", February 15, 2002
  377. https://holometer.fnal.gov/GH_FFT.pdf
  378. Examples
  379. --------
  380. Plot the window and its frequency response:
  381. >>> from scipy import signal
  382. >>> from scipy.fftpack import fft, fftshift
  383. >>> import matplotlib.pyplot as plt
  384. >>> window = signal.nuttall(51)
  385. >>> plt.plot(window)
  386. >>> plt.title("Nuttall window")
  387. >>> plt.ylabel("Amplitude")
  388. >>> plt.xlabel("Sample")
  389. >>> plt.figure()
  390. >>> A = fft(window, 2048) / (len(window)/2.0)
  391. >>> freq = np.linspace(-0.5, 0.5, len(A))
  392. >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
  393. >>> plt.plot(freq, response)
  394. >>> plt.axis([-0.5, 0.5, -120, 0])
  395. >>> plt.title("Frequency response of the Nuttall window")
  396. >>> plt.ylabel("Normalized magnitude [dB]")
  397. >>> plt.xlabel("Normalized frequency [cycles per sample]")
  398. """
  399. return general_cosine(M, [0.3635819, 0.4891775, 0.1365995, 0.0106411], sym)
  400. def blackmanharris(M, sym=True):
  401. """Return a minimum 4-term Blackman-Harris window.
  402. Parameters
  403. ----------
  404. M : int
  405. Number of points in the output window. If zero or less, an empty
  406. array is returned.
  407. sym : bool, optional
  408. When True (default), generates a symmetric window, for use in filter
  409. design.
  410. When False, generates a periodic window, for use in spectral analysis.
  411. Returns
  412. -------
  413. w : ndarray
  414. The window, with the maximum value normalized to 1 (though the value 1
  415. does not appear if `M` is even and `sym` is True).
  416. Examples
  417. --------
  418. Plot the window and its frequency response:
  419. >>> from scipy import signal
  420. >>> from scipy.fftpack import fft, fftshift
  421. >>> import matplotlib.pyplot as plt
  422. >>> window = signal.blackmanharris(51)
  423. >>> plt.plot(window)
  424. >>> plt.title("Blackman-Harris window")
  425. >>> plt.ylabel("Amplitude")
  426. >>> plt.xlabel("Sample")
  427. >>> plt.figure()
  428. >>> A = fft(window, 2048) / (len(window)/2.0)
  429. >>> freq = np.linspace(-0.5, 0.5, len(A))
  430. >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
  431. >>> plt.plot(freq, response)
  432. >>> plt.axis([-0.5, 0.5, -120, 0])
  433. >>> plt.title("Frequency response of the Blackman-Harris window")
  434. >>> plt.ylabel("Normalized magnitude [dB]")
  435. >>> plt.xlabel("Normalized frequency [cycles per sample]")
  436. """
  437. return general_cosine(M, [0.35875, 0.48829, 0.14128, 0.01168], sym)
  438. def flattop(M, sym=True):
  439. """Return a flat top window.
  440. Parameters
  441. ----------
  442. M : int
  443. Number of points in the output window. If zero or less, an empty
  444. array is returned.
  445. sym : bool, optional
  446. When True (default), generates a symmetric window, for use in filter
  447. design.
  448. When False, generates a periodic window, for use in spectral analysis.
  449. Returns
  450. -------
  451. w : ndarray
  452. The window, with the maximum value normalized to 1 (though the value 1
  453. does not appear if `M` is even and `sym` is True).
  454. Notes
  455. -----
  456. Flat top windows are used for taking accurate measurements of signal
  457. amplitude in the frequency domain, with minimal scalloping error from the
  458. center of a frequency bin to its edges, compared to others. This is a
  459. 5th-order cosine window, with the 5 terms optimized to make the main lobe
  460. maximally flat. [1]_
  461. References
  462. ----------
  463. .. [1] D'Antona, Gabriele, and A. Ferrero, "Digital Signal Processing for
  464. Measurement Systems", Springer Media, 2006, p. 70
  465. :doi:`10.1007/0-387-28666-7`.
  466. Examples
  467. --------
  468. Plot the window and its frequency response:
  469. >>> from scipy import signal
  470. >>> from scipy.fftpack import fft, fftshift
  471. >>> import matplotlib.pyplot as plt
  472. >>> window = signal.flattop(51)
  473. >>> plt.plot(window)
  474. >>> plt.title("Flat top window")
  475. >>> plt.ylabel("Amplitude")
  476. >>> plt.xlabel("Sample")
  477. >>> plt.figure()
  478. >>> A = fft(window, 2048) / (len(window)/2.0)
  479. >>> freq = np.linspace(-0.5, 0.5, len(A))
  480. >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
  481. >>> plt.plot(freq, response)
  482. >>> plt.axis([-0.5, 0.5, -120, 0])
  483. >>> plt.title("Frequency response of the flat top window")
  484. >>> plt.ylabel("Normalized magnitude [dB]")
  485. >>> plt.xlabel("Normalized frequency [cycles per sample]")
  486. """
  487. a = [0.21557895, 0.41663158, 0.277263158, 0.083578947, 0.006947368]
  488. return general_cosine(M, a, sym)
  489. def bartlett(M, sym=True):
  490. r"""
  491. Return a Bartlett window.
  492. The Bartlett window is very similar to a triangular window, except
  493. that the end points are at zero. It is often used in signal
  494. processing for tapering a signal, without generating too much
  495. ripple in the frequency domain.
  496. Parameters
  497. ----------
  498. M : int
  499. Number of points in the output window. If zero or less, an empty
  500. array is returned.
  501. sym : bool, optional
  502. When True (default), generates a symmetric window, for use in filter
  503. design.
  504. When False, generates a periodic window, for use in spectral analysis.
  505. Returns
  506. -------
  507. w : ndarray
  508. The triangular window, with the first and last samples equal to zero
  509. and the maximum value normalized to 1 (though the value 1 does not
  510. appear if `M` is even and `sym` is True).
  511. See Also
  512. --------
  513. triang : A triangular window that does not touch zero at the ends
  514. Notes
  515. -----
  516. The Bartlett window is defined as
  517. .. math:: w(n) = \frac{2}{M-1} \left(
  518. \frac{M-1}{2} - \left|n - \frac{M-1}{2}\right|
  519. \right)
  520. Most references to the Bartlett window come from the signal
  521. processing literature, where it is used as one of many windowing
  522. functions for smoothing values. Note that convolution with this
  523. window produces linear interpolation. It is also known as an
  524. apodization (which means"removing the foot", i.e. smoothing
  525. discontinuities at the beginning and end of the sampled signal) or
  526. tapering function. The Fourier transform of the Bartlett is the product
  527. of two sinc functions.
  528. Note the excellent discussion in Kanasewich. [2]_
  529. References
  530. ----------
  531. .. [1] M.S. Bartlett, "Periodogram Analysis and Continuous Spectra",
  532. Biometrika 37, 1-16, 1950.
  533. .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics",
  534. The University of Alberta Press, 1975, pp. 109-110.
  535. .. [3] A.V. Oppenheim and R.W. Schafer, "Discrete-Time Signal
  536. Processing", Prentice-Hall, 1999, pp. 468-471.
  537. .. [4] Wikipedia, "Window function",
  538. https://en.wikipedia.org/wiki/Window_function
  539. .. [5] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
  540. "Numerical Recipes", Cambridge University Press, 1986, page 429.
  541. Examples
  542. --------
  543. Plot the window and its frequency response:
  544. >>> from scipy import signal
  545. >>> from scipy.fftpack import fft, fftshift
  546. >>> import matplotlib.pyplot as plt
  547. >>> window = signal.bartlett(51)
  548. >>> plt.plot(window)
  549. >>> plt.title("Bartlett window")
  550. >>> plt.ylabel("Amplitude")
  551. >>> plt.xlabel("Sample")
  552. >>> plt.figure()
  553. >>> A = fft(window, 2048) / (len(window)/2.0)
  554. >>> freq = np.linspace(-0.5, 0.5, len(A))
  555. >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
  556. >>> plt.plot(freq, response)
  557. >>> plt.axis([-0.5, 0.5, -120, 0])
  558. >>> plt.title("Frequency response of the Bartlett window")
  559. >>> plt.ylabel("Normalized magnitude [dB]")
  560. >>> plt.xlabel("Normalized frequency [cycles per sample]")
  561. """
  562. # Docstring adapted from NumPy's bartlett function
  563. if _len_guards(M):
  564. return np.ones(M)
  565. M, needs_trunc = _extend(M, sym)
  566. n = np.arange(0, M)
  567. w = np.where(np.less_equal(n, (M - 1) / 2.0),
  568. 2.0 * n / (M - 1), 2.0 - 2.0 * n / (M - 1))
  569. return _truncate(w, needs_trunc)
  570. def hann(M, sym=True):
  571. r"""
  572. Return a Hann window.
  573. The Hann window is a taper formed by using a raised cosine or sine-squared
  574. with ends that touch zero.
  575. Parameters
  576. ----------
  577. M : int
  578. Number of points in the output window. If zero or less, an empty
  579. array is returned.
  580. sym : bool, optional
  581. When True (default), generates a symmetric window, for use in filter
  582. design.
  583. When False, generates a periodic window, for use in spectral analysis.
  584. Returns
  585. -------
  586. w : ndarray
  587. The window, with the maximum value normalized to 1 (though the value 1
  588. does not appear if `M` is even and `sym` is True).
  589. Notes
  590. -----
  591. The Hann window is defined as
  592. .. math:: w(n) = 0.5 - 0.5 \cos\left(\frac{2\pi{n}}{M-1}\right)
  593. \qquad 0 \leq n \leq M-1
  594. The window was named for Julius von Hann, an Austrian meteorologist. It is
  595. also known as the Cosine Bell. It is sometimes erroneously referred to as
  596. the "Hanning" window, from the use of "hann" as a verb in the original
  597. paper and confusion with the very similar Hamming window.
  598. Most references to the Hann window come from the signal processing
  599. literature, where it is used as one of many windowing functions for
  600. smoothing values. It is also known as an apodization (which means
  601. "removing the foot", i.e. smoothing discontinuities at the beginning
  602. and end of the sampled signal) or tapering function.
  603. References
  604. ----------
  605. .. [1] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power
  606. spectra, Dover Publications, New York.
  607. .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics",
  608. The University of Alberta Press, 1975, pp. 106-108.
  609. .. [3] Wikipedia, "Window function",
  610. https://en.wikipedia.org/wiki/Window_function
  611. .. [4] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
  612. "Numerical Recipes", Cambridge University Press, 1986, page 425.
  613. Examples
  614. --------
  615. Plot the window and its frequency response:
  616. >>> from scipy import signal
  617. >>> from scipy.fftpack import fft, fftshift
  618. >>> import matplotlib.pyplot as plt
  619. >>> window = signal.hann(51)
  620. >>> plt.plot(window)
  621. >>> plt.title("Hann window")
  622. >>> plt.ylabel("Amplitude")
  623. >>> plt.xlabel("Sample")
  624. >>> plt.figure()
  625. >>> A = fft(window, 2048) / (len(window)/2.0)
  626. >>> freq = np.linspace(-0.5, 0.5, len(A))
  627. >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
  628. >>> plt.plot(freq, response)
  629. >>> plt.axis([-0.5, 0.5, -120, 0])
  630. >>> plt.title("Frequency response of the Hann window")
  631. >>> plt.ylabel("Normalized magnitude [dB]")
  632. >>> plt.xlabel("Normalized frequency [cycles per sample]")
  633. """
  634. # Docstring adapted from NumPy's hanning function
  635. return general_hamming(M, 0.5, sym)
  636. @np.deprecate(new_name='scipy.signal.windows.hann')
  637. def hanning(*args, **kwargs):
  638. return hann(*args, **kwargs)
  639. def tukey(M, alpha=0.5, sym=True):
  640. r"""Return a Tukey window, also known as a tapered cosine window.
  641. Parameters
  642. ----------
  643. M : int
  644. Number of points in the output window. If zero or less, an empty
  645. array is returned.
  646. alpha : float, optional
  647. Shape parameter of the Tukey window, representing the fraction of the
  648. window inside the cosine tapered region.
  649. If zero, the Tukey window is equivalent to a rectangular window.
  650. If one, the Tukey window is equivalent to a Hann window.
  651. sym : bool, optional
  652. When True (default), generates a symmetric window, for use in filter
  653. design.
  654. When False, generates a periodic window, for use in spectral analysis.
  655. Returns
  656. -------
  657. w : ndarray
  658. The window, with the maximum value normalized to 1 (though the value 1
  659. does not appear if `M` is even and `sym` is True).
  660. References
  661. ----------
  662. .. [1] Harris, Fredric J. (Jan 1978). "On the use of Windows for Harmonic
  663. Analysis with the Discrete Fourier Transform". Proceedings of the
  664. IEEE 66 (1): 51-83. :doi:`10.1109/PROC.1978.10837`
  665. .. [2] Wikipedia, "Window function",
  666. https://en.wikipedia.org/wiki/Window_function#Tukey_window
  667. Examples
  668. --------
  669. Plot the window and its frequency response:
  670. >>> from scipy import signal
  671. >>> from scipy.fftpack import fft, fftshift
  672. >>> import matplotlib.pyplot as plt
  673. >>> window = signal.tukey(51)
  674. >>> plt.plot(window)
  675. >>> plt.title("Tukey window")
  676. >>> plt.ylabel("Amplitude")
  677. >>> plt.xlabel("Sample")
  678. >>> plt.ylim([0, 1.1])
  679. >>> plt.figure()
  680. >>> A = fft(window, 2048) / (len(window)/2.0)
  681. >>> freq = np.linspace(-0.5, 0.5, len(A))
  682. >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
  683. >>> plt.plot(freq, response)
  684. >>> plt.axis([-0.5, 0.5, -120, 0])
  685. >>> plt.title("Frequency response of the Tukey window")
  686. >>> plt.ylabel("Normalized magnitude [dB]")
  687. >>> plt.xlabel("Normalized frequency [cycles per sample]")
  688. """
  689. if _len_guards(M):
  690. return np.ones(M)
  691. if alpha <= 0:
  692. return np.ones(M, 'd')
  693. elif alpha >= 1.0:
  694. return hann(M, sym=sym)
  695. M, needs_trunc = _extend(M, sym)
  696. n = np.arange(0, M)
  697. width = int(np.floor(alpha*(M-1)/2.0))
  698. n1 = n[0:width+1]
  699. n2 = n[width+1:M-width-1]
  700. n3 = n[M-width-1:]
  701. w1 = 0.5 * (1 + np.cos(np.pi * (-1 + 2.0*n1/alpha/(M-1))))
  702. w2 = np.ones(n2.shape)
  703. w3 = 0.5 * (1 + np.cos(np.pi * (-2.0/alpha + 1 + 2.0*n3/alpha/(M-1))))
  704. w = np.concatenate((w1, w2, w3))
  705. return _truncate(w, needs_trunc)
  706. def barthann(M, sym=True):
  707. """Return a modified Bartlett-Hann window.
  708. Parameters
  709. ----------
  710. M : int
  711. Number of points in the output window. If zero or less, an empty
  712. array is returned.
  713. sym : bool, optional
  714. When True (default), generates a symmetric window, for use in filter
  715. design.
  716. When False, generates a periodic window, for use in spectral analysis.
  717. Returns
  718. -------
  719. w : ndarray
  720. The window, with the maximum value normalized to 1 (though the value 1
  721. does not appear if `M` is even and `sym` is True).
  722. Examples
  723. --------
  724. Plot the window and its frequency response:
  725. >>> from scipy import signal
  726. >>> from scipy.fftpack import fft, fftshift
  727. >>> import matplotlib.pyplot as plt
  728. >>> window = signal.barthann(51)
  729. >>> plt.plot(window)
  730. >>> plt.title("Bartlett-Hann window")
  731. >>> plt.ylabel("Amplitude")
  732. >>> plt.xlabel("Sample")
  733. >>> plt.figure()
  734. >>> A = fft(window, 2048) / (len(window)/2.0)
  735. >>> freq = np.linspace(-0.5, 0.5, len(A))
  736. >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
  737. >>> plt.plot(freq, response)
  738. >>> plt.axis([-0.5, 0.5, -120, 0])
  739. >>> plt.title("Frequency response of the Bartlett-Hann window")
  740. >>> plt.ylabel("Normalized magnitude [dB]")
  741. >>> plt.xlabel("Normalized frequency [cycles per sample]")
  742. """
  743. if _len_guards(M):
  744. return np.ones(M)
  745. M, needs_trunc = _extend(M, sym)
  746. n = np.arange(0, M)
  747. fac = np.abs(n / (M - 1.0) - 0.5)
  748. w = 0.62 - 0.48 * fac + 0.38 * np.cos(2 * np.pi * fac)
  749. return _truncate(w, needs_trunc)
  750. def general_hamming(M, alpha, sym=True):
  751. r"""Return a generalized Hamming window.
  752. The generalized Hamming window is constructed by multiplying a rectangular
  753. window by one period of a cosine function [1]_.
  754. Parameters
  755. ----------
  756. M : int
  757. Number of points in the output window. If zero or less, an empty
  758. array is returned.
  759. alpha : float
  760. The window coefficient, :math:`\alpha`
  761. sym : bool, optional
  762. When True (default), generates a symmetric window, for use in filter
  763. design.
  764. When False, generates a periodic window, for use in spectral analysis.
  765. Returns
  766. -------
  767. w : ndarray
  768. The window, with the maximum value normalized to 1 (though the value 1
  769. does not appear if `M` is even and `sym` is True).
  770. Notes
  771. -----
  772. The generalized Hamming window is defined as
  773. .. math:: w(n) = \alpha - \left(1 - \alpha\right) \cos\left(\frac{2\pi{n}}{M-1}\right)
  774. \qquad 0 \leq n \leq M-1
  775. Both the common Hamming window and Hann window are special cases of the
  776. generalized Hamming window with :math:`\alpha` = 0.54 and :math:`\alpha` =
  777. 0.5, respectively [2]_.
  778. See Also
  779. --------
  780. hamming, hann
  781. Examples
  782. --------
  783. The Sentinel-1A/B Instrument Processing Facility uses generalized Hamming
  784. windows in the processing of spaceborne Synthetic Aperture Radar (SAR)
  785. data [3]_. The facility uses various values for the :math:`\alpha`
  786. parameter based on operating mode of the SAR instrument. Some common
  787. :math:`\alpha` values include 0.75, 0.7 and 0.52 [4]_. As an example, we
  788. plot these different windows.
  789. >>> from scipy.signal.windows import general_hamming
  790. >>> from scipy.fftpack import fft, fftshift
  791. >>> import matplotlib.pyplot as plt
  792. >>> fig1, spatial_plot = plt.subplots()
  793. >>> spatial_plot.set_title("Generalized Hamming Windows")
  794. >>> spatial_plot.set_ylabel("Amplitude")
  795. >>> spatial_plot.set_xlabel("Sample")
  796. >>> fig2, freq_plot = plt.subplots()
  797. >>> freq_plot.set_title("Frequency Responses")
  798. >>> freq_plot.set_ylabel("Normalized magnitude [dB]")
  799. >>> freq_plot.set_xlabel("Normalized frequency [cycles per sample]")
  800. >>> for alpha in [0.75, 0.7, 0.52]:
  801. ... window = general_hamming(41, alpha)
  802. ... spatial_plot.plot(window, label="{:.2f}".format(alpha))
  803. ... A = fft(window, 2048) / (len(window)/2.0)
  804. ... freq = np.linspace(-0.5, 0.5, len(A))
  805. ... response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
  806. ... freq_plot.plot(freq, response, label="{:.2f}".format(alpha))
  807. >>> freq_plot.legend(loc="upper right")
  808. >>> spatial_plot.legend(loc="upper right")
  809. References
  810. ----------
  811. .. [1] DSPRelated, "Generalized Hamming Window Family",
  812. https://www.dsprelated.com/freebooks/sasp/Generalized_Hamming_Window_Family.html
  813. .. [2] Wikipedia, "Window function",
  814. https://en.wikipedia.org/wiki/Window_function
  815. .. [3] Riccardo Piantanida ESA, "Sentinel-1 Level 1 Detailed Algorithm
  816. Definition",
  817. https://sentinel.esa.int/documents/247904/1877131/Sentinel-1-Level-1-Detailed-Algorithm-Definition
  818. .. [4] Matthieu Bourbigot ESA, "Sentinel-1 Product Definition",
  819. https://sentinel.esa.int/documents/247904/1877131/Sentinel-1-Product-Definition
  820. """
  821. return general_cosine(M, [alpha, 1. - alpha], sym)
  822. def hamming(M, sym=True):
  823. r"""Return a Hamming window.
  824. The Hamming window is a taper formed by using a raised cosine with
  825. non-zero endpoints, optimized to minimize the nearest side lobe.
  826. Parameters
  827. ----------
  828. M : int
  829. Number of points in the output window. If zero or less, an empty
  830. array is returned.
  831. sym : bool, optional
  832. When True (default), generates a symmetric window, for use in filter
  833. design.
  834. When False, generates a periodic window, for use in spectral analysis.
  835. Returns
  836. -------
  837. w : ndarray
  838. The window, with the maximum value normalized to 1 (though the value 1
  839. does not appear if `M` is even and `sym` is True).
  840. Notes
  841. -----
  842. The Hamming window is defined as
  843. .. math:: w(n) = 0.54 - 0.46 \cos\left(\frac{2\pi{n}}{M-1}\right)
  844. \qquad 0 \leq n \leq M-1
  845. The Hamming was named for R. W. Hamming, an associate of J. W. Tukey and
  846. is described in Blackman and Tukey. It was recommended for smoothing the
  847. truncated autocovariance function in the time domain.
  848. Most references to the Hamming window come from the signal processing
  849. literature, where it is used as one of many windowing functions for
  850. smoothing values. It is also known as an apodization (which means
  851. "removing the foot", i.e. smoothing discontinuities at the beginning
  852. and end of the sampled signal) or tapering function.
  853. References
  854. ----------
  855. .. [1] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power
  856. spectra, Dover Publications, New York.
  857. .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", The
  858. University of Alberta Press, 1975, pp. 109-110.
  859. .. [3] Wikipedia, "Window function",
  860. https://en.wikipedia.org/wiki/Window_function
  861. .. [4] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
  862. "Numerical Recipes", Cambridge University Press, 1986, page 425.
  863. Examples
  864. --------
  865. Plot the window and its frequency response:
  866. >>> from scipy import signal
  867. >>> from scipy.fftpack import fft, fftshift
  868. >>> import matplotlib.pyplot as plt
  869. >>> window = signal.hamming(51)
  870. >>> plt.plot(window)
  871. >>> plt.title("Hamming window")
  872. >>> plt.ylabel("Amplitude")
  873. >>> plt.xlabel("Sample")
  874. >>> plt.figure()
  875. >>> A = fft(window, 2048) / (len(window)/2.0)
  876. >>> freq = np.linspace(-0.5, 0.5, len(A))
  877. >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
  878. >>> plt.plot(freq, response)
  879. >>> plt.axis([-0.5, 0.5, -120, 0])
  880. >>> plt.title("Frequency response of the Hamming window")
  881. >>> plt.ylabel("Normalized magnitude [dB]")
  882. >>> plt.xlabel("Normalized frequency [cycles per sample]")
  883. """
  884. # Docstring adapted from NumPy's hamming function
  885. return general_hamming(M, 0.54, sym)
  886. def kaiser(M, beta, sym=True):
  887. r"""Return a Kaiser window.
  888. The Kaiser window is a taper formed by using a Bessel function.
  889. Parameters
  890. ----------
  891. M : int
  892. Number of points in the output window. If zero or less, an empty
  893. array is returned.
  894. beta : float
  895. Shape parameter, determines trade-off between main-lobe width and
  896. side lobe level. As beta gets large, the window narrows.
  897. sym : bool, optional
  898. When True (default), generates a symmetric window, for use in filter
  899. design.
  900. When False, generates a periodic window, for use in spectral analysis.
  901. Returns
  902. -------
  903. w : ndarray
  904. The window, with the maximum value normalized to 1 (though the value 1
  905. does not appear if `M` is even and `sym` is True).
  906. Notes
  907. -----
  908. The Kaiser window is defined as
  909. .. math:: w(n) = I_0\left( \beta \sqrt{1-\frac{4n^2}{(M-1)^2}}
  910. \right)/I_0(\beta)
  911. with
  912. .. math:: \quad -\frac{M-1}{2} \leq n \leq \frac{M-1}{2},
  913. where :math:`I_0` is the modified zeroth-order Bessel function.
  914. The Kaiser was named for Jim Kaiser, who discovered a simple approximation
  915. to the DPSS window based on Bessel functions.
  916. The Kaiser window is a very good approximation to the Digital Prolate
  917. Spheroidal Sequence, or Slepian window, which is the transform which
  918. maximizes the energy in the main lobe of the window relative to total
  919. energy.
  920. The Kaiser can approximate other windows by varying the beta parameter.
  921. (Some literature uses alpha = beta/pi.) [4]_
  922. ==== =======================
  923. beta Window shape
  924. ==== =======================
  925. 0 Rectangular
  926. 5 Similar to a Hamming
  927. 6 Similar to a Hann
  928. 8.6 Similar to a Blackman
  929. ==== =======================
  930. A beta value of 14 is probably a good starting point. Note that as beta
  931. gets large, the window narrows, and so the number of samples needs to be
  932. large enough to sample the increasingly narrow spike, otherwise NaNs will
  933. be returned.
  934. Most references to the Kaiser window come from the signal processing
  935. literature, where it is used as one of many windowing functions for
  936. smoothing values. It is also known as an apodization (which means
  937. "removing the foot", i.e. smoothing discontinuities at the beginning
  938. and end of the sampled signal) or tapering function.
  939. References
  940. ----------
  941. .. [1] J. F. Kaiser, "Digital Filters" - Ch 7 in "Systems analysis by
  942. digital computer", Editors: F.F. Kuo and J.F. Kaiser, p 218-285.
  943. John Wiley and Sons, New York, (1966).
  944. .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", The
  945. University of Alberta Press, 1975, pp. 177-178.
  946. .. [3] Wikipedia, "Window function",
  947. https://en.wikipedia.org/wiki/Window_function
  948. .. [4] F. J. Harris, "On the use of windows for harmonic analysis with the
  949. discrete Fourier transform," Proceedings of the IEEE, vol. 66,
  950. no. 1, pp. 51-83, Jan. 1978. :doi:`10.1109/PROC.1978.10837`.
  951. Examples
  952. --------
  953. Plot the window and its frequency response:
  954. >>> from scipy import signal
  955. >>> from scipy.fftpack import fft, fftshift
  956. >>> import matplotlib.pyplot as plt
  957. >>> window = signal.kaiser(51, beta=14)
  958. >>> plt.plot(window)
  959. >>> plt.title(r"Kaiser window ($\beta$=14)")
  960. >>> plt.ylabel("Amplitude")
  961. >>> plt.xlabel("Sample")
  962. >>> plt.figure()
  963. >>> A = fft(window, 2048) / (len(window)/2.0)
  964. >>> freq = np.linspace(-0.5, 0.5, len(A))
  965. >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
  966. >>> plt.plot(freq, response)
  967. >>> plt.axis([-0.5, 0.5, -120, 0])
  968. >>> plt.title(r"Frequency response of the Kaiser window ($\beta$=14)")
  969. >>> plt.ylabel("Normalized magnitude [dB]")
  970. >>> plt.xlabel("Normalized frequency [cycles per sample]")
  971. """
  972. # Docstring adapted from NumPy's kaiser function
  973. if _len_guards(M):
  974. return np.ones(M)
  975. M, needs_trunc = _extend(M, sym)
  976. n = np.arange(0, M)
  977. alpha = (M - 1) / 2.0
  978. w = (special.i0(beta * np.sqrt(1 - ((n - alpha) / alpha) ** 2.0)) /
  979. special.i0(beta))
  980. return _truncate(w, needs_trunc)
  981. def gaussian(M, std, sym=True):
  982. r"""Return a Gaussian window.
  983. Parameters
  984. ----------
  985. M : int
  986. Number of points in the output window. If zero or less, an empty
  987. array is returned.
  988. std : float
  989. The standard deviation, sigma.
  990. sym : bool, optional
  991. When True (default), generates a symmetric window, for use in filter
  992. design.
  993. When False, generates a periodic window, for use in spectral analysis.
  994. Returns
  995. -------
  996. w : ndarray
  997. The window, with the maximum value normalized to 1 (though the value 1
  998. does not appear if `M` is even and `sym` is True).
  999. Notes
  1000. -----
  1001. The Gaussian window is defined as
  1002. .. math:: w(n) = e^{ -\frac{1}{2}\left(\frac{n}{\sigma}\right)^2 }
  1003. Examples
  1004. --------
  1005. Plot the window and its frequency response:
  1006. >>> from scipy import signal
  1007. >>> from scipy.fftpack import fft, fftshift
  1008. >>> import matplotlib.pyplot as plt
  1009. >>> window = signal.gaussian(51, std=7)
  1010. >>> plt.plot(window)
  1011. >>> plt.title(r"Gaussian window ($\sigma$=7)")
  1012. >>> plt.ylabel("Amplitude")
  1013. >>> plt.xlabel("Sample")
  1014. >>> plt.figure()
  1015. >>> A = fft(window, 2048) / (len(window)/2.0)
  1016. >>> freq = np.linspace(-0.5, 0.5, len(A))
  1017. >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
  1018. >>> plt.plot(freq, response)
  1019. >>> plt.axis([-0.5, 0.5, -120, 0])
  1020. >>> plt.title(r"Frequency response of the Gaussian window ($\sigma$=7)")
  1021. >>> plt.ylabel("Normalized magnitude [dB]")
  1022. >>> plt.xlabel("Normalized frequency [cycles per sample]")
  1023. """
  1024. if _len_guards(M):
  1025. return np.ones(M)
  1026. M, needs_trunc = _extend(M, sym)
  1027. n = np.arange(0, M) - (M - 1.0) / 2.0
  1028. sig2 = 2 * std * std
  1029. w = np.exp(-n ** 2 / sig2)
  1030. return _truncate(w, needs_trunc)
  1031. def general_gaussian(M, p, sig, sym=True):
  1032. r"""Return a window with a generalized Gaussian shape.
  1033. Parameters
  1034. ----------
  1035. M : int
  1036. Number of points in the output window. If zero or less, an empty
  1037. array is returned.
  1038. p : float
  1039. Shape parameter. p = 1 is identical to `gaussian`, p = 0.5 is
  1040. the same shape as the Laplace distribution.
  1041. sig : float
  1042. The standard deviation, sigma.
  1043. sym : bool, optional
  1044. When True (default), generates a symmetric window, for use in filter
  1045. design.
  1046. When False, generates a periodic window, for use in spectral analysis.
  1047. Returns
  1048. -------
  1049. w : ndarray
  1050. The window, with the maximum value normalized to 1 (though the value 1
  1051. does not appear if `M` is even and `sym` is True).
  1052. Notes
  1053. -----
  1054. The generalized Gaussian window is defined as
  1055. .. math:: w(n) = e^{ -\frac{1}{2}\left|\frac{n}{\sigma}\right|^{2p} }
  1056. the half-power point is at
  1057. .. math:: (2 \log(2))^{1/(2 p)} \sigma
  1058. Examples
  1059. --------
  1060. Plot the window and its frequency response:
  1061. >>> from scipy import signal
  1062. >>> from scipy.fftpack import fft, fftshift
  1063. >>> import matplotlib.pyplot as plt
  1064. >>> window = signal.general_gaussian(51, p=1.5, sig=7)
  1065. >>> plt.plot(window)
  1066. >>> plt.title(r"Generalized Gaussian window (p=1.5, $\sigma$=7)")
  1067. >>> plt.ylabel("Amplitude")
  1068. >>> plt.xlabel("Sample")
  1069. >>> plt.figure()
  1070. >>> A = fft(window, 2048) / (len(window)/2.0)
  1071. >>> freq = np.linspace(-0.5, 0.5, len(A))
  1072. >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
  1073. >>> plt.plot(freq, response)
  1074. >>> plt.axis([-0.5, 0.5, -120, 0])
  1075. >>> plt.title(r"Freq. resp. of the gen. Gaussian "
  1076. ... "window (p=1.5, $\sigma$=7)")
  1077. >>> plt.ylabel("Normalized magnitude [dB]")
  1078. >>> plt.xlabel("Normalized frequency [cycles per sample]")
  1079. """
  1080. if _len_guards(M):
  1081. return np.ones(M)
  1082. M, needs_trunc = _extend(M, sym)
  1083. n = np.arange(0, M) - (M - 1.0) / 2.0
  1084. w = np.exp(-0.5 * np.abs(n / sig) ** (2 * p))
  1085. return _truncate(w, needs_trunc)
  1086. # `chebwin` contributed by Kumar Appaiah.
  1087. def chebwin(M, at, sym=True):
  1088. r"""Return a Dolph-Chebyshev window.
  1089. Parameters
  1090. ----------
  1091. M : int
  1092. Number of points in the output window. If zero or less, an empty
  1093. array is returned.
  1094. at : float
  1095. Attenuation (in dB).
  1096. sym : bool, optional
  1097. When True (default), generates a symmetric window, for use in filter
  1098. design.
  1099. When False, generates a periodic window, for use in spectral analysis.
  1100. Returns
  1101. -------
  1102. w : ndarray
  1103. The window, with the maximum value always normalized to 1
  1104. Notes
  1105. -----
  1106. This window optimizes for the narrowest main lobe width for a given order
  1107. `M` and sidelobe equiripple attenuation `at`, using Chebyshev
  1108. polynomials. It was originally developed by Dolph to optimize the
  1109. directionality of radio antenna arrays.
  1110. Unlike most windows, the Dolph-Chebyshev is defined in terms of its
  1111. frequency response:
  1112. .. math:: W(k) = \frac
  1113. {\cos\{M \cos^{-1}[\beta \cos(\frac{\pi k}{M})]\}}
  1114. {\cosh[M \cosh^{-1}(\beta)]}
  1115. where
  1116. .. math:: \beta = \cosh \left [\frac{1}{M}
  1117. \cosh^{-1}(10^\frac{A}{20}) \right ]
  1118. and 0 <= abs(k) <= M-1. A is the attenuation in decibels (`at`).
  1119. The time domain window is then generated using the IFFT, so
  1120. power-of-two `M` are the fastest to generate, and prime number `M` are
  1121. the slowest.
  1122. The equiripple condition in the frequency domain creates impulses in the
  1123. time domain, which appear at the ends of the window.
  1124. References
  1125. ----------
  1126. .. [1] C. Dolph, "A current distribution for broadside arrays which
  1127. optimizes the relationship between beam width and side-lobe level",
  1128. Proceedings of the IEEE, Vol. 34, Issue 6
  1129. .. [2] Peter Lynch, "The Dolph-Chebyshev Window: A Simple Optimal Filter",
  1130. American Meteorological Society (April 1997)
  1131. http://mathsci.ucd.ie/~plynch/Publications/Dolph.pdf
  1132. .. [3] F. J. Harris, "On the use of windows for harmonic analysis with the
  1133. discrete Fourier transforms", Proceedings of the IEEE, Vol. 66,
  1134. No. 1, January 1978
  1135. Examples
  1136. --------
  1137. Plot the window and its frequency response:
  1138. >>> from scipy import signal
  1139. >>> from scipy.fftpack import fft, fftshift
  1140. >>> import matplotlib.pyplot as plt
  1141. >>> window = signal.chebwin(51, at=100)
  1142. >>> plt.plot(window)
  1143. >>> plt.title("Dolph-Chebyshev window (100 dB)")
  1144. >>> plt.ylabel("Amplitude")
  1145. >>> plt.xlabel("Sample")
  1146. >>> plt.figure()
  1147. >>> A = fft(window, 2048) / (len(window)/2.0)
  1148. >>> freq = np.linspace(-0.5, 0.5, len(A))
  1149. >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
  1150. >>> plt.plot(freq, response)
  1151. >>> plt.axis([-0.5, 0.5, -120, 0])
  1152. >>> plt.title("Frequency response of the Dolph-Chebyshev window (100 dB)")
  1153. >>> plt.ylabel("Normalized magnitude [dB]")
  1154. >>> plt.xlabel("Normalized frequency [cycles per sample]")
  1155. """
  1156. if np.abs(at) < 45:
  1157. warnings.warn("This window is not suitable for spectral analysis "
  1158. "for attenuation values lower than about 45dB because "
  1159. "the equivalent noise bandwidth of a Chebyshev window "
  1160. "does not grow monotonically with increasing sidelobe "
  1161. "attenuation when the attenuation is smaller than "
  1162. "about 45 dB.")
  1163. if _len_guards(M):
  1164. return np.ones(M)
  1165. M, needs_trunc = _extend(M, sym)
  1166. # compute the parameter beta
  1167. order = M - 1.0
  1168. beta = np.cosh(1.0 / order * np.arccosh(10 ** (np.abs(at) / 20.)))
  1169. k = np.r_[0:M] * 1.0
  1170. x = beta * np.cos(np.pi * k / M)
  1171. # Find the window's DFT coefficients
  1172. # Use analytic definition of Chebyshev polynomial instead of expansion
  1173. # from scipy.special. Using the expansion in scipy.special leads to errors.
  1174. p = np.zeros(x.shape)
  1175. p[x > 1] = np.cosh(order * np.arccosh(x[x > 1]))
  1176. p[x < -1] = (2 * (M % 2) - 1) * np.cosh(order * np.arccosh(-x[x < -1]))
  1177. p[np.abs(x) <= 1] = np.cos(order * np.arccos(x[np.abs(x) <= 1]))
  1178. # Appropriate IDFT and filling up
  1179. # depending on even/odd M
  1180. if M % 2:
  1181. w = np.real(fftpack.fft(p))
  1182. n = (M + 1) // 2
  1183. w = w[:n]
  1184. w = np.concatenate((w[n - 1:0:-1], w))
  1185. else:
  1186. p = p * np.exp(1.j * np.pi / M * np.r_[0:M])
  1187. w = np.real(fftpack.fft(p))
  1188. n = M // 2 + 1
  1189. w = np.concatenate((w[n - 1:0:-1], w[1:n]))
  1190. w = w / max(w)
  1191. return _truncate(w, needs_trunc)
  1192. def slepian(M, width, sym=True):
  1193. """Return a digital Slepian (DPSS) window.
  1194. Used to maximize the energy concentration in the main lobe. Also called
  1195. the digital prolate spheroidal sequence (DPSS).
  1196. .. note:: Deprecated in SciPy 1.1.
  1197. `slepian` will be removed in a future version of SciPy, it is
  1198. replaced by `dpss`, which uses the standard definition of a
  1199. digital Slepian window.
  1200. Parameters
  1201. ----------
  1202. M : int
  1203. Number of points in the output window. If zero or less, an empty
  1204. array is returned.
  1205. width : float
  1206. Bandwidth
  1207. sym : bool, optional
  1208. When True (default), generates a symmetric window, for use in filter
  1209. design.
  1210. When False, generates a periodic window, for use in spectral analysis.
  1211. Returns
  1212. -------
  1213. w : ndarray
  1214. The window, with the maximum value always normalized to 1
  1215. See Also
  1216. --------
  1217. dpss
  1218. References
  1219. ----------
  1220. .. [1] D. Slepian & H. O. Pollak: "Prolate spheroidal wave functions,
  1221. Fourier analysis and uncertainty-I," Bell Syst. Tech. J., vol.40,
  1222. pp.43-63, 1961. https://archive.org/details/bstj40-1-43
  1223. .. [2] H. J. Landau & H. O. Pollak: "Prolate spheroidal wave functions,
  1224. Fourier analysis and uncertainty-II," Bell Syst. Tech. J. , vol.40,
  1225. pp.65-83, 1961. https://archive.org/details/bstj40-1-65
  1226. Examples
  1227. --------
  1228. Plot the window and its frequency response:
  1229. >>> from scipy import signal
  1230. >>> from scipy.fftpack import fft, fftshift
  1231. >>> import matplotlib.pyplot as plt
  1232. >>> window = signal.slepian(51, width=0.3)
  1233. >>> plt.plot(window)
  1234. >>> plt.title("Slepian (DPSS) window (BW=0.3)")
  1235. >>> plt.ylabel("Amplitude")
  1236. >>> plt.xlabel("Sample")
  1237. >>> plt.figure()
  1238. >>> A = fft(window, 2048) / (len(window)/2.0)
  1239. >>> freq = np.linspace(-0.5, 0.5, len(A))
  1240. >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
  1241. >>> plt.plot(freq, response)
  1242. >>> plt.axis([-0.5, 0.5, -120, 0])
  1243. >>> plt.title("Frequency response of the Slepian window (BW=0.3)")
  1244. >>> plt.ylabel("Normalized magnitude [dB]")
  1245. >>> plt.xlabel("Normalized frequency [cycles per sample]")
  1246. """
  1247. warnings.warn('slepian is deprecated and will be removed in a future '
  1248. 'version, use dpss instead', DeprecationWarning)
  1249. if _len_guards(M):
  1250. return np.ones(M)
  1251. M, needs_trunc = _extend(M, sym)
  1252. # our width is the full bandwidth
  1253. width = width / 2
  1254. # to match the old version
  1255. width = width / 2
  1256. m = np.arange(M, dtype='d')
  1257. H = np.zeros((2, M))
  1258. H[0, 1:] = m[1:] * (M - m[1:]) / 2
  1259. H[1, :] = ((M - 1 - 2 * m) / 2)**2 * np.cos(2 * np.pi * width)
  1260. _, win = linalg.eig_banded(H, select='i', select_range=(M-1, M-1))
  1261. win = win.ravel() / win.max()
  1262. return _truncate(win, needs_trunc)
  1263. def cosine(M, sym=True):
  1264. """Return a window with a simple cosine shape.
  1265. Parameters
  1266. ----------
  1267. M : int
  1268. Number of points in the output window. If zero or less, an empty
  1269. array is returned.
  1270. sym : bool, optional
  1271. When True (default), generates a symmetric window, for use in filter
  1272. design.
  1273. When False, generates a periodic window, for use in spectral analysis.
  1274. Returns
  1275. -------
  1276. w : ndarray
  1277. The window, with the maximum value normalized to 1 (though the value 1
  1278. does not appear if `M` is even and `sym` is True).
  1279. Notes
  1280. -----
  1281. .. versionadded:: 0.13.0
  1282. Examples
  1283. --------
  1284. Plot the window and its frequency response:
  1285. >>> from scipy import signal
  1286. >>> from scipy.fftpack import fft, fftshift
  1287. >>> import matplotlib.pyplot as plt
  1288. >>> window = signal.cosine(51)
  1289. >>> plt.plot(window)
  1290. >>> plt.title("Cosine window")
  1291. >>> plt.ylabel("Amplitude")
  1292. >>> plt.xlabel("Sample")
  1293. >>> plt.figure()
  1294. >>> A = fft(window, 2048) / (len(window)/2.0)
  1295. >>> freq = np.linspace(-0.5, 0.5, len(A))
  1296. >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
  1297. >>> plt.plot(freq, response)
  1298. >>> plt.axis([-0.5, 0.5, -120, 0])
  1299. >>> plt.title("Frequency response of the cosine window")
  1300. >>> plt.ylabel("Normalized magnitude [dB]")
  1301. >>> plt.xlabel("Normalized frequency [cycles per sample]")
  1302. >>> plt.show()
  1303. """
  1304. if _len_guards(M):
  1305. return np.ones(M)
  1306. M, needs_trunc = _extend(M, sym)
  1307. w = np.sin(np.pi / M * (np.arange(0, M) + .5))
  1308. return _truncate(w, needs_trunc)
  1309. def exponential(M, center=None, tau=1., sym=True):
  1310. r"""Return an exponential (or Poisson) window.
  1311. Parameters
  1312. ----------
  1313. M : int
  1314. Number of points in the output window. If zero or less, an empty
  1315. array is returned.
  1316. center : float, optional
  1317. Parameter defining the center location of the window function.
  1318. The default value if not given is ``center = (M-1) / 2``. This
  1319. parameter must take its default value for symmetric windows.
  1320. tau : float, optional
  1321. Parameter defining the decay. For ``center = 0`` use
  1322. ``tau = -(M-1) / ln(x)`` if ``x`` is the fraction of the window
  1323. remaining at the end.
  1324. sym : bool, optional
  1325. When True (default), generates a symmetric window, for use in filter
  1326. design.
  1327. When False, generates a periodic window, for use in spectral analysis.
  1328. Returns
  1329. -------
  1330. w : ndarray
  1331. The window, with the maximum value normalized to 1 (though the value 1
  1332. does not appear if `M` is even and `sym` is True).
  1333. Notes
  1334. -----
  1335. The Exponential window is defined as
  1336. .. math:: w(n) = e^{-|n-center| / \tau}
  1337. References
  1338. ----------
  1339. S. Gade and H. Herlufsen, "Windows to FFT analysis (Part I)",
  1340. Technical Review 3, Bruel & Kjaer, 1987.
  1341. Examples
  1342. --------
  1343. Plot the symmetric window and its frequency response:
  1344. >>> from scipy import signal
  1345. >>> from scipy.fftpack import fft, fftshift
  1346. >>> import matplotlib.pyplot as plt
  1347. >>> M = 51
  1348. >>> tau = 3.0
  1349. >>> window = signal.exponential(M, tau=tau)
  1350. >>> plt.plot(window)
  1351. >>> plt.title("Exponential Window (tau=3.0)")
  1352. >>> plt.ylabel("Amplitude")
  1353. >>> plt.xlabel("Sample")
  1354. >>> plt.figure()
  1355. >>> A = fft(window, 2048) / (len(window)/2.0)
  1356. >>> freq = np.linspace(-0.5, 0.5, len(A))
  1357. >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
  1358. >>> plt.plot(freq, response)
  1359. >>> plt.axis([-0.5, 0.5, -35, 0])
  1360. >>> plt.title("Frequency response of the Exponential window (tau=3.0)")
  1361. >>> plt.ylabel("Normalized magnitude [dB]")
  1362. >>> plt.xlabel("Normalized frequency [cycles per sample]")
  1363. This function can also generate non-symmetric windows:
  1364. >>> tau2 = -(M-1) / np.log(0.01)
  1365. >>> window2 = signal.exponential(M, 0, tau2, False)
  1366. >>> plt.figure()
  1367. >>> plt.plot(window2)
  1368. >>> plt.ylabel("Amplitude")
  1369. >>> plt.xlabel("Sample")
  1370. """
  1371. if sym and center is not None:
  1372. raise ValueError("If sym==True, center must be None.")
  1373. if _len_guards(M):
  1374. return np.ones(M)
  1375. M, needs_trunc = _extend(M, sym)
  1376. if center is None:
  1377. center = (M-1) / 2
  1378. n = np.arange(0, M)
  1379. w = np.exp(-np.abs(n-center) / tau)
  1380. return _truncate(w, needs_trunc)
  1381. def dpss(M, NW, Kmax=None, sym=True, norm=None, return_ratios=False):
  1382. """
  1383. Compute the Discrete Prolate Spheroidal Sequences (DPSS).
  1384. DPSS (or Slepian sequences) are often used in multitaper power spectral
  1385. density estimation (see [1]_). The first window in the sequence can be
  1386. used to maximize the energy concentration in the main lobe, and is also
  1387. called the Slepian window.
  1388. Parameters
  1389. ----------
  1390. M : int
  1391. Window length.
  1392. NW : float
  1393. Standardized half bandwidth corresponding to ``2*NW = BW/f0 = BW*N*dt``
  1394. where ``dt`` is taken as 1.
  1395. Kmax : int | None, optional
  1396. Number of DPSS windows to return (orders ``0`` through ``Kmax-1``).
  1397. If None (default), return only a single window of shape ``(M,)``
  1398. instead of an array of windows of shape ``(Kmax, M)``.
  1399. sym : bool, optional
  1400. When True (default), generates a symmetric window, for use in filter
  1401. design.
  1402. When False, generates a periodic window, for use in spectral analysis.
  1403. norm : {2, 'approximate', 'subsample'} | None, optional
  1404. If 'approximate' or 'subsample', then the windows are normalized by the
  1405. maximum, and a correction scale-factor for even-length windows
  1406. is applied either using ``M**2/(M**2+NW)`` ("approximate") or
  1407. a FFT-based subsample shift ("subsample"), see Notes for details.
  1408. If None, then "approximate" is used when ``Kmax=None`` and 2 otherwise
  1409. (which uses the l2 norm).
  1410. return_ratios : bool, optional
  1411. If True, also return the concentration ratios in addition to the
  1412. windows.
  1413. Returns
  1414. -------
  1415. v : ndarray, shape (Kmax, N) or (N,)
  1416. The DPSS windows. Will be 1D if `Kmax` is None.
  1417. r : ndarray, shape (Kmax,) or float, optional
  1418. The concentration ratios for the windows. Only returned if
  1419. `return_ratios` evaluates to True. Will be 0D if `Kmax` is None.
  1420. Notes
  1421. -----
  1422. This computation uses the tridiagonal eigenvector formulation given
  1423. in [2]_.
  1424. The default normalization for ``Kmax=None``, i.e. window-generation mode,
  1425. simply using the l-infinity norm would create a window with two unity
  1426. values, which creates slight normalization differences between even and odd
  1427. orders. The approximate correction of ``M**2/float(M**2+NW)`` for even
  1428. sample numbers is used to counteract this effect (see Examples below).
  1429. For very long signals (e.g., 1e6 elements), it can be useful to compute
  1430. windows orders of magnitude shorter and use interpolation (e.g.,
  1431. `scipy.interpolate.interp1d`) to obtain tapers of length `M`,
  1432. but this in general will not preserve orthogonality between the tapers.
  1433. .. versionadded:: 1.1
  1434. References
  1435. ----------
  1436. .. [1] Percival DB, Walden WT. Spectral Analysis for Physical Applications:
  1437. Multitaper and Conventional Univariate Techniques.
  1438. Cambridge University Press; 1993.
  1439. .. [2] Slepian, D. Prolate spheroidal wave functions, Fourier analysis, and
  1440. uncertainty V: The discrete case. Bell System Technical Journal,
  1441. Volume 57 (1978), 1371430.
  1442. .. [3] Kaiser, JF, Schafer RW. On the Use of the I0-Sinh Window for
  1443. Spectrum Analysis. IEEE Transactions on Acoustics, Speech and
  1444. Signal Processing. ASSP-28 (1): 105-107; 1980.
  1445. Examples
  1446. --------
  1447. We can compare the window to `kaiser`, which was invented as an alternative
  1448. that was easier to calculate [3]_ (example adapted from
  1449. `here <https://ccrma.stanford.edu/~jos/sasp/Kaiser_DPSS_Windows_Compared.html>`_):
  1450. >>> import numpy as np
  1451. >>> import matplotlib.pyplot as plt
  1452. >>> from scipy.signal import windows, freqz
  1453. >>> N = 51
  1454. >>> fig, axes = plt.subplots(3, 2, figsize=(5, 7))
  1455. >>> for ai, alpha in enumerate((1, 3, 5)):
  1456. ... win_dpss = windows.dpss(N, alpha)
  1457. ... beta = alpha*np.pi
  1458. ... win_kaiser = windows.kaiser(N, beta)
  1459. ... for win, c in ((win_dpss, 'k'), (win_kaiser, 'r')):
  1460. ... win /= win.sum()
  1461. ... axes[ai, 0].plot(win, color=c, lw=1.)
  1462. ... axes[ai, 0].set(xlim=[0, N-1], title=r'$\\alpha$ = %s' % alpha,
  1463. ... ylabel='Amplitude')
  1464. ... w, h = freqz(win)
  1465. ... axes[ai, 1].plot(w, 20 * np.log10(np.abs(h)), color=c, lw=1.)
  1466. ... axes[ai, 1].set(xlim=[0, np.pi],
  1467. ... title=r'$\\beta$ = %0.2f' % beta,
  1468. ... ylabel='Magnitude (dB)')
  1469. >>> for ax in axes.ravel():
  1470. ... ax.grid(True)
  1471. >>> axes[2, 1].legend(['DPSS', 'Kaiser'])
  1472. >>> fig.tight_layout()
  1473. >>> plt.show()
  1474. And here are examples of the first four windows, along with their
  1475. concentration ratios:
  1476. >>> M = 512
  1477. >>> NW = 2.5
  1478. >>> win, eigvals = windows.dpss(M, NW, 4, return_ratios=True)
  1479. >>> fig, ax = plt.subplots(1)
  1480. >>> ax.plot(win.T, linewidth=1.)
  1481. >>> ax.set(xlim=[0, M-1], ylim=[-0.1, 0.1], xlabel='Samples',
  1482. ... title='DPSS, M=%d, NW=%0.1f' % (M, NW))
  1483. >>> ax.legend(['win[%d] (%0.4f)' % (ii, ratio)
  1484. ... for ii, ratio in enumerate(eigvals)])
  1485. >>> fig.tight_layout()
  1486. >>> plt.show()
  1487. Using a standard :math:`l_{\\infty}` norm would produce two unity values
  1488. for even `M`, but only one unity value for odd `M`. This produces uneven
  1489. window power that can be counteracted by the approximate correction
  1490. ``M**2/float(M**2+NW)``, which can be selected by using
  1491. ``norm='approximate'`` (which is the same as ``norm=None`` when
  1492. ``Kmax=None``, as is the case here). Alternatively, the slower
  1493. ``norm='subsample'`` can be used, which uses subsample shifting in the
  1494. frequency domain (FFT) to compute the correction:
  1495. >>> Ms = np.arange(1, 41)
  1496. >>> factors = (50, 20, 10, 5, 2.0001)
  1497. >>> energy = np.empty((3, len(Ms), len(factors)))
  1498. >>> for mi, M in enumerate(Ms):
  1499. ... for fi, factor in enumerate(factors):
  1500. ... NW = M / float(factor)
  1501. ... # Corrected using empirical approximation (default)
  1502. ... win = windows.dpss(M, NW)
  1503. ... energy[0, mi, fi] = np.sum(win ** 2) / np.sqrt(M)
  1504. ... # Corrected using subsample shifting
  1505. ... win = windows.dpss(M, NW, norm='subsample')
  1506. ... energy[1, mi, fi] = np.sum(win ** 2) / np.sqrt(M)
  1507. ... # Uncorrected (using l-infinity norm)
  1508. ... win /= win.max()
  1509. ... energy[2, mi, fi] = np.sum(win ** 2) / np.sqrt(M)
  1510. >>> fig, ax = plt.subplots(1)
  1511. >>> hs = ax.plot(Ms, energy[2], '-o', markersize=4,
  1512. ... markeredgecolor='none')
  1513. >>> leg = [hs[-1]]
  1514. >>> for hi, hh in enumerate(hs):
  1515. ... h1 = ax.plot(Ms, energy[0, :, hi], '-o', markersize=4,
  1516. ... color=hh.get_color(), markeredgecolor='none',
  1517. ... alpha=0.66)
  1518. ... h2 = ax.plot(Ms, energy[1, :, hi], '-o', markersize=4,
  1519. ... color=hh.get_color(), markeredgecolor='none',
  1520. ... alpha=0.33)
  1521. ... if hi == len(hs) - 1:
  1522. ... leg.insert(0, h1[0])
  1523. ... leg.insert(0, h2[0])
  1524. >>> ax.set(xlabel='M (samples)', ylabel=r'Power / $\\sqrt{M}$')
  1525. >>> ax.legend(leg, ['Uncorrected', r'Corrected: $\\frac{M^2}{M^2+NW}$',
  1526. ... 'Corrected (subsample)'])
  1527. >>> fig.tight_layout()
  1528. """ # noqa: E501
  1529. if _len_guards(M):
  1530. return np.ones(M)
  1531. if norm is None:
  1532. norm = 'approximate' if Kmax is None else 2
  1533. known_norms = (2, 'approximate', 'subsample')
  1534. if norm not in known_norms:
  1535. raise ValueError('norm must be one of %s, got %s'
  1536. % (known_norms, norm))
  1537. if Kmax is None:
  1538. singleton = True
  1539. Kmax = 1
  1540. else:
  1541. singleton = False
  1542. Kmax = operator.index(Kmax)
  1543. if not 0 < Kmax <= M:
  1544. raise ValueError('Kmax must be greater than 0 and less than M')
  1545. if NW >= M/2.:
  1546. raise ValueError('NW must be less than M/2.')
  1547. if NW <= 0:
  1548. raise ValueError('NW must be positive')
  1549. M, needs_trunc = _extend(M, sym)
  1550. W = float(NW) / M
  1551. nidx = np.arange(M)
  1552. # Here we want to set up an optimization problem to find a sequence
  1553. # whose energy is maximally concentrated within band [-W,W].
  1554. # Thus, the measure lambda(T,W) is the ratio between the energy within
  1555. # that band, and the total energy. This leads to the eigen-system
  1556. # (A - (l1)I)v = 0, where the eigenvector corresponding to the largest
  1557. # eigenvalue is the sequence with maximally concentrated energy. The
  1558. # collection of eigenvectors of this system are called Slepian
  1559. # sequences, or discrete prolate spheroidal sequences (DPSS). Only the
  1560. # first K, K = 2NW/dt orders of DPSS will exhibit good spectral
  1561. # concentration
  1562. # [see https://en.wikipedia.org/wiki/Spectral_concentration_problem]
  1563. # Here we set up an alternative symmetric tri-diagonal eigenvalue
  1564. # problem such that
  1565. # (B - (l2)I)v = 0, and v are our DPSS (but eigenvalues l2 != l1)
  1566. # the main diagonal = ([N-1-2*t]/2)**2 cos(2PIW), t=[0,1,2,...,N-1]
  1567. # and the first off-diagonal = t(N-t)/2, t=[1,2,...,N-1]
  1568. # [see Percival and Walden, 1993]
  1569. d = ((M - 1 - 2 * nidx) / 2.) ** 2 * np.cos(2 * np.pi * W)
  1570. e = nidx[1:] * (M - nidx[1:]) / 2.
  1571. # only calculate the highest Kmax eigenvalues
  1572. w, windows = linalg.eigh_tridiagonal(
  1573. d, e, select='i', select_range=(M - Kmax, M - 1))
  1574. w = w[::-1]
  1575. windows = windows[:, ::-1].T
  1576. # By convention (Percival and Walden, 1993 pg 379)
  1577. # * symmetric tapers (k=0,2,4,...) should have a positive average.
  1578. fix_even = (windows[::2].sum(axis=1) < 0)
  1579. for i, f in enumerate(fix_even):
  1580. if f:
  1581. windows[2 * i] *= -1
  1582. # * antisymmetric tapers should begin with a positive lobe
  1583. # (this depends on the definition of "lobe", here we'll take the first
  1584. # point above the numerical noise, which should be good enough for
  1585. # sufficiently smooth functions, and more robust than relying on an
  1586. # algorithm that uses max(abs(w)), which is susceptible to numerical
  1587. # noise problems)
  1588. thresh = max(1e-7, 1. / M)
  1589. for i, w in enumerate(windows[1::2]):
  1590. if w[w * w > thresh][0] < 0:
  1591. windows[2 * i + 1] *= -1
  1592. # Now find the eigenvalues of the original spectral concentration problem
  1593. # Use the autocorr sequence technique from Percival and Walden, 1993 pg 390
  1594. if return_ratios:
  1595. dpss_rxx = _fftautocorr(windows)
  1596. r = 4 * W * np.sinc(2 * W * nidx)
  1597. r[0] = 2 * W
  1598. ratios = np.dot(dpss_rxx, r)
  1599. if singleton:
  1600. ratios = ratios[0]
  1601. # Deal with sym and Kmax=None
  1602. if norm != 2:
  1603. windows /= windows.max()
  1604. if M % 2 == 0:
  1605. if norm == 'approximate':
  1606. correction = M**2 / float(M**2 + NW)
  1607. else:
  1608. s = np.fft.rfft(windows[0])
  1609. shift = -(1 - 1./M) * np.arange(1, M//2 + 1)
  1610. s[1:] *= 2 * np.exp(-1j * np.pi * shift)
  1611. correction = M / s.real.sum()
  1612. windows *= correction
  1613. # else we're already l2 normed, so do nothing
  1614. if needs_trunc:
  1615. windows = windows[:, :-1]
  1616. if singleton:
  1617. windows = windows[0]
  1618. return (windows, ratios) if return_ratios else windows
  1619. def _fftautocorr(x):
  1620. """Compute the autocorrelation of a real array and crop the result."""
  1621. N = x.shape[-1]
  1622. use_N = fftpack.next_fast_len(2*N-1)
  1623. x_fft = np.fft.rfft(x, use_N, axis=-1)
  1624. cxy = np.fft.irfft(x_fft * x_fft.conj(), n=use_N)[:, :N]
  1625. # Or equivalently (but in most cases slower):
  1626. # cxy = np.array([np.convolve(xx, yy[::-1], mode='full')
  1627. # for xx, yy in zip(x, x)])[:, N-1:2*N-1]
  1628. return cxy
  1629. _win_equiv_raw = {
  1630. ('barthann', 'brthan', 'bth'): (barthann, False),
  1631. ('bartlett', 'bart', 'brt'): (bartlett, False),
  1632. ('blackman', 'black', 'blk'): (blackman, False),
  1633. ('blackmanharris', 'blackharr', 'bkh'): (blackmanharris, False),
  1634. ('bohman', 'bman', 'bmn'): (bohman, False),
  1635. ('boxcar', 'box', 'ones',
  1636. 'rect', 'rectangular'): (boxcar, False),
  1637. ('chebwin', 'cheb'): (chebwin, True),
  1638. ('cosine', 'halfcosine'): (cosine, False),
  1639. ('exponential', 'poisson'): (exponential, True),
  1640. ('flattop', 'flat', 'flt'): (flattop, False),
  1641. ('gaussian', 'gauss', 'gss'): (gaussian, True),
  1642. ('general gaussian', 'general_gaussian',
  1643. 'general gauss', 'general_gauss', 'ggs'): (general_gaussian, True),
  1644. ('hamming', 'hamm', 'ham'): (hamming, False),
  1645. ('hanning', 'hann', 'han'): (hann, False),
  1646. ('kaiser', 'ksr'): (kaiser, True),
  1647. ('nuttall', 'nutl', 'nut'): (nuttall, False),
  1648. ('parzen', 'parz', 'par'): (parzen, False),
  1649. ('slepian', 'slep', 'optimal', 'dpss', 'dss'): (slepian, True),
  1650. ('triangle', 'triang', 'tri'): (triang, False),
  1651. ('tukey', 'tuk'): (tukey, True),
  1652. }
  1653. # Fill dict with all valid window name strings
  1654. _win_equiv = {}
  1655. for k, v in _win_equiv_raw.items():
  1656. for key in k:
  1657. _win_equiv[key] = v[0]
  1658. # Keep track of which windows need additional parameters
  1659. _needs_param = set()
  1660. for k, v in _win_equiv_raw.items():
  1661. if v[1]:
  1662. _needs_param.update(k)
  1663. def get_window(window, Nx, fftbins=True):
  1664. """
  1665. Return a window.
  1666. Parameters
  1667. ----------
  1668. window : string, float, or tuple
  1669. The type of window to create. See below for more details.
  1670. Nx : int
  1671. The number of samples in the window.
  1672. fftbins : bool, optional
  1673. If True (default), create a "periodic" window, ready to use with
  1674. `ifftshift` and be multiplied by the result of an FFT (see also
  1675. `fftpack.fftfreq`).
  1676. If False, create a "symmetric" window, for use in filter design.
  1677. Returns
  1678. -------
  1679. get_window : ndarray
  1680. Returns a window of length `Nx` and type `window`
  1681. Notes
  1682. -----
  1683. Window types:
  1684. `boxcar`, `triang`, `blackman`, `hamming`, `hann`, `bartlett`,
  1685. `flattop`, `parzen`, `bohman`, `blackmanharris`, `nuttall`,
  1686. `barthann`, `kaiser` (needs beta), `gaussian` (needs standard
  1687. deviation), `general_gaussian` (needs power, width), `slepian`
  1688. (needs width), `dpss` (needs normalized half-bandwidth),
  1689. `chebwin` (needs attenuation), `exponential` (needs decay scale),
  1690. `tukey` (needs taper fraction)
  1691. If the window requires no parameters, then `window` can be a string.
  1692. If the window requires parameters, then `window` must be a tuple
  1693. with the first argument the string name of the window, and the next
  1694. arguments the needed parameters.
  1695. If `window` is a floating point number, it is interpreted as the beta
  1696. parameter of the `kaiser` window.
  1697. Each of the window types listed above is also the name of
  1698. a function that can be called directly to create a window of
  1699. that type.
  1700. Examples
  1701. --------
  1702. >>> from scipy import signal
  1703. >>> signal.get_window('triang', 7)
  1704. array([ 0.125, 0.375, 0.625, 0.875, 0.875, 0.625, 0.375])
  1705. >>> signal.get_window(('kaiser', 4.0), 9)
  1706. array([ 0.08848053, 0.29425961, 0.56437221, 0.82160913, 0.97885093,
  1707. 0.97885093, 0.82160913, 0.56437221, 0.29425961])
  1708. >>> signal.get_window(4.0, 9)
  1709. array([ 0.08848053, 0.29425961, 0.56437221, 0.82160913, 0.97885093,
  1710. 0.97885093, 0.82160913, 0.56437221, 0.29425961])
  1711. """
  1712. sym = not fftbins
  1713. try:
  1714. beta = float(window)
  1715. except (TypeError, ValueError):
  1716. args = ()
  1717. if isinstance(window, tuple):
  1718. winstr = window[0]
  1719. if len(window) > 1:
  1720. args = window[1:]
  1721. elif isinstance(window, string_types):
  1722. if window in _needs_param:
  1723. raise ValueError("The '" + window + "' window needs one or "
  1724. "more parameters -- pass a tuple.")
  1725. else:
  1726. winstr = window
  1727. else:
  1728. raise ValueError("%s as window type is not supported." %
  1729. str(type(window)))
  1730. try:
  1731. winfunc = _win_equiv[winstr]
  1732. except KeyError:
  1733. raise ValueError("Unknown window type.")
  1734. params = (Nx,) + args + (sym,)
  1735. else:
  1736. winfunc = kaiser
  1737. params = (Nx, beta, sym)
  1738. return winfunc(*params)