R/taper.R
In gjmvanboxtel/DSPutils: Some Useful Signal Processing Utilities in R
# taper.R
# Copyright (C) 2019 Geert van Boxtel,
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
# See also: http://www.gnu.org/licenses/gpl-2.0.txt
#
# Version history
# 20190114 GvB Initial setup for package DSPutils
#
#---------------------------------------------------------------------------------------------------------------------
taper <- function (N=512, type="hann", args) {
# N must be a scalar, not be missing and a positive integer
# is.wholenumber in internal.R
if (length(N) != 1 || is.na(N) || !is.wholenumber(N) || N < 2)
stop ('N must be a positive integer > 1')
# match 'type' argument against possible choices
choices <- c("rectangular", "boxchar", "dirichlet", "triangular", "bartlett", "parzen", "welch",
"sine", "cosine", "hann", "hanning", "hamming","blackman",
"nuttall", "blackman-nuttall", "blackman-harris", "SRS-flat-top",
"gaussian", "confined-gaussian", "approximate-confined-gaussian",
"tukey", "planck-taper", "dpss", "slepian",
"kaiser", "dolph-chebyshev", "chebyshev", "exponential", "poisson",
"bartlett-hann", "planck-bessel", "hann-poisson", "lanczos", "sinc")
type <- match.arg(type, choices)
# Note on terminology: N represents the width, in samples, of a discrete-time, symmetrical window function w[n],
# 0 ≤ n ≤ N − 1. When N is an odd number, the non-flat windows have a singular maximum point. When N is even, they have
# a double maximum. It is sometimes useful to express w[n] as a sequence of samples of the lagged version of a zero-phase
# function: w[n]= w0 * (n-((N-1)/2)), 0 <= n <= N-1. This is done in different ways for the various windows (see below)
# 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.
if (type == "rectangular" || type == "boxchar" || type == "dirichlet")
w <- rep(1, N)
# 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 Fejér window. All three definitions converge at large N.
# Here only the first version is implemented. Subtracting ((N-1)/2)) is omitted because it is already present in n
if (type == "triangular") {
n <- seq(-(N - 1), (N - 1), by = 2)
w <- 1 - abs(n/N)
}
# bartlett: see above, L = N-1
if (type == "bartlett") {
n <- seq(-(N - 1), (N - 1), by = 2)
w <- 1 - abs(n/(N-1))
}
# parzen: The Parzen window, also known as the de la Vallée 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
if (type == "parzen") {
n <- seq(-(N-1)/2, (N-1)/2, length.out=N)
n1 <- n[abs(n) <= (N-1)/4]
n2 <- n[n > (N-1)/4]
n3 <- n[n < (-(N-1)/4)]
w1 <- 1 - 6*(abs(n1)/(N/2))**2 + 6*(abs(n1)/(N/2))**3
w2 <- 2 * (1 - abs(n2)/(N/2))**3
w3 <- 2 * (1-abs(n3)/(N/2))**3
w <- c(w3, w1, w2)
}
# 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.
if (type == "welch") {
n <- 0:(N-1)
w <- 1 - ((n-((N-1)/2)) / ((N-1)/2))^2
}
# w[n] = sin(pi*n/(N-1)) = cos((pi*n)/(N-1)-(pi/2))
# The corresponding w0[n] function is a cosine without the (pi/2) phase offset. So the sine window is sometimes
# called cosine window
if (type == "sine" || type == "cosine") {
n <- 0:(N-1)
w <- sin ((n*pi) / (N-1))
}
# 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*pi*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.
if (type == "hann" || type == "hanning") {
n <- 0:(N-1)
w <- 0.5 * (1 - cos ((2*pi*n)/(N-1)))
}
# Hamming (default). 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*pi*n/(N-1)), with α = 0.54, β = 1 − α = 0.46, instead of both constants being equal to 1/2 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.
if (type == "hamming") {
n <- 0:(N-1)
w <- (25/46) - (21/46) * cos ((2*pi*n)/(N-1))
}
# Blackman windows are defined as: w[n] = a0 - a1*cos((2*pi*n)/(N-1)) + a2*cos((4*pi*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 ≈ 0.42659, a1 = 9240/18608 ≈ 0.49656, and a2 = 1430/18608 ≈ 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 additional 'args' argument.
if (type == "blackman") {
n <- 0:(N-1)
a0 <- (7938/18608)
a1 <- (9240/18608)
a2 <- (1430/18608)
if (!missing(args)) {
args <- match.arg (args, c("exact","truncated"))
if (args == "truncated") {
a0 <- 0.42
a1 <- 0.5
a2 <- 0.08
}
}
w <- a0 - a1*cos((2*pi*n)/(N-1)) + a2*cos((4*pi*n)/(N-1))
}
# 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*pi*n)/(N−1)) + a2*cos((4*pi*n)/(N−1)) − a3*cos((6*pi*n)/(N-1)), where a0 = 0.355768, a1 = 0.487396,
# a2 = 0.144232 ; a3 = 0.012604
if (type == "nuttall") {
n <- 0:(N-1)
a0 <- 0.355768
a1 <- 0.487396
a2 <- 0.144232
a3 <- 0.012604
w <- a0 - a1*cos((2*pi*n)/(N-1)) + a2*cos((4*pi*n)/(N-1)) - a3*cos((6*pi*n)/(N-1))
}
# Blackman-Nuttall: w[n] = a0 − a1*cos((2*pi*n)/(N−1)) + a2*cos((4*pi*n)/(N−1)) − a3*cos((6*pi*n)/(N-1)), where
# a0 = 0.3635819 ; a1 = 0.4891775 ; a2 = 0.1365995 ; a3 = 0.0106411
if (type == "blackman-nuttall") {
n <- 0:(N-1)
a0 <- 0.3635819
a1 <- 0.4891775
a2 <- 0.1365995
a3 <- 0.0106411
w <- a0 - a1*cos((2*pi*n)/(N-1)) + a2*cos((4*pi*n)/(N-1)) - a3*cos((6*pi*n)/(N-1))
}
# 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*pi*n)/(N−1)) + a2*cos((4*pi*n)/(N−1)) − a3*cos((6*pi*n)/(N-1)), where
# a0 = 0.35875, a1 = 0.48829, a2 = 0.14128, a3 = 0.01168
if (type == "blackman-harris") {
n <- 0:(N-1)
a0 <- 0.35875
a1 <- 0.48829
a2 <- 0.14128
a3 <- 0.01168
w <- a0 - a1*cos((2*pi*n)/(N-1)) + a2*cos((4*pi*n)/(N-1)) - a3*cos((6*pi*n)/(N-1))
}
# 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*pi*n)/(N−1)) + a2*cos((4*pi*n)/(N−1)) − a3*cos((6*pi*n)/(N-1)) + a4*cos((8*pi*n)/(N-1)), with
# a0 = 1, a1 = 1.93, a2 = 1.29, a3 = 0.388, a4 = 0.028
if (type == "SRS-flat-top") {
n <- 0:(N-1)
a0 <- 1
a1 <- 1.93
a2 <- 1.29
a3 <- 0.388
a4 <- 0.028
w <- a0 - a1*cos((2*pi*n)/(N-1)) + a2*cos((4*pi*n)/(N-1)) - a3*cos((6*pi*n)/(N-1)) + a4*cos((8*pi*n)/(N-1))
w <- w/max(w)
}
# 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) / (sigma * (N-1)/ 2))^2), sigma <= 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 <= 0.5. Default: 0.4
if (type == "gaussian") {
if (missing (args)) args <- 0.4
if (is.na(args) || !is.numeric(args) || length(args) != 1 || abs(args)>0.5) { # use || to stop on first FALSE
stop ("args (sigma) should be a numerical value <= 0.5")
}
n <- 0:(N-1)
w <- exp (-.5 * ((n - (N-1)/2) / (args * (N-1)/ 2))^2)
}
# 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., & Hägele, D. (2014). Discrete-time windows with minimal RMS bandwidth for given RMS temporal width.
# Signal Processing, 102, 240--246.
# Note GvB: this implementation follows matlab code obtained from the second author
# This window is vey suitable for time-frequency analysis
if (type == 'confined-gaussian') {
if (missing (args)) args <- 1
if (is.na(args) || !is.numeric(args) || length(args) != 1 ) { # use || to stop on first FALSE
stop ("args (alpha) should be a numerical value")
}
alpha <- args*(10/N)^4/4
T <- matrix(0,N,N)
P <- matrix(0,N,N)
for (k in 1:N) {
T[k,k] = (k - (N+1)/2)^2
for (l in 1:N) {
if (k != l) {
P[k,l] = 2*(-1)^(k-l)/(k-l)^2
} else {
P[k,l] = pi^2/3;
}
}
}
w <- eigen(P + alpha*T)
w = abs(w$vector[,ncol(w$vectors)]) # eigenvector with lowest eigenvalue is last one in matrix
w <- w/max(w)
}
# 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.14 N.
# Note GvB: code based on matlab script provided by 2nd author
# The additional parameter alpha is passed via 'args'. Default: 2.5 (as in Matlab script)
if (type == 'approximate-confined-gaussian') {
if (missing (args)) args <- 2.5
if (is.na(args) || !is.numeric(args) || length(args) != 1 ) { # use || to stop on first FALSE
stop ("args (alpha) should be a numerical value")
}
#define Gaussian function
G <- function (x) exp(-(1/2)*(args*(x-(N-1)/2)/((N-1)/2))^2)
n <- 0:(N-1)
if (args < 7) {
w <- G(n) - G(-1/2)*(G(n+N) + G(n-N))/(G(-1/2+N)+G(1/2-N))
} else {
w <- G(n)
}
}
# 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(pi*(2n/(α(N − 1)) − 1))], 0 ⩽ n < (α(N − 1 ))/2
# = 1, α(N − 1 ))/2 ⩽ n ⩽ (N − 1)(1-(a/2))
# = 1/2 * [1 + cos(pi*(2n/(α(N − 1)) − (2/α) + 1))], (N − 1)(1-(a/2)) < n ⩽ (N − 1)
# At α = 0 it becomes rectangular, and at α = 1 it becomes a Hann window.
# The value for α is passed through the 'args' parameter
if (type == "tukey") {
if (missing (args)) args <- 0.5
if (is.na(args) || !is.numeric(args) || length(args) != 1 || args < 0) { # use || to stop on first FALSE
stop ("args (alpha) should be a positive numerical value")
}
# check that 0 < α < 1
if (args > 1) args <- 1
else if (args < 0) args <- 0
# generate window
if (args == 0) {
w <- rep (1, N)
} else if (args == 1) {
# Hann window
n <- 0:(N-1)
w <- 0.5 * (1 - cos ((2*pi*n)/(N-1)))
} else {
# cosine-tapered window
n <- seq(0, 1, length.out=N)[1:(N/2)]
w <- (1 + cos(pi*(2*n/args-1)))/2
w[(floor(args*(N-1)/2)+2):length(w)] <- 1
w <- c(w, rep(1, N%%2), rev(w))
}
}
# 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); 0 ⩽ n < ϵ(N − 1)
# w[n] = 1, ϵ(N − 1) < n < (1 − ϵ)(N − 1)
# w[n] = 1 / (exp(Z+)+1); (1 − ϵ)(N − 1) < n ⩽ (N − 1)
# w[n] = 0 otherwise
# where
# Z ± (n;ϵ) = 2ϵ[ (1 / (1 ± ( 2 n / ( N − 1 ) − 1 )) )+ (1 / (1 − 2 ϵ ± ( 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 <= 0.5, as greater values
# result in no flat region.
if (type == "planck-taper") {
if (missing (args)) args <- 0.1
if (is.na(args) || !is.numeric(args) || length(args) != 1 || args < 0 || args > 0.5) { # use || to stop on first FALSE
stop ("args (epsilon) should be a positive numerical value <= 0.5")
}
t2 <- min(N/2, floor(args * (N-1)))
w1 <- seq(1, t2)
w1 <- t2/w1 + t2/(w1-t2-1)
w1 <- 1/(exp(w1)+1)
w2 <- rep(1, N-2 - length(w1)*2)
w <- c(0, w1,w2,rev(w1),0)
}
# 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.
if (type == "dpss" || type == "slepian") {
if (missing (args)) args <- 0.3
if (is.na(args) || !is.numeric(args) || length(args) != 1 || args < 0) { # use || to stop on first FALSE
stop ("args (Wc) should be a positive numerical value.")
if (N * args > 27.38) stop ('Cannot reliably obtain Slepian sequences for N * Wc > 27.38')
}
n <- 1:(N-1)
s <- sin(args*n)/n
c0 <- c(args,s)
A <- stats::toeplitz(c0)
eig <- eigen(A)
w <- eig$vectors[,1]
w <- w/max(w)
}
# The Kaiser, or Kaiser-Bessel, window is a simple approximation of the DPSS window using Bessel functions, discovered by
# James Kaiser.
# w[n] = (I0 (pi*α*sqrt(1 − (2*pi/(N − 1) − 1 )^2 ) ))/(I0 ( pi* α ))
# 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*sqsrt( 1 + α^2), in units of DFT bins, and a typical value of α is 3.
# Sometimes the formula for w[n] is written in terms of a parameter β =def pi*α.
# zero-phase version:
# w0[n] = (I0 pi*α*sqrt(1 − ((2*n)/(N − 1) )^2) )
if (type == "kaiser") {
if (missing (args)) args <- 3
if (is.na(args) || !is.numeric(args) || length(args) != 1 || args < 0) # use || to stop on first FALSE
stop ("args (alpha) should be a positive numerical value.")
n <- 0:(N-1)
w <- besselI (args*pi * sqrt (1 - (2*(n)/(N-1) -1)^2), 0) / besselI(args*pi, 0)
}
# 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((pi*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] ⋅ e^(i*2*pi*k*n)/N, −N/2 ≤ n ≤ N / 2.
# Variations:
# - Due to the equiripple condition, the time-domain window has discontinuities at the edges. An approximation that avoids them,
# by allowing the equiripples to drop off at the edges, is a Taylor window.
# - An alternative to the inverse DFT definition is also available.
# This implementation is copied from R package signal. Argument args = at (default: 100)
if (type == "dolph-chebyshev" || type == "chebyshev") {
if (missing (args)) args <- 100
if (is.na(args) || !is.numeric(args) || length(args) != 1 || args < 0) # use || to stop on first FALSE
stop ("args (at) should be a positive numerical value.")
# auxiliary function to computer Chebyshev polynomial
cheb <- function(n, x) {
if (!(is.numeric(n) && (n == round(n)) && (n >= 0)))
stop("n has to be a positive integer")
T <- numeric(length(x))
ind <- x <= 1
if (any(ind))
T[ind] <- cos(n * acos(as.complex(x[ind])))
ind <- x > 1
myacosh <- function(x) log(x + sqrt(x^2 - 1)) # workaround for a win32 bug in acosh
if (any(ind))
T[ind] <- cosh(n * myacosh(as.complex(x[ind])))
Re(T)
}
# Chebyshev window
gamma <- 10^(-args/20)
beta <- cosh(1/(N - 1) * acosh(1/gamma))
k <- 0:(N - 1)
x <- beta * cos(pi * k/N)
p <- cheb(N - 1, x)
if (N%%2) {
w <- Re(stats::fft(p))
M <- (N + 1)/2
w <- w[1:M]/w[1]
w <- c(w[M:2], w)
} else {
p <- p * exp((0+1i) * pi/N * (0:(N - 1)))
w <- Re(stats::fft(p))
M <- N/2 + 1
w <- w/w[2]
w <- c(w[M:2], w[2:M])
}
}
# The Ultraspherical window's µ parameter determines whether its Fourier transform's side-lobe amplitudes decrease, are level,
# or (shown here) increase with frequency. The Ultraspherical window was introduced in 1984 by Roy Streit and has application
# in antenna array design,[53] non-recursive filter design,[52] and spectrum analysis. Like other adjustable windows,
# the Ultraspherical window has parameters that can be used to control its Fourier transform main-lobe width and relative
# side-lobe amplitude. Uncommon to other windows, it has an additional parameter which can be used to set the rate at which
# side-lobes decrease (or increase) in amplitude.
# The window can be expressed in the time-domain as follows:
# w[n] = 1/N [ C(N −1)~μ (x0) + ∑ (k=1)~(N−1)/2 C (N−1)~μ ( x0 cos((k*pi)/N) ) cos((2*n*pi*k)/N) ]
# where C N~μ is the Ultraspherical polynomial of degree N, and x0 and μ control the side-lobe patterns.
# Certain specific values of μ yield other well-known windows: μ = 0 and μ = 1 give the Dolph–Chebyshev and Saramäki
# windows respectively.
# 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 ≃ 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
if (type == "exponential" || type == "poisson") {
if (missing (args)) args <- N/2
if (is.na(args) || !is.numeric(args) || length(args) != 1 || args < 0) # use || to stop on first FALSE
stop ("args (tau) should be a positive numerical value.")
n <- 0:((N-1)/2)
w <- exp(-( n - (N - 1)/2 )*(1/args))
w <- c(rev(w),w)
w <- w/max(w)
if (N%%2) w <- w[-((N+1)/2)]
}
# Bartlett-Hann
# w[n] = a0 − a1 | n/(N − 1) − 1/2 | − a2 cos ( (2 pi n) / (N − 1) )
# a0 = 0.62, a1 = 0.48, a2 = 0.38
if (type == "bartlett-hann") {
n <- 0:(N-1)
a0 <- 0.62
a1 <- 0.48
a2 <- 0.38
w <- a0 - a1 * abs((n/(N-1)) - 1/2) - a2 * cos ((2*pi*n)/(N-1))
}
# 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.[58] 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)
if (type == "planck-bessel") {
if (missing (args)) args <- c(0.1,3)
if (is.na(args) || !is.numeric(args) || length(args) != 2) { # use || to stop on first FALSE
stop ("args should be a vector of 2 numeric values")
}
eps <- args[1]
alpha <- args[2]
if (eps < 0 || eps > 0.5) stop ('The value for epsilon (args[1]) should be a positive numer <= 0.5')
if (alpha < 0) stop ('The value for alpha (args[2] should be a positive numer')
# planck-taper
t2 <- min(N/2, floor(eps * (N-1)))
w1 <- seq(1, t2)
w1 <- t2/w1 + t2/(w1-t2-1)
w1 <- 1/(exp(w1)+1)
w2 <- rep(1, N-2 - length(w1)*2)
pt <- c(0, w1,w2,rev(w1),0)
# kaiser-bessel
n <- 0:(N-1)
kb <- besselI (alpha*pi * sqrt (1 - (2*(n)/(N-1) -1)^2), 0) / besselI(alpha*pi, 0)
w <- pt * kb
}
# 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*pi*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.
if (type == "hann-poisson") {
if (missing (args)) args <- 3
if (is.na(args) || !is.numeric(args) || length(args) != 1 || args < 0) # use || to stop on first FALSE
stop ("args (alpha) should be a positive numerical value.")
# hann
n <- 0:(N-1)
h <- 0.5 * (1 - cos ((2*pi*n)/(N-1)))
# poisson
n <- 0:((N-1)/2)
p <- exp(-( n - (N - 1)/2 )*(1/args))
p <- c(rev(p),p)
p <- p/max(p)
if (N%%2) p <- p[-((N+1)/2)]
w <- h*p
}
# Lanczos window: 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
if (type == "lanczos" || type =="sinc") {
sinc <- function(x) sin(pi*x)/(pi*x)
n <- 0:(N-1)
w <- sinc(((2*n)/(N-1))-1)
w[is.na(w)] <- 1
}
if (is.null(w) | any(is.na(w)) | length(w) != N)
stop ('Invalid parameters')
return (w)
}
gjmvanboxtel/DSPutils documentation built on May 18, 2019, 2:35 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# taper.R
# Copyright (C) 2019 Geert van Boxtel,
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
# See also: http://www.gnu.org/licenses/gpl-2.0.txt
#
# Version history
# 20190114 GvB Initial setup for package DSPutils
#
#---------------------------------------------------------------------------------------------------------------------
taper <- function (N=512, type="hann", args) {
# N must be a scalar, not be missing and a positive integer
# is.wholenumber in internal.R
if (length(N) != 1 || is.na(N) || !is.wholenumber(N) || N < 2)
stop ('N must be a positive integer > 1')
# match 'type' argument against possible choices
choices <- c("rectangular", "boxchar", "dirichlet", "triangular", "bartlett", "parzen", "welch",
"sine", "cosine", "hann", "hanning", "hamming","blackman",
"nuttall", "blackman-nuttall", "blackman-harris", "SRS-flat-top",
"gaussian", "confined-gaussian", "approximate-confined-gaussian",
"tukey", "planck-taper", "dpss", "slepian",
"kaiser", "dolph-chebyshev", "chebyshev", "exponential", "poisson",
"bartlett-hann", "planck-bessel", "hann-poisson", "lanczos", "sinc")
type <- match.arg(type, choices)
# Note on terminology: N represents the width, in samples, of a discrete-time, symmetrical window function w[n],
# 0 ≤ n ≤ N − 1. When N is an odd number, the non-flat windows have a singular maximum point. When N is even, they have
# a double maximum. It is sometimes useful to express w[n] as a sequence of samples of the lagged version of a zero-phase
# function: w[n]= w0 * (n-((N-1)/2)), 0 <= n <= N-1. This is done in different ways for the various windows (see below)
# 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.
if (type == "rectangular" || type == "boxchar" || type == "dirichlet")
w <- rep(1, N)
# 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 Fejér window. All three definitions converge at large N.
# Here only the first version is implemented. Subtracting ((N-1)/2)) is omitted because it is already present in n
if (type == "triangular") {
n <- seq(-(N - 1), (N - 1), by = 2)
w <- 1 - abs(n/N)
}
# bartlett: see above, L = N-1
if (type == "bartlett") {
n <- seq(-(N - 1), (N - 1), by = 2)
w <- 1 - abs(n/(N-1))
}
# parzen: The Parzen window, also known as the de la Vallée 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
if (type == "parzen") {
n <- seq(-(N-1)/2, (N-1)/2, length.out=N)
n1 <- n[abs(n) <= (N-1)/4]
n2 <- n[n > (N-1)/4]
n3 <- n[n < (-(N-1)/4)]
w1 <- 1 - 6*(abs(n1)/(N/2))**2 + 6*(abs(n1)/(N/2))**3
w2 <- 2 * (1 - abs(n2)/(N/2))**3
w3 <- 2 * (1-abs(n3)/(N/2))**3
w <- c(w3, w1, w2)
}
# 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.
if (type == "welch") {
n <- 0:(N-1)
w <- 1 - ((n-((N-1)/2)) / ((N-1)/2))^2
}
# w[n] = sin(pi*n/(N-1)) = cos((pi*n)/(N-1)-(pi/2))
# The corresponding w0[n] function is a cosine without the (pi/2) phase offset. So the sine window is sometimes
# called cosine window
if (type == "sine" || type == "cosine") {
n <- 0:(N-1)
w <- sin ((n*pi) / (N-1))
}
# 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*pi*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.
if (type == "hann" || type == "hanning") {
n <- 0:(N-1)
w <- 0.5 * (1 - cos ((2*pi*n)/(N-1)))
}
# Hamming (default). 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*pi*n/(N-1)), with α = 0.54, β = 1 − α = 0.46, instead of both constants being equal to 1/2 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.
if (type == "hamming") {
n <- 0:(N-1)
w <- (25/46) - (21/46) * cos ((2*pi*n)/(N-1))
}
# Blackman windows are defined as: w[n] = a0 - a1*cos((2*pi*n)/(N-1)) + a2*cos((4*pi*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 ≈ 0.42659, a1 = 9240/18608 ≈ 0.49656, and a2 = 1430/18608 ≈ 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 additional 'args' argument.
if (type == "blackman") {
n <- 0:(N-1)
a0 <- (7938/18608)
a1 <- (9240/18608)
a2 <- (1430/18608)
if (!missing(args)) {
args <- match.arg (args, c("exact","truncated"))
if (args == "truncated") {
a0 <- 0.42
a1 <- 0.5
a2 <- 0.08
}
}
w <- a0 - a1*cos((2*pi*n)/(N-1)) + a2*cos((4*pi*n)/(N-1))
}
# 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*pi*n)/(N−1)) + a2*cos((4*pi*n)/(N−1)) − a3*cos((6*pi*n)/(N-1)), where a0 = 0.355768, a1 = 0.487396,
# a2 = 0.144232 ; a3 = 0.012604
if (type == "nuttall") {
n <- 0:(N-1)
a0 <- 0.355768
a1 <- 0.487396
a2 <- 0.144232
a3 <- 0.012604
w <- a0 - a1*cos((2*pi*n)/(N-1)) + a2*cos((4*pi*n)/(N-1)) - a3*cos((6*pi*n)/(N-1))
}
# Blackman-Nuttall: w[n] = a0 − a1*cos((2*pi*n)/(N−1)) + a2*cos((4*pi*n)/(N−1)) − a3*cos((6*pi*n)/(N-1)), where
# a0 = 0.3635819 ; a1 = 0.4891775 ; a2 = 0.1365995 ; a3 = 0.0106411
if (type == "blackman-nuttall") {
n <- 0:(N-1)
a0 <- 0.3635819
a1 <- 0.4891775
a2 <- 0.1365995
a3 <- 0.0106411
w <- a0 - a1*cos((2*pi*n)/(N-1)) + a2*cos((4*pi*n)/(N-1)) - a3*cos((6*pi*n)/(N-1))
}
# 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*pi*n)/(N−1)) + a2*cos((4*pi*n)/(N−1)) − a3*cos((6*pi*n)/(N-1)), where
# a0 = 0.35875, a1 = 0.48829, a2 = 0.14128, a3 = 0.01168
if (type == "blackman-harris") {
n <- 0:(N-1)
a0 <- 0.35875
a1 <- 0.48829
a2 <- 0.14128
a3 <- 0.01168
w <- a0 - a1*cos((2*pi*n)/(N-1)) + a2*cos((4*pi*n)/(N-1)) - a3*cos((6*pi*n)/(N-1))
}
# 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*pi*n)/(N−1)) + a2*cos((4*pi*n)/(N−1)) − a3*cos((6*pi*n)/(N-1)) + a4*cos((8*pi*n)/(N-1)), with
# a0 = 1, a1 = 1.93, a2 = 1.29, a3 = 0.388, a4 = 0.028
if (type == "SRS-flat-top") {
n <- 0:(N-1)
a0 <- 1
a1 <- 1.93
a2 <- 1.29
a3 <- 0.388
a4 <- 0.028
w <- a0 - a1*cos((2*pi*n)/(N-1)) + a2*cos((4*pi*n)/(N-1)) - a3*cos((6*pi*n)/(N-1)) + a4*cos((8*pi*n)/(N-1))
w <- w/max(w)
}
# 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) / (sigma * (N-1)/ 2))^2), sigma <= 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 <= 0.5. Default: 0.4
if (type == "gaussian") {
if (missing (args)) args <- 0.4
if (is.na(args) || !is.numeric(args) || length(args) != 1 || abs(args)>0.5) { # use || to stop on first FALSE
stop ("args (sigma) should be a numerical value <= 0.5")
}
n <- 0:(N-1)
w <- exp (-.5 * ((n - (N-1)/2) / (args * (N-1)/ 2))^2)
}
# 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., & Hägele, D. (2014). Discrete-time windows with minimal RMS bandwidth for given RMS temporal width.
# Signal Processing, 102, 240--246.
# Note GvB: this implementation follows matlab code obtained from the second author
# This window is vey suitable for time-frequency analysis
if (type == 'confined-gaussian') {
if (missing (args)) args <- 1
if (is.na(args) || !is.numeric(args) || length(args) != 1 ) { # use || to stop on first FALSE
stop ("args (alpha) should be a numerical value")
}
alpha <- args*(10/N)^4/4
T <- matrix(0,N,N)
P <- matrix(0,N,N)
for (k in 1:N) {
T[k,k] = (k - (N+1)/2)^2
for (l in 1:N) {
if (k != l) {
P[k,l] = 2*(-1)^(k-l)/(k-l)^2
} else {
P[k,l] = pi^2/3;
}
}
}
w <- eigen(P + alpha*T)
w = abs(w$vector[,ncol(w$vectors)]) # eigenvector with lowest eigenvalue is last one in matrix
w <- w/max(w)
}
# 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.14 N.
# Note GvB: code based on matlab script provided by 2nd author
# The additional parameter alpha is passed via 'args'. Default: 2.5 (as in Matlab script)
if (type == 'approximate-confined-gaussian') {
if (missing (args)) args <- 2.5
if (is.na(args) || !is.numeric(args) || length(args) != 1 ) { # use || to stop on first FALSE
stop ("args (alpha) should be a numerical value")
}
#define Gaussian function
G <- function (x) exp(-(1/2)*(args*(x-(N-1)/2)/((N-1)/2))^2)
n <- 0:(N-1)
if (args < 7) {
w <- G(n) - G(-1/2)*(G(n+N) + G(n-N))/(G(-1/2+N)+G(1/2-N))
} else {
w <- G(n)
}
}
# 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(pi*(2n/(α(N − 1)) − 1))], 0 ⩽ n < (α(N − 1 ))/2
# = 1, α(N − 1 ))/2 ⩽ n ⩽ (N − 1)(1-(a/2))
# = 1/2 * [1 + cos(pi*(2n/(α(N − 1)) − (2/α) + 1))], (N − 1)(1-(a/2)) < n ⩽ (N − 1)
# At α = 0 it becomes rectangular, and at α = 1 it becomes a Hann window.
# The value for α is passed through the 'args' parameter
if (type == "tukey") {
if (missing (args)) args <- 0.5
if (is.na(args) || !is.numeric(args) || length(args) != 1 || args < 0) { # use || to stop on first FALSE
stop ("args (alpha) should be a positive numerical value")
}
# check that 0 < α < 1
if (args > 1) args <- 1
else if (args < 0) args <- 0
# generate window
if (args == 0) {
w <- rep (1, N)
} else if (args == 1) {
# Hann window
n <- 0:(N-1)
w <- 0.5 * (1 - cos ((2*pi*n)/(N-1)))
} else {
# cosine-tapered window
n <- seq(0, 1, length.out=N)[1:(N/2)]
w <- (1 + cos(pi*(2*n/args-1)))/2
w[(floor(args*(N-1)/2)+2):length(w)] <- 1
w <- c(w, rep(1, N%%2), rev(w))
}
}
# 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); 0 ⩽ n < ϵ(N − 1)
# w[n] = 1, ϵ(N − 1) < n < (1 − ϵ)(N − 1)
# w[n] = 1 / (exp(Z+)+1); (1 − ϵ)(N − 1) < n ⩽ (N − 1)
# w[n] = 0 otherwise
# where
# Z ± (n;ϵ) = 2ϵ[ (1 / (1 ± ( 2 n / ( N − 1 ) − 1 )) )+ (1 / (1 − 2 ϵ ± ( 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 <= 0.5, as greater values
# result in no flat region.
if (type == "planck-taper") {
if (missing (args)) args <- 0.1
if (is.na(args) || !is.numeric(args) || length(args) != 1 || args < 0 || args > 0.5) { # use || to stop on first FALSE
stop ("args (epsilon) should be a positive numerical value <= 0.5")
}
t2 <- min(N/2, floor(args * (N-1)))
w1 <- seq(1, t2)
w1 <- t2/w1 + t2/(w1-t2-1)
w1 <- 1/(exp(w1)+1)
w2 <- rep(1, N-2 - length(w1)*2)
w <- c(0, w1,w2,rev(w1),0)
}
# 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.
if (type == "dpss" || type == "slepian") {
if (missing (args)) args <- 0.3
if (is.na(args) || !is.numeric(args) || length(args) != 1 || args < 0) { # use || to stop on first FALSE
stop ("args (Wc) should be a positive numerical value.")
if (N * args > 27.38) stop ('Cannot reliably obtain Slepian sequences for N * Wc > 27.38')
}
n <- 1:(N-1)
s <- sin(args*n)/n
c0 <- c(args,s)
A <- stats::toeplitz(c0)
eig <- eigen(A)
w <- eig$vectors[,1]
w <- w/max(w)
}
# The Kaiser, or Kaiser-Bessel, window is a simple approximation of the DPSS window using Bessel functions, discovered by
# James Kaiser.
# w[n] = (I0 (pi*α*sqrt(1 − (2*pi/(N − 1) − 1 )^2 ) ))/(I0 ( pi* α ))
# 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*sqsrt( 1 + α^2), in units of DFT bins, and a typical value of α is 3.
# Sometimes the formula for w[n] is written in terms of a parameter β =def pi*α.
# zero-phase version:
# w0[n] = (I0 pi*α*sqrt(1 − ((2*n)/(N − 1) )^2) )
if (type == "kaiser") {
if (missing (args)) args <- 3
if (is.na(args) || !is.numeric(args) || length(args) != 1 || args < 0) # use || to stop on first FALSE
stop ("args (alpha) should be a positive numerical value.")
n <- 0:(N-1)
w <- besselI (args*pi * sqrt (1 - (2*(n)/(N-1) -1)^2), 0) / besselI(args*pi, 0)
}
# 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((pi*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] ⋅ e^(i*2*pi*k*n)/N, −N/2 ≤ n ≤ N / 2.
# Variations:
# - Due to the equiripple condition, the time-domain window has discontinuities at the edges. An approximation that avoids them,
# by allowing the equiripples to drop off at the edges, is a Taylor window.
# - An alternative to the inverse DFT definition is also available.
# This implementation is copied from R package signal. Argument args = at (default: 100)
if (type == "dolph-chebyshev" || type == "chebyshev") {
if (missing (args)) args <- 100
if (is.na(args) || !is.numeric(args) || length(args) != 1 || args < 0) # use || to stop on first FALSE
stop ("args (at) should be a positive numerical value.")
# auxiliary function to computer Chebyshev polynomial
cheb <- function(n, x) {
if (!(is.numeric(n) && (n == round(n)) && (n >= 0)))
stop("n has to be a positive integer")
T <- numeric(length(x))
ind <- x <= 1
if (any(ind))
T[ind] <- cos(n * acos(as.complex(x[ind])))
ind <- x > 1
myacosh <- function(x) log(x + sqrt(x^2 - 1)) # workaround for a win32 bug in acosh
if (any(ind))
T[ind] <- cosh(n * myacosh(as.complex(x[ind])))
Re(T)
}
# Chebyshev window
gamma <- 10^(-args/20)
beta <- cosh(1/(N - 1) * acosh(1/gamma))
k <- 0:(N - 1)
x <- beta * cos(pi * k/N)
p <- cheb(N - 1, x)
if (N%%2) {
w <- Re(stats::fft(p))
M <- (N + 1)/2
w <- w[1:M]/w[1]
w <- c(w[M:2], w)
} else {
p <- p * exp((0+1i) * pi/N * (0:(N - 1)))
w <- Re(stats::fft(p))
M <- N/2 + 1
w <- w/w[2]
w <- c(w[M:2], w[2:M])
}
}
# The Ultraspherical window's µ parameter determines whether its Fourier transform's side-lobe amplitudes decrease, are level,
# or (shown here) increase with frequency. The Ultraspherical window was introduced in 1984 by Roy Streit and has application
# in antenna array design,[53] non-recursive filter design,[52] and spectrum analysis. Like other adjustable windows,
# the Ultraspherical window has parameters that can be used to control its Fourier transform main-lobe width and relative
# side-lobe amplitude. Uncommon to other windows, it has an additional parameter which can be used to set the rate at which
# side-lobes decrease (or increase) in amplitude.
# The window can be expressed in the time-domain as follows:
# w[n] = 1/N [ C(N −1)~μ (x0) + ∑ (k=1)~(N−1)/2 C (N−1)~μ ( x0 cos((k*pi)/N) ) cos((2*n*pi*k)/N) ]
# where C N~μ is the Ultraspherical polynomial of degree N, and x0 and μ control the side-lobe patterns.
# Certain specific values of μ yield other well-known windows: μ = 0 and μ = 1 give the Dolph–Chebyshev and Saramäki
# windows respectively.
# 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 ≃ 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
if (type == "exponential" || type == "poisson") {
if (missing (args)) args <- N/2
if (is.na(args) || !is.numeric(args) || length(args) != 1 || args < 0) # use || to stop on first FALSE
stop ("args (tau) should be a positive numerical value.")
n <- 0:((N-1)/2)
w <- exp(-( n - (N - 1)/2 )*(1/args))
w <- c(rev(w),w)
w <- w/max(w)
if (N%%2) w <- w[-((N+1)/2)]
}
# Bartlett-Hann
# w[n] = a0 − a1 | n/(N − 1) − 1/2 | − a2 cos ( (2 pi n) / (N − 1) )
# a0 = 0.62, a1 = 0.48, a2 = 0.38
if (type == "bartlett-hann") {
n <- 0:(N-1)
a0 <- 0.62
a1 <- 0.48
a2 <- 0.38
w <- a0 - a1 * abs((n/(N-1)) - 1/2) - a2 * cos ((2*pi*n)/(N-1))
}
# 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.[58] 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)
if (type == "planck-bessel") {
if (missing (args)) args <- c(0.1,3)
if (is.na(args) || !is.numeric(args) || length(args) != 2) { # use || to stop on first FALSE
stop ("args should be a vector of 2 numeric values")
}
eps <- args[1]
alpha <- args[2]
if (eps < 0 || eps > 0.5) stop ('The value for epsilon (args[1]) should be a positive numer <= 0.5')
if (alpha < 0) stop ('The value for alpha (args[2] should be a positive numer')
# planck-taper
t2 <- min(N/2, floor(eps * (N-1)))
w1 <- seq(1, t2)
w1 <- t2/w1 + t2/(w1-t2-1)
w1 <- 1/(exp(w1)+1)
w2 <- rep(1, N-2 - length(w1)*2)
pt <- c(0, w1,w2,rev(w1),0)
# kaiser-bessel
n <- 0:(N-1)
kb <- besselI (alpha*pi * sqrt (1 - (2*(n)/(N-1) -1)^2), 0) / besselI(alpha*pi, 0)
w <- pt * kb
}
# 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*pi*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.
if (type == "hann-poisson") {
if (missing (args)) args <- 3
if (is.na(args) || !is.numeric(args) || length(args) != 1 || args < 0) # use || to stop on first FALSE
stop ("args (alpha) should be a positive numerical value.")
# hann
n <- 0:(N-1)
h <- 0.5 * (1 - cos ((2*pi*n)/(N-1)))
# poisson
n <- 0:((N-1)/2)
p <- exp(-( n - (N - 1)/2 )*(1/args))
p <- c(rev(p),p)
p <- p/max(p)
if (N%%2) p <- p[-((N+1)/2)]
w <- h*p
}
# Lanczos window: 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
if (type == "lanczos" || type =="sinc") {
sinc <- function(x) sin(pi*x)/(pi*x)
n <- 0:(N-1)
w <- sinc(((2*n)/(N-1))-1)
w[is.na(w)] <- 1
}
if (is.null(w) | any(is.na(w)) | length(w) != N)
stop ('Invalid parameters')
return (w)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.