taper: Compute tapering (window) function
In gjmvanboxtel/DSPutils: Some Useful Signal Processing Utilities in R

Description Usage Arguments Details Value Author(s) References Examples

This function will return the results of a tapering function that can be used prior to computing the spectrum of a EEG data channel.

1	taper (N, type="hann", args)

N

Width, in samples, of a discrete-time, symmetrical window function w[n], 0 ≤q n ≤q N - 1. N must be a positive integer. Default: 512.

type

Type of tapering function; must be one of the following:

"rectangular"|"boxchar"|"dirichlet": The rectangular window (sometimes known as the boxcar or Dirichlet window) is the simplest window, equivalent to replacing all but N values of a data sequence by zeros, making it appear as though the waveform suddenly turns on and off: w[n] = 1.
"triangular": Triangular windows are given by: w[n] = 1-abs((n - ((N-1)/2)) / (L/2)), where L can be N, N+1, or N-1. The last one is also known as Bartlett window or Fejer window. All three definitions converge at large N. Specifying "triangular" results in a window with L = N.
"bartlett": Same as "triangular" (see above), but with L = N-1.
"parzen": The Parzen window, also known as the de la Vallee Poussin window, is the 4th order B-spline window given by:
w[n] = 1 - 6*(n/(N-2))^2*(1-(abs(n)/(N/2))), for 0 <= |n| <= (N/4)
w[n] = 2 * (1 - (abs(n)/(N/2)))^3, for N/4 < |n| <= N/2.
"welch": The Welch window consists of a single parabolic section: w[n] = 1 - ((n-(N-1)/2))/((N-1/2))^2. The defining quadratic polynomial reaches a value of zero at the samples just outside the span of the window.
"sine"|"cosine": w[n] = sin(π * n / (N-1)) = cos((π * n) / (N-1) - (pi/2)). The corresponding w0[n] function is a cosine without the (π / 2) phase offset. So the sine window is sometimes called cosine window.
"hann"|"hanning": (Default) The Hann window, named after Julius von Hann, is sometimes referred to as Hanning, presumably due to its linguistic and formulaic similarities to Hamming window. It is also known as raised cosine, because the zero-phase version, w0[n], is one lobe of an elevated cosine function. On the interval n in [0,N-1], the Hann window function is: w[n] = 0.5 * (1 - cos ((2 * π * n)/(N-1))). This function is a member of both the cosine-sum and power-of-sine families. Unlike the Hamming window, the end points of the Hann window just touch zero. The resulting side-lobes roll off at about 18 dB per octave.
"hamming": The window with these particular coefficients was proposed by Richard W. Hamming. The window is optimized to minimize the maximum (nearest) side lobe, giving it a height of about one-fifth that of the Hann window. w[n] = α - β cos(2 * π * n / (N-1)), with α = 0.54, β = 1 - α = 0.46, instead of both constants being equal to 0.5 in the Hann window. The constants are approximations of values α = 25/46 and β = 21/46, which cancel the first sidelobe of the Hann window by placing a zero at frequency 5 π / (N - 1). Approximation of the constants to two decimal places substantially lowers the level of sidelobes, to a nearly equiripple condition. In the equiripple sense, the optimal values for the coefficients are α = 0.53836 and β = 0.46164.
"blackman": Blackman windows are defined as: w[n] = a0 - a1*cos((2 * π * n) / (N-1)) + a2*cos((4 * π *n) / (N-1)), where a0 = ((1 - α)/2), a1 = 1/2, a2 = α/2. By common convention, the unqualified term Blackman window refers to Blackman's "not very serious proposal" of α = 0.16 (a0 = 0.42, a1 = 0.5, a2 = 0.08), which closely approximates the "exact Blackman", with a0 = 7938/18608 \approx 0.42659, a1 = 9240/18608 \approx 0.49656, and a2 = 1430/18608 \approx 0.076849. These exact values place zeros at the third and fourth sidelobes, but result in a discontinuity at the edges and a 6 dB/oct fall-off. The truncated coefficients do not null the sidelobes as well, but have an improved 18 dB/oct fall-off.
Choice of the "exact" (default) or the "truncated" version is done through the optional 'args' argument.
"nuttall": Considering n as a real number, the Nuttall window function and its first derivative are continuous everywhere.That is, the function goes to 0 at n = 0, unlike the Blackman–Nuttall and Blackman–Harris windows, which have a small positive value at zero (at "step" from the zero outside the window), like the Hamming window. The Blackman window defined via α is also continuous with continuous derivative at the edge, but the described "exact Blackman window" is not.
w[n] = a0 - a1*cos((2 * π * n) / (N-1)) + a2*cos((4 * π * n) / (N-1)) - a3*cos((6 * π * n) / (N-1)), where
a0 = 0.355768, a1 = 0.487396, a2 = 0.144232, a3 = 0.012604.
"blackman-nuttal": w[n] = a0 - a1*cos((2 * π * n) / (N-1)) + a2*cos((4 * π * n) / (N-1)) - a3*cos((6 * π * n) / (N-1)), where a0 = 0.3635819, a1 = 0.4891775, a2 = 0.1365995, a3 = 0.0106411.
"blackman-harris": A generalization of the Hamming family, produced by adding more shifted sinc functions, meant to minimize side-lobe levels; w[n] = a0 - a1*cos((2 * π * n) / (N-1)) + a2*cos((4 * π * n) / (N-1)) - a3*cos((6 * π * n) / (N-1)), where a0 = 0.35875, a1 = 0.48829, a2 = 0.14128, a3 = 0.01168.
"SRS-flat-top": A flat top window is a partially negative-valued window that has minimal scalloping loss in the frequency domain. Such windows have been made available in spectrum analyzers for the measurement of amplitudes of sinusoidal frequency components. Drawbacks of the broad bandwidth are poor frequency resolution and high noise bandwidth.
Flat top windows can be designed using low-pass filter design methods, or they may be of the usual cosine-sum variety. An example of the latter is the flat top window available in the Stanford Research Systems (SRS) SR785 spectrum analyzer: w[n] = a0 - a1*cos((2 * π * n) / (N-1)) + a2*cos((4 * π * n) / (N-1)) - a3*cos((6 * π * n) / (N-1)) + a4*cos((8 * π * n) / (N-1)), with a0 = 1, a1 = 1.93, a2 = 1.29, a3 = 0.388, a4 = 0.028.
"gaussian": The Fourier transform of a Gaussian is also a Gaussian (it is an eigenfunction of the Fourier Transform). Since the Gaussian function extends to infinity, it must either be truncated at the ends of the window, or itself windowed with another zero-ended window. Since the log of a Gaussian produces a parabola, this can be used for nearly exact quadratic interpolation in frequency estimation.
exp (-.5 * ((n - (N-1)/2) / (σ * (N-1)/ 2))^2), σ ≤q 0.5
The standard deviation of the Gaussian function is σ (N-1) / 2 sampling periods. Sigma is specified through the additional 'args' argument. Sigma must be ≤q 0.5. Default: 0.4
"confined-gaussian": The confined Gaussian window yields the smallest possible root mean square frequency width σ ω for a given temporal width σ t. These windows optimize the RMS time-frequency bandwidth products. They are computed as the minimum eigenvectors of a parameter-dependent matrix. The confined Gaussian window family contains the cosine window and the Gaussian window in the limiting cases of large and small σ t, respectively.
Starosielec, S., & Haegele, D. (2014). Discrete-time windows with minimal RMS bandwidth for given RMS temporal width. Signal Processing, 102, 240-246.
Note: this implementation follows matlab code obtained from the second author. Default for alpha (via args): 1, which seems to give good results).
This window is vey suitable for time-frequency analysis.
"approximate-confined-gaussian": A confined Gaussian window of temporal width σ t is well approximated by:
w[n] = G[n] - ((G(-1/2)*[G(n + N) + G(n - N) ])/ (G (-1/2 + N) + G (-1/2 - N))) with the Gaussian:
G[x] = exp(- ((x -(N - 1)/2)/(2 * σ t)))^2
The temporal width of the approximate window is asymptotically equal to σ t for σ t < 0.14N.
Note: code based on matlab script provided by 2nd author.
The additional parameter α is passed via 'args'. Default: 2.5 (as in the Matlab script)
"tukey": The Tukey window, also known as the tapered cosine window, can be regarded as a cosine lobe of width α N/2 that is convolved with a rectangular window of width (1 - α / 2) N.
w[n] = 1/2 * [1 + cos( π *(2n/(α (N - 1)) - 1))], for 0 ≤q n < (α (N - 1 ))/2
w[n] = 1, for α (N - 1 ))/2 ≤q n ≤q (N - 1)(1-(a/2))
w[n] = 1/2 * [1 + cos( π * (2n/(α (N - 1)) - (2/α ) + 1))], for (N - 1)(1-(a/2)) < n ≤q (N - 1)
At α = 0 it becomes rectangular, and at α = 1 it becomes a Hann window.
The value for α is passed through the 'args' parameter.
"planck-taper": The so-called "Planck-taper" window is a bump function that has been widely used in the theory of partitions of unity in manifolds. It is smooth (a C^∞ function) everywhere, but is exactly zero outside of a compact region, exactly one over an interval within that region, and varies smoothly and monotonically between those limits. Its use as a window function in signal processing was first suggested in the context of gravitational-wave astronomy, inspired by the Planck distribution. It is defined as a piecewise function:
w[n] = 1 / (exp(Z+)+1) for 0 ≤q n < ε (N - 1)
w[n] = 1 for ε (N - 1) < n < (1 - ε )(N - 1)
w[n] = 1 / (exp(Z+)+1) for (1 - ε )(N - 1) < n ≤q (N - 1)
w[n] = 0 otherwise
where
Z \pm (n;ε) = 2 ε [ (1 / (1 \pm ( 2 n / ( N - 1 ) - 1 )) )+ (1 / (1 - 2 ε \pm ( 2 n / ( N - 1 ) - 1 ) ))].
The amount of tapering (the region over which the function is exactly 1) is controlled by the parameter ε, with smaller values giving sharper transitions (Detault: 0.1, often used is 0.1*fs). The parameter ε should be ≤q 0.5, as greater values result in no flat region.
"dpss"|"slepian": The DPSS (discrete prolate spheroidal sequence) or Slepian window maximizes the energy concentration in the main lobe, and is used in multitaper spectral analysis, which averages out noise in the spectrum and reduces information loss at the edges of the window. The main lobe ends at a frequency bin given by the parameter α.
From https://ccrma.stanford.edu/~jos/sasp/Slepian_DPSS_Window.html: The discrete-time counterpart is Digital Prolate Spheroidal Sequences (DPSS), which may be defined as the eigenvectors of the following symmetric Toeplitz matrix constructed from a sampled sinc function:
S[k,l] = ((sin(omega_cT(k-l))/(k-l)), k,l=0,1,2,...M-1
M denotes the desired window length in samples, omega_c is the desired main-lobe cut-off frequency in radians per second, and T is the sampling period in seconds. The main-lobe bandwidth is thus omega_c rad/sec, counting both positive and negative frequencies.) The digital Slepian window (or DPSS window) is then given by the eigenvector corresponding to the largest eigenvalue.
The omega_c (Wc) parameter is passed through the additional argument args (default: 0.1. N*Wc should be <= 27.38.
"kaiser": The Kaiser, or Kaiser-Bessel, window is a simple approximation of the DPSS window using Bessel functions, discovered by James Kaiser.
w[n] = (I0 (π * α * √{(1 - (2 * π / (N - 1) - 1 )^2 )} ))/(I0 ( π * α ))
where I0 is the zero-th order modified Bessel function of the first kind. Variable parameter α determines the tradeoff between main lobe width and side lobe levels of the spectral leakage pattern. The main lobe width, in between the nulls, is given by 2 * √{(1 + α ^2)}, in units of DFT bins, and a typical value of α is 3.
"dolph-chebyshev"|"chebyshev": Minimizes the Chebyshev norm of the side-lobes for a given main lobe width. The zero-phase Dolph–Chebyshev window function w0[n] is usually defined in terms of its real-valued discrete Fourier transform, W0[k]:
W0[k] = cos{N*cos^-1[β *cos(( π * k) /N)]} / cosh[N * cosh^ - 1 (β)], β = cosh [1 / N cosh ^-1(10 ^ α )], where the parameter α sets the Chebyshev norm of the sidelobes to -20α decibels. The window function can be calculated from W0[k] by an inverse discrete Fourier transform (DFT):
w0[n] = 1/N Σ k = 0, N - 1 W0[k] \cdot e^(i*2* π *k*n)/N, -N/2 ≤q n ≤q N / 2.
This implementation is copied from R package signal. Optional argument args = at (default: 100).
"exponential"|"poisson": The Poisson window, or more generically the exponential window increases exponentially towards the center of the window and decreases exponentially in the second half. Since the exponential function never reaches zero, the values of the window at its limits are non-zero (it can be seen as the multiplication of an exponential function by a rectangular window). It is defined by w[n] = e^(-| n - (N - 1)/2 |(1/τ)), where τ is the time constant of the function. The exponential function decays as e \simeq 2.71828 or approximately 8.69 dB per time constant. This means that for a targeted decay of D dB over half of the window length, the time constant τ is given by τ = (N/2)*(8.69/D).
The time constant is passed as the additional argument. Default: N/2.
"bartlett-hann": w[n] = a0 - a1 | n/(N - 1) - 1/2 | - a2 cos ( (2 π n) / (N - 1) ), a0 = 0.62, a1 = 0.48, a2 = 0.38.
"planck-bessel": A Planck-taper window multiplied by a Kaiser window which is defined in terms of a modified Bessel function. This hybrid window function was introduced to decrease the peak side-lobe level of the Planck-taper window while still exploiting its good asymptotic decay. It has two tunable parameters, ε from the Planck-taper and α from the Kaiser window, so it can be adjusted to fit the requirements of a given signal.
args: c(epsilon, alpha), defaults c(0.1,3).
"hann-poisson": A Hann window multiplied by a Poisson window, which has no side-lobes, in the sense that its Fourier transform drops off forever away from the main lobe. It can thus be used in hill climbing algorithms like Newton's method. The Hann–Poisson window is defined by: w[n] = 1/2 * ( 1 - cos((2* π *n)/(N-1))) * e^((-α|N-1-2n)|/(N-1)) = hav ( (2 π n)/ (N - 1) ) * e^ (-α|N-1-2n|/(N-1)) where α is a parameter that controls the slope of the exponential.
"lanczos"|"sinc": w[n] = sinc(2n/(N-1)-1)
- used in Lanczos resampling
- for the Lanczos window, sinc(x) is defined as sin( π x)/( π x)
- also known as a sinc window, because: w0[n] = sinc(2n/(N-1)) is the main lobe of a normalized sinc function.

args

Further arguments needed for adjustable windows (gaussian, confined-gaussian, approximate-confined-gaussian, generalized-normal, tukey, planck-taper, dpss, slepian, kaiser, dolph-chebyshev, ultraspherical, exponential, poisson, bartlett-hann, planck-bessel, hann-poisson, lanczos)

In signal processing, a window function (also known as an apodization function or tapering function) is a mathematical function that is zero-valued outside of some chosen interval. For instance, a function that is constant inside the interval and zero elsewhere is called a rectangular window, which describes the shape of its graphical representation. When another function or waveform/data-sequence is multiplied by a window function, the product is also zero-valued outside the interval: all that is left is the part where they overlap, the "view through the window".

Tapering functions are applied to time series data such as the EEG before spectral analysis to avoid "spectral leakage", which arises because the FFT can only be performed on limited chunks of data resulting in discontinuities at the edges of those chunks. Multiplying the chunk of data by a window weighting function may taper the edges to (near) zero, thus reducing the discontinuities and the spectral leaking. The "cost" of this procedure is that in the frequency domain the spectral peaks are widened, and the side lobes are not zero. Tapering functions differ in the trade-off between spectral peak width and side lobe attentuation.

Hamming and Hann tapering functions are reasonable choices for many applications, and are often used in EEG analysis. They provide a compromise between the spectral peak spreading across frequency bins and side lobe attenuation. If very low side lobes are needed, a Blackman-Harris window may be required (at the expense of a broader main lobe).

A vector of length N containing the windowing coefficients. When N is an odd number, the non-flat windows have a singular maximum point. When N is even, they have a double maximum.

Geert van Boxtel

The tapering functions, their optional arguments, and their descriptions above are derived from https://en.wikipedia.org/wiki/Window_function

N <- 128
w <- taper (N, "cheby")

## Not run: 
op <- par (mfrow=c(2,1))
plot(w, type="l", main="Chebyshev window", xlab="", ylab="")
plot(10*log10(abs(fft(postpad(w,1024))))[1:512], type="l", 
    main="Frequency response", xlab="", ylab="")
par (op)

## End(Not run)