Package 'gsignal'

Description

Compute the power spectral density of an autoregressive model.

Usage

ar_psd(
  a,
  v = 1,
  freq = 256,
  fs = 1,
  range = ifelse(is.numeric(a), "half", "whole"),
  method = ifelse(length(freq) == 1 && bitwAnd(freq, freq - 1) == 0, "fft", "poly")
)

## S3 method for class 'ar_psd'
plot(
  x,
  yscale = c("linear", "log", "dB"),
  xlab = NULL,
  ylab = NULL,
  main = NULL,
  ...
)

## S3 method for class 'ar_psd'
print(
  x,
  yscale = c("linear", "log", "dB"),
  xlab = NULL,
  ylab = NULL,
  main = NULL,
  ...
)

Arguments

Details

This function calculates the power spectrum of the autoregressive model

                       M
x(n) = sqrt(v).e(n) + SUM a(k).x(n-k)
                      k=1

where x(n) is the output of the model and e(n) is white noise.

Value

An object of class "ar_psd" , which is a list containing two elements, freq and psd containing the frequency values and the estimates of power-spectral density, respectively.

Author(s)

Peter V. Lanspeary, [email protected].
Conversion to R by Geert van Boxtel, [email protected]

Examples

a <- c(1, -2.7607, 3.8106, -2.6535, 0.9238)
psd <- ar_psd(a)

Description

Calculate the coefficients of an autoregressive model using the whitening lattice-filter method of Burg (1968)[1].

Usage

arburg(x, p, criterion = NULL)

Arguments

Details

The inverse of the autoregressive model is a moving-average filter which reduces x to white noise. The power spectrum of the AR model is an estimate of the maximum entropy power spectrum of the data. The function ar_psd calculates the power spectrum of the AR model.

For data input x(n) and white noise e(n), the autoregressive model is

                          p+1
    x(n) = sqrt(v).e(n) + SUM a(k).x(n-k)
                          k=1
 

arburg does not remove the mean from the data. You should remove the mean from the data if you want a power spectrum. A non-zero mean can produce large errors in a power-spectrum estimate. See detrend

Value

A list containing the following elements:

a

vector or matrix containing (p+1) autoregression coefficients. If x is a matrix, then each row of a corresponds to a column of x. a has p + 1 columns.

e

white noise input variance, returned as a vector. If x is a matrix, then each element of e corresponds to a column of x.

k

Reflection coefficients defining the lattice-filter embodiment of the model returned as vector or a matrix. If x is a matrix, then each column of k corresponds to a column of x. k has p rows.

Note

AIC, AICc, KIC and AKICc are based on information theory. They attempt to balance the complexity (or length) of the model against how well the model fits the data. AIC and KIC are biased estimates of the asymmetric and the symmetric Kullback-Leibler divergence, respectively. AICc and AKICc attempt to correct the bias. See reference [2].

Author(s)

Peter V. Lanspeary, [email protected].
Conversion to R by Geert van Boxtel, [email protected].

References

[1] Burg, J.P. (1968) A new analysis technique for time series data, NATO advanced study Institute on Signal Processing with Emphasis on Underwater Acoustics, Enschede, Netherlands, Aug. 12-23, 1968.
[2] Seghouane, A. and Bekara, M. (2004). A small sample model selection criterion based on Kullback’s symmetric divergence. IEEE Trans. Sign. Proc., 52(12), pp 3314-3323,

See Also

ar_psd

Examples

A <- Arma(1, c(1, -2.7607, 3.8106, -2.6535, 0.9238))
y <- filter(A, 0.2 * rnorm(1024))
coefs <- arburg(y, 4)

Description

Create an ARMA model representing a filter or system model, or convert other forms to an ARMA model.

Usage

Arma(b, a)

as.Arma(x, ...)

## S3 method for class 'Arma'
as.Arma(x, ...)

## S3 method for class 'Ma'
as.Arma(x, ...)

## S3 method for class 'Sos'
as.Arma(x, ...)

## S3 method for class 'Zpg'
as.Arma(x, ...)

Arguments

Details

The ARMA model is defined by:

$a(L)y(t) = b(L)x(t)$

The ARMA model can define an analog or digital model. The AR and MA polynomial coefficients follow the convention in 'Matlab' and 'Octave' where the coefficients are in decreasing order of the polynomial (the opposite of the definitions for filterfilter and polyroot). For an analog model,

  H(s) = (b[1]*s^(m-1) + b[2]*s^(m-2) + ... + b[m]) / (a[1]*s^(n-1) +
  a[2]*s^(n-2) + ... + a[n])

For a z-plane digital model,

  H(z) = (b[1] + b[2]*z^(-1) + … + b[m]*z^(-m+1)) / (a[1] + a[2]*z^(-1) + … +
  a[n]*z^(-n+1))

as.Arma converts from other forms, including Zpg and Ma.

Value

A list of class 'Arma' with the following list elements:

b

moving average (MA) polynomial coefficients

a

autoregressive (AR) polynomial coefficients

Author(s)

Tom Short, [email protected],
adapted by Geert van Boxtel, [email protected].

See Also

See also Zpg, Ma, filter, and various filter-generation functions like butter and cheby1 that return Arma models.

Examples

filt <- Arma(b = c(1, 2, 1)/3, a = c(1, 1))
zplane(filt)

Description

compute autoregressive all-pole model parameters using the Yule-Walker method.

Usage

aryule(x, p)

Arguments

Details

aryule uses the Levinson-Durbin recursion on the biased estimate of the sample autocorrelation sequence to compute the parameters.

Value

A list containing the following elements:

a

vector or matrix containing (p + 1) autoregression coefficients. If x is a matrix, then each row of a corresponds to a column of x. a has p + 1 columns.

e

white noise input variance, returned as a vector. If x is a matrix, then each element of e corresponds to a column of x.

k

Reflection coefficients defining the lattice-filter embodiment of the model returned as vector or a matrix. If x is a matrix, then each column of k corresponds to a column of x. k has p rows.

Note

The power spectrum of the resulting filter can be plotted with pyulear(x, p), or you can plot it directly with ar_psd(a,v,...).

Author(s)

Paul Kienzle, [email protected],
Peter V. Lanspeary, [email protected].
Conversion to R by Geert van Boxtel, [email protected].

See Also

ar_psd, arburg

Examples

a <- Arma(1, c(1, -2.7607, 3.8106, -2.6535, 0.9238))
y <- filter(a, rnorm(1024))
coefs <- aryule(y, 4)

Description

Return the filter coefficients of a modified Bartlett-Hann window.

Usage

barthannwin(n)

Arguments

Details

Like Bartlett, Hann, and Hamming windows, the Bartlett-Hann window has a mainlobe at the origin and asymptotically decaying sidelobes on both sides. It is a linear combination of weighted Bartlett and Hann windows with near sidelobes lower than both Bartlett and Hann and with far sidelobes lower than both Bartlett and Hamming windows. The mainlobe width of the modified Bartlett-Hann window is not increased relative to either Bartlett or Hann window mainlobes.

Value

Modified Bartlett-Hann window, returned as a vector. If you specify a one-point window (n = 1), the value 1 is returned.

Author(s)

Andreas Weingessel, [email protected]. Conversion to R by Geert van Boxtel, [email protected].

See Also

bartlett, hann, hamming

Examples

t <- barthannwin(64)
plot (t, type = "l", xlab = "Samples", ylab =" Amplitude")

Description

Return the filter coefficients of a Bartlett (triangular) window.

Usage

bartlett(n)

Arguments

Details

The Bartlett window is very similar to a triangular window as returned by the triang function. However, the Bartlett window always has zeros at the first and last samples, while the triangular window is nonzero at those points.

Value

Bartlett window, returned as a vector. If you specify a one-point window (n = 1), the value 1 is returned.

Author(s)

Andreas Weingessel, [email protected]. Conversion to R by Geert van Boxtel, [email protected].

See Also

triang

Examples

bw <- bartlett(64)
plot (bw, type = "l", xlab = "Samples", ylab =" Amplitude")

Description

Return the poles and gain of a Bessel analog low-pass filter prototype.

Usage

besselap(n)

Arguments

Details

The transfer function is

                     k
 H(s) = -----------------------------
         (s-p(1))(s-p(2))...(s-p(n))

besselap normalizes the poles and gain so that at low frequency and high frequency the Bessel prototype is asymptotically equivalent to the Butterworth prototype of the same order. The magnitude of the filter is less than $1/\sqrt{2}$ at the unity cutoff frequency $\Omega_c = 1$.

Analog Bessel filters are characterized by a group delay that is maximally flat at zero frequency and almost constant throughout the passband. The group delay at zero frequency is

   /  (2n!) \ 2
   | ------ |
   \ 2^n n! /

Value

List of class Zpg containing poles and gain of the filter

Author(s)

Thomas Sailer, [email protected].
Conversion to R by Geert van Boxtel, [email protected].

References

https://en.wikipedia.org/wiki/Bessel_polynomials

Examples

## 6th order Bessel low-pass analog filter
zp <- besselap(6)
w <- seq(0, 4, length.out = 128)
freqs(zp, w)

Description

Compute the transfer function coefficients of an analog Bessel filter.

Usage

besself(n, w, type = c("low", "high", "stop", "pass"))

Arguments

Details

Bessel filters are characterized by an almost constant group delay across the entire passband, thus preserving the wave shape of filtered signals in the passband.

Lowpass Bessel filters have a monotonically decreasing magnitude response, as do lowpass Butterworth filters. Compared to the Butterworth, Chebyshev, and elliptic filters, the Bessel filter has the slowest rolloff and requires the highest order to meet an attenuation specification.

Value

List of class 'Zpg' containing poles and gain of the filter.

Note

As the important characteristic of a Bessel filter is its maximally-flat group delay, and not the amplitude response, it is inappropriate to use the bilinear transform to convert the analog Bessel filter into a digital form (since this preserves the amplitude response but not the group delay) [1].

Author(s)

Paul Kienzle, [email protected],
Doug Stewart, [email protected],
Thomas Sailer, [email protected].
Conversion to R by Geert van Boxtel, [email protected].

References

[1] https://en.wikipedia.org/wiki/Bessel_filter

Examples

w <- seq(0, 4, length.out = 128)

## 5th order Bessel low-pass analog filter
zp <- besself(5, 1.0)
freqs(zp, w)

## 5th order Bessel high-pass analog filter
zp <- besself(5, 1.0, 'high')
freqs(zp, w)

## 5th order Bessel band-pass analog filter
zp <- besself(5, c(1, 2), 'pass')
freqs(zp, w)

## 5th order Bessel band-stop analog filter
zp <- besself(5, c(1, 2), 'stop')
freqs(zp, w)

Description

Transform a s-plane (analog) filter specification into a z-plane (digital) specification.

Usage

bilinear(Sz, ...)

## S3 method for class 'Zpg'
bilinear(Sz, T = 2 * tan(1/2), ...)

## S3 method for class 'Arma'
bilinear(Sz, T = 2 * tan(1/2), ...)

## Default S3 method:
bilinear(Sz, Sp, Sg, T = 2 * tan(1/2), ...)

Arguments

Details

Given a piecewise flat filter design, you can transform it from the s-plane to the z-plane while maintaining the band edges by means of the bilinear transform. This maps the left hand side of the s-plane into the interior of the unit circle. The mapping is highly non-linear, so you must design your filter with band edges in the s-plane positioned at $2/T tan(wT / 2)$ so that they will be positioned at w after the bilinear transform is complete.

The bilinear transform is:

$z = (1 + sT / 2) / (1 - sT / 2)$

$s = (T / 2) (z - 1) / (z + 1)$

Please note that a pole and a zero at the same place exactly cancel. This is significant since the bilinear transform creates numerous extra poles and zeros, most of which cancel. Those which do not cancel have a “fill-in” effect, extending the shorter of the sets to have the same number of as the longer of the sets of poles and zeros (or at least split the difference in the case of the band pass filter). There may be other opportunistic cancellations, but it will not check for them.

Also note that any pole on the unit circle or beyond will result in an unstable filter. Because of cancellation, this will only happen if the number of poles is smaller than the number of zeros. The analytic design methods all yield more poles than zeros, so this will not be a problem.

Value

For the default case or for bilinear.Zpg, an object of class 'Zpg', containing the list elements:

z

complex vector of the zeros of the transformed model

p

complex vector of the poles of the transformed model

g

gain of the transformed model

For bilinear.Arma, an object of class 'Arma', containing the list elements:

b

moving average (MA) polynomial coefficients

a

autoregressive (AR) polynomial coefficients

Author(s)

Paul Kienzle [email protected]. Conversion to R by Tom Short, adapted by Geert van Boxtel [email protected].

References

https://en.wikipedia.org/wiki/Bilinear_transform

Examples

## 6th order Bessel low-pass analog filter
zp <- besselap(6)
w <- seq(0, 4, length.out = 128)
freqs(zp, w)
zzp <- bilinear(zp)
freqz(zzp)

Description

Reorder the elements of the input vector in bit-reversed order.

Usage

bitrevorder(x, index.return = FALSE)

Arguments

Details

This function is equivalent to calling digitrevorder(x, 2), and is useful for prearranging filter coefficients so that bit-reversed ordering does not have to be performed as part of an fft or ifft computation.

Value

The bit-reversed input vector. If index.return = TRUE, then a list containing the bit-reversed input vector (y), and the digit-reversed indices (i).

Author(s)

Mike Miller.
Port to to by Geert van Boxtel, [email protected].

See Also

digitrevorder, fft, ifft

Examples

x <- 0:15
v <- bitrevorder(x)
dec2bin <- function(x, l)
  substr(paste(as.integer(rev(intToBits(x))), collapse = ""),
  32 - l + 1, 32)
x_bin <- sapply(x, dec2bin, 4)
v_bin <- sapply(v, dec2bin, 4)
data.frame(x, x_bin, v, v_bin)

Description

Return the filter coefficients of a Blackman window.

Usage

blackman(n, method = c("symmetric", "periodic"))

Arguments

Details

The Blackman window is a member of the family of cosine sum windows.

Value

Blackman window, returned as a vector.

Author(s)

Andreas Weingessel, [email protected].
Conversion to R by Geert van Boxtel, [email protected].

Examples

h <- blackman(64)
plot (h, type = "l", xlab = "Samples", ylab =" Amplitude")

bs = blackman(64,'symmetric')
bp = blackman(63,'periodic')
plot (bs, type = "l", xlab = "Samples", ylab =" Amplitude")
lines(bp, col="red")

Description

Return the filter coefficients of a minimum four-term Blackman-Harris window.

Usage

blackmanharris(n, method = c("symmetric", "periodic"))

Arguments

Details

The Blackman window is a member of the family of cosine sum windows. It is a generalization of the Hamming family, produced by adding more shifted sinc functions, meant to minimize side-lobe levels.

Value

Blackman-Harris window, returned as a vector.

Author(s)

Sylvain Pelissier, [email protected].
Conversion to R by Geert van Boxtel, [email protected].

Examples

b <- blackmanharris(64)
plot (b, type = "l", xlab = "Samples", ylab =" Amplitude")

bs = blackmanharris(64,'symmetric')
bp = blackmanharris(63,'periodic')
plot (bs, type = "l", xlab = "Samples", ylab =" Amplitude")
lines(bp, col="red")

Description

Return the filter coefficients of a Blackman-Nuttal window.

Usage

blackmannuttall(n, method = c("symmetric", "periodic"))

Arguments

Details

The Blackman-Nuttall window is a member of the family of cosine sum windows.

Value

Blackman-Nuttall window, returned as a vector.

Author(s)

Muthiah Annamalai, [email protected].
Conversion to R by Geert van Boxtel, [email protected].

Examples

b <- blackmannuttall(64)
plot (b, type = "l", xlab = "Samples", ylab =" Amplitude")

bs = blackmannuttall(64,'symmetric')
bp = blackmannuttall(63,'periodic')
plot (bs, type = "l", xlab = "Samples", ylab =" Amplitude")
lines(bp, col="red")

Description

Return the filter coefficients of a Bohman window.

Usage

bohmanwin(n)

Arguments

Details

A Bohman window is the convolution of two half-duration cosine lobes. In the time domain, it is the product of a triangular window and a single cycle of a cosine with a term added to set the first derivative to zero at the boundary.

Value

Bohman window, returned as a vector. If you specify a one-point window (n = 1), the value 1 is returned.

Author(s)

Sylvain Pelissier, [email protected].
Conversion to R by Geert van Boxtel, [email protected].

See Also

triang

Examples

b <- bohmanwin(64)
plot (b, type = "l", xlab = "Samples", ylab =" Amplitude")

Description

Return the filter coefficients of a boxcar (rectangular) window.

Usage

boxcar(n)

Arguments

Details

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. Other windows are designed to moderate these sudden changes, which reduces scalloping loss and improves dynamic range.

Value

rectangular window, returned as a vector.

Author(s)

Paul Kienzle, [email protected].
Conversion to R by Geert van Boxtel, [email protected].

See Also

triang

Examples

b <- boxcar(64)
plot (b, type = "l", xlab = "Samples", ylab =" Amplitude")

Description

Partition a signal vector into nonoverlapping, overlapping, or underlapping data segments.

Usage

buffer(x, n, p = 0, opt, zopt = FALSE)

Arguments

Details

y <- buffer(x, n) partitions a signal vector x of length L into nonoverlapping data segments of length n. Each data segment occupies one column of matrix output y, which has n rows and ceil(L / n) columns. If L is not evenly divisible by n, the last column is zero-padded to length n.

y <- buffer(x, n, p) overlaps or underlaps successive frames in the output matrix by p samples.

• For 0 < p < n (overlap), buffer repeats the final p samples of each segment at the beginning of the following segment. See the example where x = 1:30, n = 7, and an overlap of p = 3. In this case, the first segment starts with p zeros (the default initial condition), and the number of columns in y is ceil(L / (n - p)).

• For p < 0 (underlap), buffer skips p samples between consecutive segments. See the example where x = 1:30, n = 7, and p = -3. The number of columns in y is ceil(L / (n - p)).

In y <- buffer(x, n, p, opt), opt specifies a vector of samples to precede x[1] in an overlapping buffer, or the number of initial samples to skip in an underlapping buffer.

• For 0 < p < n (overlap), opt specifies a vector of length p to insert before x[1] in the buffer. This vector can be considered an initial condition, which is needed when the current buffering operation is one in a sequence of consecutive buffering operations. To maintain the desired segment overlap from one buffer to the next, opt should contain the final p samples of the previous buffer in the sequence. Set opt to "nodelay" to skip the initial condition and begin filling the buffer immediately with x[1]. In this case, L must be length(p) or longer. See the example where x = 1:30, n = 7, p = 3, and opt = "nodelay".

• For p < 0 (underlap), opt is an integer value in the range 0 : -p specifying the number of initial input samples, x[1:opt], to skip before adding samples to the buffer. The first value in the buffer is therefore x[opt + 1].

The opt option is especially useful when the current buffering operation is one in a sequence of consecutive buffering operations. To maintain the desired frame underlap from one buffer to the next, opt should equal the difference between the total number of points to skip between frames (p) and the number of points that were available to be skipped in the previous input to buffer. If the previous input had fewer than p points that could be skipped after filling the final frame of that buffer, the remaining opt points need to be removed from the first frame of the current buffer. See Continuous Buffering for an example of how this works in practice.

buf <- buffer(..., zopt = TRUE) returns the last p samples of a overlapping buffer in output buf$opt. In an underlapping buffer, buf$opt is the difference between the total number of points to skip between frames (-p) and the number of points in x that were available to be skipped after filling the last frame:

• For 0 < p < n (overlap), buf$opt contains the final p samples in the last frame of the buffer. This vector can be used as the initial condition for a subsequent buffering operation in a sequence of consecutive buffering operations. This allows the desired frame overlap to be maintained from one buffer to the next. See Continuous Buffering below.

• For p < 0 (underlap), buf$opt is the difference between the total number of points to skip between frames (-p) and the number of points in x that were available to be skipped after filling the last frame: buf$opt = m*(n-p) + opt - L where opt on the right is the input argument to buffer, and buf$opt on the left is the output argument. Note that for an underlapping buffer output buf$opt is always zero when output buf$z contains data.
  The opt output for an underlapping buffer is especially useful when the current buffering operation is one in a sequence of consecutive buffering operations. The buf$opt output from each buffering operation specifies the number of samples that need to be skipped at the start of the next buffering operation to maintain the desired frame underlap from one buffer to the next. If fewer than p points were available to be skipped after filling the final frame of the current buffer, the remaining opt points need to be removed from the first frame of the next buffer.

In a sequence of buffering operations, the buf$opt output from each operation should be used as the opt input to the subsequent buffering operation. This ensures that the desired frame overlap or underlap is maintained from buffer to buffer, as well as from frame to frame within the same buffer. See Continuous Buffering below for an example of how this works in practice.

Continuous Buffering

In a continuous buffering operation, the vector input to the buffer function represents one frame in a sequence of frames that make up a discrete signal. These signal frames can originate in a frame-based data acquisition process, or within a frame-based algorithm like the FFT.
As an example, you might acquire data from an A/D card in frames of 64 samples. In the simplest case, you could rebuffer the data into frames of 16 samples; buffer with n = 16 creates a buffer of four frames from each 64-element input frame. The result is that the signal of frame size 64 has been converted to a signal of frame size 16; no samples were added or removed.
In the general case where the original signal frame size, L, is not equally divisible by the new frame size, n, the overflow from the last frame needs to be captured and recycled into the following buffer. You can do this by iteratively calling buffer on input x with the zopt parameter set to TRUE. This simply captures any buffer overflow in buf$z, and prepends the data to the subsequent input in the next call to buffer.
Note that continuous buffering cannot be done without the zopt parameter being set to TRUE, because the last frame of y (buf$y in this case) is zero padded, which adds new samples to the signal.
Continuous buffering in the presence of overlap and underlap is handled with the opt parameter, which is used as both an input (opt and output (buf$opt) to buffer. The two examples on this page demonstrate how the opt parameter should be used.

Value

If zopt equals FALSE (the default), this function returns a single numerical array containing the buffered data (y). If zopt equals TRUE, then a list containing 3 variables is returned: y: the buffered data, z: the over or underlap (if any), opt: the over- or underlap that might be used for a future call to buffer to allow continuous buffering.

Author(s)

David Bateman, [email protected].
Conversion to R by Geert van Boxtel, [email protected]

Examples

## Examples without continuous buffering
y <- buffer(1:10, 5)
y <- buffer(1:10, 4)
y <- buffer(1:30, 7, 3)
y <- buffer(1:30, 7, -3)
y <- buffer(1:30, 7, 3, 'nodelay')

## Continuous buffering examples
# with overlap:
data <- buffer(1:1100, 11)
n <- 4
p <- 1
buf <- list(y = NULL, z = NULL, opt = -5)
for (i in 1:ncol(data)) {
  x <- data[,i]
  buf <- buffer(x = c(buf$z,x), n, p, opt=buf$opt, zopt = TRUE)
}
# with underlap:
data <- buffer(1:1100, 11)
n <- 4
p <- -2
buf <- list(y = NULL, z = NULL, opt = 1)
for (i in 1:ncol(data)) {
  x <- data[,i]
  buf <- buffer(x = c(buf$z,x), n, p, opt=buf$opt, zopt = TRUE)
}

Description

Return the poles and gain of an analog Butterworth lowpass filter prototype.

Usage

buttap(n)

Arguments

Details

This function exists for compatibility with 'Matlab' and 'Octave' only, and is equivalent to butter(n, 1, "low", "s").

Value

List of class Zpg containing poles and gain of the filter.

Author(s)

Carne Draug, [email protected].
Conversion to R by Geert van Boxtel, [email protected].

Examples

## 9th order Butterworth low-pass analog filter
zp <- buttap(9)
w <- seq(0, 4, length.out = 128)
freqs(zp, w)

Description

Compute the transfer function coefficients of a Butterworth filter.

Usage

butter(n, ...)

## S3 method for class 'FilterSpecs'
butter(n, ...)

## Default S3 method:
butter(
  n,
  w,
  type = c("low", "high", "stop", "pass"),
  plane = c("z", "s"),
  output = c("Arma", "Zpg", "Sos"),
  ...
)

Arguments

Details

Butterworth filters have a magnitude response that is maximally flat in the passband and monotonic overall. This smoothness comes at the price of decreased rolloff steepness. Elliptic and Chebyshev filters generally provide steeper rolloff for a given filter order.

Because butter is generic, it can be extended to accept other inputs, using buttord to generate filter criteria for example.

Value

Depending on the value of the output parameter, a list of class Arma, Zpg, or Sos containing the filter coefficients

Author(s)

Paul Kienzle, [email protected],
Doug Stewart, [email protected],
Alexander Klein, [email protected],
John W. Eaton.
Conversion to R by Tom Short,
adapted by Geert van Boxtel [email protected].

References

https://en.wikipedia.org/wiki/Butterworth_filter

See Also

Arma, Zpg, Sos, filter, cheby1, ellip, buttord.

Examples

## 50 Hz notch filter
fs <- 256
bf <- butter(4, c(48, 52) / (fs / 2), "stop")
freqz(bf, fs = fs)

## EEG alpha rhythm (8 - 12 Hz) bandpass filter
fs <- 128
fpass <- c(8, 12)
wpass <- fpass / (fs / 2)
but <- butter(5, wpass, "pass")
freqz(but, fs = fs)

## filter to remove vocals from songs, 25 dB attenuation in stop band
## (not optimal with a Butterworth filter)
fs <- 44100
specs <- buttord(230/(fs/2), 450/(fs/2), 1, 25)
bf <- butter(specs)
freqz(bf, fs = fs)
zplane(bf)

Description

Compute the minimum filter order of a Butterworth filter with the desired response characteristics.

Usage

buttord(Wp, Ws, Rp, Rs, plane = c("z", "s"))

Arguments

Details

Deriving the order and cutoff is based on:

       2                  (2N)       (-R / 10)
 |H(W)|  = 1/[1 + (W / Wc)    ] = 10

With some algebra, you can solve simultaneously for Wc and N given Ws, Rs and Wp,Rp. Rounding N to the next greater integer, one can recalculate the allowable range for Wc (filter characteristic touching the pass band edge or the stop band edge).

For other types of filter, before making the above calculation, the requirements must be transformed to lowpass requirements. After the calculation, Wc must be transformed back to the original filter type.

Value

A list of class FilterSpecs with the following list elements:

n

filter order

Wc

cutoff frequency

type

filter type, normally one of "low", "high", "stop", or "pass".

Author(s)

Paul Kienzle,
adapted by Charles Praplan.
Conversion to R by Tom Short,
adapted by Geert van Boxtel, [email protected].

See Also

butter, FilterSpecs

Examples

## low-pass 30 Hz filter
fs <- 128
butspec <- buttord(30/(fs/2), 40/(fs/2), 0.5, 40)
but <- butter(butspec)
freqz(but, fs = fs)

Description

Return the complex cepstrum of the input vector.

Usage

cceps(x)

Arguments

Details

Cepstral analysis is a nonlinear signal processing technique that is applied most commonly in speech and image processing, or as a tool to investigate periodic structures within frequency spectra, for instance resulting from echos/reflections in the signal or to the occurrence of harmonic frequencies (partials, overtones).

The cepstrum is used in many variants. Most important are the power cepstrum, the complex cepstrum, and real cepstrum. The function cceps implements the complex cepstrum by computing the inverse of the log-transformed FFT, i.e.,

$cceps(x) <- ifft(log(fft(x)))$

However, because taking the logarithm of a complex number can lead to unexpected results, the phase of fft(x) needs to be unwrapped before taking the log.

Value

Complex cepstrum, returned as a vector.

Note

This function returns slightly different results in comparison with the 'Matlab' and 'Octave' equivalents. The 'Octave' version does not apply phase unwrapping, but has an optional correction procedure in case of zero phase at $\pi$ radians. The present implementation does apply phase unwrapping so that the correction procedure is unnecessary. The 'Matlab' implementation also applies phase unwrapping, and a circular shift if necessary to avoid zero phase at $\pi$ radians. The circular shift is not done here. In addition, the 'Octave' version shifts the zero frequency to the center of the series, which neither the 'Matlab' nor the present implementation do.

Author(s)

Geert van Boxtel, [email protected].

References

https://en.wikipedia.org/wiki/Cepstrum

See Also

rceps

Examples

## Generate a sine of frequency 45 Hz, sampled at 100 Hz.
fs <- 100
t <- seq(0, 1.27, 1/fs)
s1 <- sin(2 * pi * 45 * t)
## Add an echo with half the amplitude and 0.2 s later.
s2 <- s1 + 0.5 * c(rep(0L, 20), s1[1:108])
## Compute the complex cepstrum of the signal. Notice the echo at 0.2 s.
cep <- cceps(s2)
plot(t, cep, type="l")

Description

Compute the modulo-n circular convolution.

Usage

cconv(a, b, n = length(a) + length(b) - 1)

Arguments

Details

Linear and circular convolution are fundamentally different operations. Linear convolution of an n-point vector x, and an l-point vector y, has length n + l - 1, and can be computed by the function conv, which uses filter. The circular convolution, by contrast, is equal to the inverse discrete Fourier transform (DFT) of the product of the vectors' DFTs.

For the circular convolution of x and y to be equivalent to their linear convolution, the vectors must be padded with zeros to length at least n + l - 1 before taking the DFT. After inverting the product of the DFTs, only the first n + l - 1 elements should be retained.

For long sequences circular convolution may be more efficient than linear convolution. You can also use cconv to compute the circular cross-correlation of two sequences.

Value

Circular convolution of input vectors, returned as a vector.

Author(s)

Leonardo Araujo.
Conversion to R by Geert van Boxtel, [email protected].

See Also

conv, convolve

Examples

a <- c(1, 2, -1, 1)
b <- c(1, 1, 2, 1, 2, 2, 1, 1)
c <- cconv(a, b)       # Circular convolution
cref = conv(a, b)      # Linear convolution
all.equal(max(c - cref), 0)

cconv(a, b, 6)

Description

Return the value of the Chebyshev polynomial at specific points.

Usage

cheb(n, x)

Arguments

Details

The Chebyshev polynomials are defined by the equations:

  Tn(x) = cos(n . acos(x),    |x|<= 1
  Tn(x) = cosh(n . acosh(x),  |x|> 1

If x is a vector, the output is a vector of the same size, where each element is calculated as $y(i) = Tn(x(i))$.

Value

Polynomial of order x, evaluated at point(s) x.

Author(s)

André Carezia, [email protected].
Conversion to R by Geert van Boxtel, [email protected].

Examples

cp <- cheb(5, 1)
cp <- cheb(5, c(2,3))

Description

Return the poles and gain of an analog Chebyshev Type I lowpass filter prototype.

Usage

cheb1ap(n, Rp)

Arguments

Details

This function exists for compatibility with 'Matlab' and 'OCtave' only, and is equivalent to cheby1(n, Rp, 1, "low", "s").

Value

List of class Zpg containing the poles and gain of the filter.

Author(s)

Carne Draug, [email protected].
Conversion to R by Geert van Boxtel, [email protected].

Examples

## 9th order Chebyshev type I low-pass analog filter
zp <- cheb1ap(9, .1)
w <- seq(0, 4, length.out = 128)
freqs(zp, w)

Description

Compute Chebyshev type-I filter order and cutoff for the desired response characteristics.

Usage

cheb1ord(Wp, Ws, Rp, Rs, plane = c("z", "s"))

Arguments

Value

A list of class 'FilterSpecs' with the following list elements:

n

filter order

Wc

cutoff frequency

type

filter type, normally one of "low", "high", "stop", or "pass".

Author(s)

Paul Kienzle, Laurent S. Mazet, Charles Praplan.
Conversion to R by Tom Short, adapted by Geert van Boxtel, [email protected].

See Also

cheby1

Examples

## low-pass 30 Hz filter
fs <- 128
spec <- cheb1ord(30/(fs/2), 40/(fs/2), 0.5, 40)
cf <- cheby1(spec)
freqz(cf, fs = fs)

Description

Return the poles and gain of an analog Chebyshev Type II lowpass filter prototype.

Usage

cheb2ap(n, Rs)

Arguments

Details

This function exists for compatibility with 'Matlab' and 'Octave' only, and is equivalent to cheby2(n, Rp, 1, "low", "s").

Value

list of class Zpg containing poles and gain of the filter

Author(s)

Carne Draug, [email protected].
Conversion to R by Geert van Boxtel, [email protected].

Examples

## 9th order Chebyshev type II low-pass analog filter
zp <- cheb2ap(9, 30)
w <- seq(0, 4, length.out = 128)
freqs(zp, w)

Description

Compute Chebyshev type-II filter order and cutoff for the desired response characteristics.

Usage

cheb2ord(Wp, Ws, Rp, Rs, plane = c("z", "s"))

Arguments

Value

A list of class 'FilterSpecs' with the following list elements:

n

filter order

Wc

cutoff frequency

type

filter type, normally one of "low", "high", "stop", or "pass".

Author(s)

Paul Kienzle, Charles Praplan.
Conversion to R by Geert van Boxtel, [email protected].

See Also

cheby1

Examples

## low-pass 30 Hz filter
fs <- 128
spec <- cheb2ord(30/(fs/2), 40/(fs/2), 0.5, 40)
cf <- cheby2(spec)
freqz(cf, fs = fs)

Description

Return the filter coefficients of a Dolph-Chebyshev window.

Usage

chebwin(n, at = 100)

Arguments

Details

The window is described in frequency domain by the expression:

                Cheb(m - 1, Beta * cos(\pi * k / m))
         W(k) = ------------------------------------
                       Cheb(m - 1, Beta)

with

  Beta = cosh(1 / (m - 1) * acosh(10^(at / 20))

and and $Cheb(m, x)$ denoting the $m$-th order Chebyshev polynomial calculated at the point $x$.

Note that the denominator in W(k) above is not computed, and after the inverse Fourier transform the window is scaled by making its maximum value unitary.

Value

Chebyshev window, returned as a vector. If you specify a one-point window (n = 1), the value 1 is returned.

Author(s)

André Carezia, [email protected].
Conversion to R by Geert van Boxtel, [email protected].

Examples

cw <- chebwin(64)
plot (cw, type = "l", xlab = "Samples", ylab =" Amplitude")

Description

Compute the transfer function coefficients of a Chebyshev Type I filter.

Usage

cheby1(n, ...)

## S3 method for class 'FilterSpecs'
cheby1(n, ...)

## Default S3 method:
cheby1(
  n,
  Rp,
  w,
  type = c("low", "high", "stop", "pass"),
  plane = c("z", "s"),
  output = c("Arma", "Zpg", "Sos"),
  ...
)

Arguments

Details

Chebyshev filters are analog or digital filters having a steeper roll-off than Butterworth filters, and have passband ripple (type I) or stopband ripple (type II).

Because cheby1 is generic, it can be extended to accept other inputs, using cheb1ord to generate filter criteria for example.

Value

Depending on the value of the output parameter, a list of class Arma, Zpg, or Sos containing the filter coefficients

Author(s)

Paul Kienzle, [email protected],
Doug Stewart, [email protected].
Conversion to R Tom Short,
adapted by Geert van Boxtel, [email protected].

References

https://en.wikipedia.org/wiki/Chebyshev_filter

See Also

Arma, filter, butter, ellip, cheb1ord, FilterSpecs

Examples

## compare the frequency responses of 5th-order
## Butterworth and Chebyshev filters.
bf <- butter(5, 0.1)
cf <- cheby1(5, 3, 0.1)
bfr <- freqz(bf)
cfr <- freqz(cf)
plot(bfr$w / pi, 20 * log10(abs(bfr$h)), type = "l", ylim = c(-40, 0),
  xlim = c(0, .5), xlab = "Frequency", ylab = c("dB"))
lines(cfr$w / pi, 20 * log10(abs(cfr$h)), col = "red")

# compare type I and type II Chebyshev filters.
c1fr <- freqz(cheby1(5, .5, 0.5))
c2fr <- freqz(cheby2(5, 20, 0.5))
plot(c1fr$w / pi, abs(c1fr$h), type = "l", ylim = c(0, 1),
  xlab = "Frequency", ylab = c("Magnitude"))
  lines(c2fr$w / pi, abs(c2fr$h), col = "red")

Description

Compute the transfer function coefficients of a Chebyshev Type II filter.

Usage

cheby2(n, ...)

## S3 method for class 'FilterSpecs'
cheby2(n, ...)

## Default S3 method:
cheby2(
  n,
  Rs,
  w,
  type = c("low", "high", "stop", "pass"),
  plane = c("z", "s"),
  output = c("Arma", "Zpg", "Sos"),
  ...
)

Arguments

Details

Chebyshev filters are analog or digital filters having a steeper roll-off than Butterworth filters, and have passband ripple (type I) or stopband ripple (type II).

Because cheby2 is generic, it can be extended to accept other inputs, using cheb2ord to generate filter criteria for example.

Value

Depending on the value of the output parameter, a list of class Arma, Zpg, or Sos containing the filter coefficients

Author(s)

Paul Kienzle, [email protected],
Doug Stewart, [email protected].
Conversion to R Tom Short,
adapted by Geert van Boxtel, [email protected].

References

https://en.wikipedia.org/wiki/Chebyshev_filter

See Also

Arma, filter, butter, ellip, cheb2ord

Examples

## compare the frequency responses of 5th-order
## Butterworth and Chebyshev filters.
bf <- butter(5, 0.1)
cf <- cheby2(5, 20, 0.1)
bfr <- freqz(bf)
cfr <- freqz(cf)
plot(bfr$w / pi, 20 * log10(abs(bfr$h)), type = "l", ylim = c(-40, 0),
  xlim = c(0, .5), xlab = "Frequency", ylab = c("dB"))
lines(cfr$w / pi, 20 * log10(abs(cfr$h)), col = "red")

# compare type I and type II Chebyshev filters.
c1fr <- freqz(cheby1(5, .5, 0.5))
c2fr <- freqz(cheby2(5, 20, 0.5))
plot(c1fr$w / pi, abs(c1fr$h), type = "l", ylim = c(0, 1.1),
  xlab = "Frequency", ylab = c("Magnitude"))
lines(c2fr$w / pi, abs(c2fr$h), col = "red")

Description

Evaluate a chirp signal (frequency swept cosine wave).

Usage

chirp(
  t,
  f0,
  t1 = 1,
  f1 = 100,
  shape = c("linear", "quadratic", "logarithmic"),
  phase = 0
)

Arguments

Details

A chirp is a signal in which the frequency changes with time, commonly used in sonar, radar, and laser. The name is a reference to the chirping sound made by birds.

The chirp can have one of three shapes:

"linear"

Specifies an instantaneous frequency sweep $f_i(t)$ given by $f_i(t) = f_0 + \beta t$, where $\beta = (f_1 - f_0) / t_1$ and the default value for $f_0$ is 0. The coefficient $\beta$ ensures that the desired frequency breakpoint $f_1$ at time $t_1$ is maintained.

"quadratic"

Specifies an instantaneous frequency sweep $f_i(t)$ given by $f_i(t) = f_0 + \beta t^2$, where $\beta = (f_1 - f_0) / t_1^2$ and the default value for $f_0$ is 0. If $f_0 > f_1$ (downsweep), the default shape is convex. If $f_0 < f_1$ (upsweep), the default shape is concave.

"logarithmic"

Specifies an instantaneous frequency sweep $f_i(t)$ given by $f_i(t) = f_0 \times \beta t$, where $\beta = \left( \frac {f_1}{f_0} \right) ^ \frac{1}{t1}$ and the default value for $f_0$ is $10^{-6}$.

Value

Chirp signal, returned as an array of the same length as t.

Author(s)

Paul Kienzle, [email protected],
Mike Miller.
Conversion to R by Geert van Boxtel, [email protected].

Examples

# Shows linear sweep of 100 Hz/sec starting at zero for 5 sec
# since the sample rate is 1000 Hz, this should be a diagonal
# from bottom left to top right.
t <- seq(0, 5, 0.001)
y <- chirp (t)
specgram (y, 256, 1000)

# Shows a quadratic chirp of 400 Hz at t=0 and 100 Hz at t=10
# Time goes from -2 to 15 seconds.
specgram(chirp(seq(-2, 15, by = 0.001), 400, 10, 100, "quadratic"))

# Shows a logarithmic chirp of 200 Hz at t = 0 and 500 Hz at t = 2
# Time goes from 0 to 5 seconds at 8000 Hz.
specgram(chirp(seq(0, 5, by = 1/8000), 200, 2, 500, "logarithmic"),
         fs = 8000)

Description

Constrained least square band-pass FIR filter design without specified transition bands.

Usage

cl2bp(m = 30, w1, w2, up, lo, L = 2048)

Arguments

Details

This is a fast implementation of the algorithm cited below. Compared to remez, it offers implicit specification of transition bands, a higher likelihood of convergence, and an error criterion combining features of both L2 and Chebyshev approaches

Value

The FIR filter coefficients, a vector of length 2 * m + 1, of class Ma.

Author(s)

Ivan Selesnick, Rice University, 1995, downloaded from https://www.ece.rice.edu/dsp/software/cl2.shtml.
Conversion to R by Geert van Boxtel, [email protected].

References

Selesnick, I.W., Lang, M., and Burrus, C.S. (1998) A modified algorithm for constrained least square design of multiband FIR filters without specified transition bands. IEEE Trans. on Signal Processing, 46(2), 497-501.
https://www.ece.rice.edu/dsp/software/cl2.shtml

See Also

Ma, filter, remez

Examples

w1 <- 0.3 * pi
w2 <- 0.6 * pi
up <- c(0.02, 1.02, 0.02)
lo <- c(-0.02, 0.98, -0.02)
h  <- cl2bp(30, w1, w2, up, lo, 2^11)
freqz(h)

Description

Calculate boundary indexes of clusters of 1’s.

Usage

clustersegment(x)

Arguments

Details

The function calculates the initial index and end index of sequences of 1s rising and falling phases of the signal in x. The clusters are sought in the rows of the array x. The function works by finding the indexes of jumps between consecutive values in the rows of x.

Value

A list of size nr, where nr is the number of rows in x. Each element of the list contains a matrix with two rows. The first row is the initial index of a sequence of 1’s and the second row is the end index of that sequence.

Author(s)

Juan Pablo Carbajal, [email protected].
Conversion to R by Geert van Boxtel, [email protected].

Examples

(x <- c(0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1))
(ranges <- clustersegment(x))
# The first sequence of 1's in x lies in the interval
(r <- ranges[1,1]:ranges[2,1])

x <- matrix(as.numeric(runif(30) > 0.4), 3, 10)
ranges <- clustersegment(x)

x <- c(0, 1.2, 3, -8, 0)
ranges <- clustersegment(x)

Description

Compute the complex Morlet wavelet on a grid.

Usage

cmorwavf(lb = -8, ub = 8, n = 1000, fb = 5, fc = 1)

Arguments

Details

The Morlet (or Gabor) wavelet is a wavelet composed of a complex exponential (carrier) multiplied by a Gaussian window (envelope). The wavelet exists as a complex version or a purely real-valued version. Some distinguish between the "real Morlet" versus the "complex Morlet". Others consider the complex version to be the "Gabor wavelet", while the real-valued version is the "Morlet wavelet". This function returns the complex Morlet wavelet, with time-decay parameter fb, and center frequency fc. The general expression for the complex Morlet wavelet is

 Psi(x) = ((pi * fb)^-0.5) *  exp(2 * pi * i * fc * x) * exp(-(x^2) / fb)

x is evaluated on an n-point regular grid in the interval (lb, ub).

fb controls the decay in the time domain and the corresponding energy spread (bandwidth) in the frequency domain. fb is the inverse of the variance in the frequency domain. Increasing fb makes the wavelet energy more concentrated around the center frequency and results in slower decay of the wavelet in the time domain. Decreasing fb results in faster decay of the wavelet in the time domain and less energy spread in the frequency domain. The value of fb does not affect the center frequency. When converting from scale to frequency, only the center frequency affects the frequency values. The energy spread or bandwidth parameter affects how localized the wavelet is in the frequency domain. See the examples.

Value

A list containing 2 variables; x, the grid on which the complex Morlet wavelet was evaluated, and psi ($\Psi$), the evaluated wavelet on the grid x.

Author(s)

Sylvain Pelissier, [email protected].
Conversion to R by Geert van Boxtel, [email protected].

Examples

## Construct a complex-valued Morlet wavelet with a bandwidth parameter
##  of 1.5 and a center frequency of 1. Set the effective support to [-8,8]
## and the length of the wavelet to 1000.
cmw <- cmorwavf(-8, 8, 1000, 1.5, 1)

# Plot the real and imaginary parts of the wavelet.
op <- par(mfrow = c(2, 1))
plot(cmw$x, Re(cmw$psi), type = "l", main = "Real Part")
plot(cmw$x, Im(cmw$psi), type = "l", main = "Imaginary Part")
par(op)

## This example shows how the complex Morlet wavelet shape in the frequency
## domain is affected by the value of the bandwidth parameter (fb). Both
## wavelets have a center frequency of 1. One wavelet has an fb value of
## 0.5 and the other wavelet has a value of 8.

op <- par(mfrow = c(2,1))
cmw1 <- cmorwavf(fb = 0.5)
cmw2 <- cmorwavf(fb = 8)

# time domain plot
plot(cmw1$x, Re(cmw1$psi), type = "l", xlab = "Time", ylab = "",
     main = "Time domain, Real part")
lines(cmw2$x, Re(cmw2$psi), col = "red")
legend("topright", legend = c("fb = 0.5", "fb = 8"), lty= 1, col = c(1,2))

# frequency domain plot
f <- seq(-5, 5, .01)
Fc <- 1
Fb1 <- 0.5
Fb2 <- 8
PSI1 <- exp(-pi^2 * Fb1 * (f-Fc)^2)
PSI2 <- exp(-pi^2 * Fb2 * (f-Fc)^2)
plot(f, PSI1, type="l", xlab = "Frequency", ylab = "",
     main = "Frequency domain")
lines(f, PSI2, col = "red")
legend("topright", legend = c("fb = 0.5", "fb = 8"),
       lty= 1, col = c(1,2))
par(op)

## The fb bandwidth parameter for the complex Morlet wavelet is the
## inverse of the variance in frequency. Therefore, increasing Fb results
## in a narrower concentration of energy around the center frequency.

## alternative to the above frequency plot:
fs <- length(cmw1$x) / sum(abs(range(cmw1$x)))
hz  <- seq(0, fs/2, len=floor(length(cmw1$psi)/2)+1)
PSI1 <- fft(cmw1$psi) / length(cmw1$psi)
PSI2 <- fft(cmw2$psi) / length(cmw2$psi)
plot(hz, 2 * abs(PSI1)[1:length(hz)], type="l", xlab = "Frequency",
     ylab = "", main = "Frequency domain", xlim=c(0,5))
lines(hz, 2 * abs(PSI2)[1:length(hz)], col = 2)
legend("topright", legend = c("fb = 0.5", "fb = 8"), lty= 1, col = c(1,2))

Description

Convolve two vectors a and b.

Usage

conv(a, b, shape = c("full", "same", "valid"))

Arguments

Details

The convolution of two vectors, a and b, represents the area of overlap under the points as B slides across a. Algebraically, convolution is the same operation as multiplying polynomials whose coefficients are the elements of a and b.

The function conv uses the filter function, NOT fft, which may be faster for large vectors.

Value

Output vector with length equal to length (a) + length (b) - 1. When the parameter shape is set to "valid", the length of the output is max(length(a) - length(b) + 1, 0), except when length(b) is zero. In that case, the length of the output vector equals length(a).

When a and b are the coefficient vectors of two polynomials, the convolution represents the coefficient vector of the product polynomial.

Author(s)

Tony Richardson, [email protected], adapted by John W. Eaton.
Conversion to R by Geert van Boxtel, [email protected].

Examples

u <- rep(1L, 3)
v <- c(1, 1, 0, 0, 0, 1, 1)
w <- conv(u, v)

## Create vectors u and v containing the coefficients of the polynomials
## x^2 + 1 and 2x + 7.
u <- c(1, 0, 1)
v <- c(2, 7)
## Use convolution to multiply the polynomials.
w <- conv(u, v)
## w contains the polynomial coefficients for 2x^3 + 7x^2 + 2x + 7.

## Central part of convolution
u <- c(-1, 2, 3, -2, 0, 1, 2)
v <- c(2, 4, -1, 1)
w <- conv(u, v, 'same')

Description

Compute the two-dimensional convolution of two matrices.

Usage

conv2(a, b, shape = c("full", "same", "valid"))

Arguments

Value

Convolution of input matrices, returned as a matrix.

Author(s)

Geert van Boxtel, [email protected].

See Also

conv, convolve

Examples

a <- matrix(1:16, 4, 4)
b <- matrix(1:9, 3,3)
cnv <- conv2(a, b)
cnv <- conv2(a, b, "same")
cnv <- conv2(a, b, "valid")

Description

Returns the convolution matrix for a filter kernel.

Usage

convmtx(h, n)

Arguments

Details

Computing a convolution using conv when the signals are vectors is generally more efficient than using convmtx. For multichannel signals, however, when a large number of vectors are to be convolved with the same filter kernel, convmtx might be more efficient.

The code cm <- convmtx(h, n) computes the convolution matrix of the filter kernel h with a vector of length n. Then, cm x gives the convolution of h and x.

Value

Convolution matrix of input h for a vector of length n. If h is a vector of length m, then the convolution matrix has m + n - 1 rows and n columns.

Author(s)

David Bateman [email protected].
Conversion to R by Geert van Boxtel [email protected].

See Also

conv

Examples

N <- 1000
a <- runif(N)
b <- runif(N)
cm <- convmtx(b, N)
d <- cm %*% a

cref = conv(a, b)
all.equal(max(d - cref), 0)

Description

Sort complex numbers into complex conjugate pairs ordered by increasing real part.

Usage

cplxpair(z, tol = 100 * .Machine$double.eps, MARGIN = 2)

Arguments

Details

The negative imaginary complex numbers are placed first within each pair. All real numbers (those with abs(Im (z) / z) < tol) are placed after the complex pairs.

An error is signaled if some complex numbers could not be paired and if all complex numbers are not exact conjugates (to within tol).

Value

Vector, matrix or array containing ordered complex conjugate pairs by increasing real parts.

Note

There is no defined order for pairs with identical real parts but differing imaginary parts.

Author(s)

Geert van Boxtel, [email protected].

See Also

cplxreal

Examples

r <- rbind(t(cplxpair(exp(2i * pi * 0:4 / 5))),
           t(exp(2i * pi *c(3, 2, 4, 1, 0) / 5)))

Description

Sort numbers into into complex-conjugate-valued and real-valued elements.

Usage

cplxreal(z, tol = 100 * .Machine$double.eps, MARGIN = 2)

Arguments

Details

An error is signaled if some complex numbers could not be paired and if all complex numbers are not exact conjugates (to within tol). Note that here is no defined order for pairs with identical real parts but differing imaginary parts.

Value

A list containing two variables:

zc

Vector, matrix or array containing ordered complex conjugate pairs by increasing real parts. Only the positive imaginary complex numbers of each complex conjugate pair are returned.

zr

Vector, matrix or array containing ordered real numbers.

Author(s)

Geert van Boxtel, [email protected].

See Also

cplxpair

Examples

r <- cplxreal(c(1, 1 + 3i, 2 - 5i, 1-3i, 2 + 5i, 4, 3))

Description

Estimates the cross power spectral density (CPSD) of discrete-time signals.

Usage

cpsd(
  x,
  window = nextpow2(sqrt(NROW(x))),
  overlap = 0.5,
  nfft = ifelse(isScalar(window), window, length(window)),
  fs = 1,
  detrend = c("long-mean", "short-mean", "long-linear", "short-linear", "none")
)

csd(
  x,
  window = nextpow2(sqrt(NROW(x))),
  overlap = 0.5,
  nfft = ifelse(isScalar(window), window, length(window)),
  fs = 1,
  detrend = c("long-mean", "short-mean", "long-linear", "short-linear", "none")
)

Arguments

Details

cpsd estimates the cross power spectral density function using Welch’s overlapped averaged periodogram method [1].

Value

A list containing the following elements:

freq

vector of frequencies at which the spectral variables are estimated. If x is numeric, power from negative frequencies is added to the positive side of the spectrum, but not at zero or Nyquist (fs/2) frequencies. This keeps power equal in time and spectral domains. If x is complex, then the whole frequency range is returned.

cross

NULL for univariate series. For multivariate series, a matrix containing the squared coherence between different series. Column $i + (j - 1) * (j - 2)/2$ of coh contains the cross-spectral estimates between columns $i$ and $j$ of $x$, where $i < j$.

Note

The function cpsd (and its deprecated alias csd) is a wrapper for the function pwelch, which is more complete and more flexible.

Author(s)

Peter V. Lanspeary, [email protected].
Conversion to R by Geert van Boxtel, [email protected].

References

[1] Welch, P.D. (1967). The use of Fast Fourier Transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE Transactions on Audio and Electroacoustics, AU-15 (2): 70–73.

Examples

fs <- 1000
f <- 250
t <- seq(0, 1 - 1/fs, 1/fs)
s1 <- sin(2 * pi * f * t) + runif(length(t))
s2 <- sin(2 * pi * f * t - pi / 3) + runif(length(t))
rv <- cpsd(cbind(s1, s2), fs = fs)
plot(rv$freq, 10 * log10(rv$cross), type="l", xlab = "Frequency",
     ylab = "Cross Spectral Density (dB)")

Description

Compute the Chirp Z-transform along a spiral contour on the z-plane.

Usage

czt(x, m = NROW(x), w = exp(complex(real = 0, imaginary = -2 * pi/m)), a = 1)

Arguments

Details

The chirp Z-transform (CZT) is a generalization of the discrete Fourier transform (DFT). While the DFT samples the Z plane at uniformly-spaced points along the unit circle, the chirp Z-transform samples along spiral arcs in the Z-plane, corresponding to straight lines in the S plane. The DFT, real DFT, and zoom DFT can be calculated as special cases of the CZT[1]. For the specific case of the DFT, a = 0, m = NCOL(x), and w = 2 * pi / m[2, p. 656].

Value

Chirp Z-transform, returned as a vector or matrix.

Author(s)

Daniel Gunyan.
Conversion to R by Geert van Boxtel, [email protected].

References

[1] https://en.wikipedia.org/wiki/Chirp_Z-transform
[2]Oppenheim, A.V., Schafer, R.W., and Buck, J.R. (1999). Discrete-Time Signal Processing, 2nd edition. Prentice-Hall.

Examples

fs <- 1000                                           # sampling frequency
secs <- 10                                           # number of seconds
t <- seq(0, secs, 1/fs)                              # time series
x <- sin(100 * 2 * pi * t) + runif(length(t))        # 100 Hz signal + noise
m <- 32                                              # n of points desired
f0 <- 75; f1 <- 175;                                 # desired freq range
w <- exp(-1i * 2 * pi * (f1 - f0) / ((m - 1) * fs))  # freq step of f1-f0/m
a <- exp(1i * 2 * pi * f0 / fs);                     # starting at freq f0
y <- czt(x, m, w, a)

# compare DFT and FFT
fs <- 1000
h <- as.numeric(fir1(100, 125/(fs / 2), type = "low"))
m <- 1024
y <- stats::fft(postpad(h, m))

f1 <- 75; f2 <- 175;
w <- exp(-1i * 2 * pi * (f2 - f1) / (m * fs))
a <- exp(1i * 2 * pi * f1 / fs)
z <- czt(h, m, w, a)

fn <- seq(0, m - 1, 1) / m
fy <- fs * fn
fz = (f2 - f1) * fn + f1
plot(fy, 10 * log10(abs(y)), type = "l", xlim = c(50, 200),
  xlab = "Frequency", ylab = "Magnitude (dB")
lines(fz, 10 * log10(abs(z)), col = "red")
legend("topright", legend = c("FFT", "CZT"), col=1:2, lty = 1)

Description

Compute the unitary discrete cosine transform of a signal.

Usage

dct(x, n = NROW(x))

Arguments

Details

The discrete cosine transform (DCT) is closely related to the discrete Fourier transform. You can often reconstruct a sequence very accurately from only a few DCT coefficients. This property is useful for applications requiring data reduction.

The DCT has four standard variants. This function implements the DCT-II according to the definition in [1], which is the most common variant, and the original variant first proposed for image processing.

Value

Discrete cosine transform, returned as a vector or matrix.

Note

The transform is faster if x is real-valued and has even length.

Author(s)

Paul Kienzle, [email protected].
Conversion to R by Geert van Boxtel, [email protected].

References

[1] https://en.wikipedia.org/wiki/Discrete_cosine_transform

See Also

idct

Examples

x <- matrix(seq_len(100) + 50 * cos(seq_len(100) * 2 * pi / 40))
X <- dct(x)

# Find which cosine coefficients are significant (approx.)
# zero the rest
nsig <- which(abs(X) < 1)
N <- length(X) - length(nsig) + 1
X[nsig] <- 0

# Reconstruct the signal and compare it to the original signal.
xx <- idct(X)
plot(x, type = "l")
lines(xx, col = "red")
legend("bottomright", legend = c("Original", paste("Reconstructed, N =", N)),
       lty = 1, col = 1:2)

Description

Compute the two-dimensional discrete cosine transform of a matrix.

Usage

dct2(x, m = NROW(x), n = NCOL(x))

Arguments

Details

The discrete cosine transform (DCT) is closely related to the discrete Fourier transform. It is a separable linear transformation; that is, the two-dimensional transform is equivalent to a one-dimensional DCT performed along a single dimension followed by a one-dimensional DCT in the other dimension.

Value

m-by-n numeric discrete cosine transformed matrix.

Author(s)

Paul Kienzle, [email protected].
Conversion to R by Geert van Boxtel, [email protected].

See Also

idct2

Examples

A <- matrix(runif(100), 10, 10)
B <- dct2(A)

Description

Compute the discrete cosine transform matrix.

Usage

dctmtx(n)

Arguments

Details

A DCT transformation matrix is useful for doing things like JPEG image compression, in which an 8x8 DCT matrix is applied to non-overlapping blocks throughout an image and only a sub-block on the top left of each block is kept. During restoration, the remainder of the block is filled with zeros and the inverse transform is applied to the block.

The two-dimensional DCT of A can be computed as D %*% A %*% t(D). This computation is sometimes faster than using dct2, especially if you are computing a large number of small DCTs, because D needs to be determined only once. For example, in JPEG compression, the DCT of each 8-by-8 block is computed. To perform this computation, use dctmtx to determine D of input image A, and then calculate each DCT using D %*% A %*% t(D) (where A is each 8-by-8 block). This is faster than calling dct2 for each individual block.

Value

Discrete cosine transform, returned as a vector or matrix.

Author(s)

Paul Kienzle, [email protected].
Conversion to R by Geert van Boxtel, [email protected].

See Also

dct, dct2, idct, idct2

Examples

D <- dctmtx(8)

Description

Downsample a signal by an integer factor.

Usage

decimate(x, q, ftype = c("iir", "fir"), n = ifelse(ftype == "iir", 8, 30))

Arguments

Value

downsampled signal, returned as a vector or matrix.

Author(s)

Paul Kienzle, [email protected].
Conversion to R by Geert van Boxtel, [email protected].

See Also

cheby1, fir1

Examples

t <- seq(0, 2, 0.01)
x <- chirp(t, 2, .5, 10, 'quadratic') + sin(2 * pi * t * 0.4)
w <- 1:121
plot(t[w] * 1000, x[w], type = "h", col = "green")
points(t[w] * 1000, x[w], col = "green")
y = decimate(x, 4)
lines(t[seq(1, 121, 4)] * 1000, y[1:31], type = "h", col = "red")
points(t[seq(1, 121, 4)] * 1000, y[1:31], col = "red")

Description

detrend removes the polynomial trend of order p from the data x.

Usage

detrend(x, p = 1)

Arguments

Value

The detrended data, of same type and dimensions as x

Author(s)

Kurt Hornik, [email protected].
Conversion to R by Geert van Boxtel, [email protected].

Examples

t <- 0:20
x <- 3 * sin(t) + t
y <- detrend(x)
plot(t, x, type = "l", ylim = c(-5, 25), xlab = "", ylab = "")
lines(t, y, col = "red")
lines(t, x - y, lty = 2)
legend('topleft', legend = c('Input Data', 'Detrended Data', 'Trend'),
 col = c(1, 2 ,1), lty = c(1, 1, 2))

Description

Compute the discrete Fourier transform matrix

Usage

dftmtx(n)

Arguments

Details

A discrete Fourier transform matrix is a complex matrix whose matrix product with a vector computes the discrete Fourier transform of the vector. dftmtx takes the FFT of the identity matrix to generate the transform matrix. For a column vector x, y <- dftmtx(n) * x is the same as y <- fft(x, postpad(x, n). The inverse discrete Fourier transform matrix is inv <- Conj(dftmtx(n)) / n.

In general this is less efficient than calling the fft and ifft functions directly.

Value

Fourier transform matrix.

Author(s)

David Bateman, [email protected].
Conversion to R by Geert van Boxtel, [email protected].

See Also

fft, ifft

Examples

x <- seq_len(256)
y1 <- stats::fft(x)
n <- length(x)
y2 <- drop(x %*% dftmtx(n))
mx <- max(abs(y1 - y2))

Description

Reorder the elements of the input vector in digit-reversed order.

Usage

digitrevorder(x, r, index.return = FALSE)

Arguments

Details

This function is useful for pre-ordering a vector of filter coefficients for use in frequency-domain filtering algorithms, in which the fft and ifft transforms are computed without digit-reversed ordering for improved run-time efficiency.

Value

The digit-reversed input vector. If index.return = TRUE, then a list containing the digit-reversed input vector (y, and the digit-reversed indices (i).

Author(s)

Mike Miller.
Conversion to R by Geert van Boxtel, [email protected].

See Also

bitrevorder, fft, ifft

Examples

res <- digitrevorder(0:8, 3)

Description

Compute the Dirichlet or periodic sinc function.

Usage

diric(x, n)

Arguments

Details

y <- diric(x, n) returns the Dirichlet Function of degree n evaluated at the elements of the input array x.

The Dirichlet function, or periodic sinc function, has period $2 \pi$ for odd $N$ and period $4 \pi$ for even $N$. Its maximum value is 1 for all N, and its minimum value is -1 for even N. The magnitude of the function is 1 / N times the magnitude of the discrete-time Fourier transform of the N-point rectangular window.

Value

Output array, returned as a real-valued scalar, vector, matrix, or multidimensional array of the same size as x.

Author(s)

Sylvain Pelissier, [email protected].
Conversion to R by Geert van Boxtel [email protected].

Examples

## Compute and plot the Dirichlet function between -2pi and 2pi for N = 7
## and N = 8. The function has a period of 2pi for odd N and 4pi for even N.
x <- seq(-2*pi, 2*pi, len = 301)
d7 <- diric(x, 7)
d8 <- diric(x, 8)
op <- par(mfrow = c(2,1))
plot(x/pi, d7, type="l", main = "Dirichlet function",
     xlab = "", ylab = "N = 7")
plot(x/pi, d8, type="l", ylab = "N = 8", xlab = expression(x / pi))
par(op)

Description

Downsample a signal by an integer factor.

Usage

downsample(x, n, phase = 0)

Arguments

Details

For most signals you will want to use decimate instead since it prefilters the high frequency components of the signal and avoids aliasing effects.

Value

Downsampled signal, returned as a vector or matrix.

Author(s)

Paul Kienzle, [email protected].
Conversion to R by Geert van Boxtel, [email protected].

See Also

decimate, resample

Examples

x <- seq_len(10)
xd <- downsample(x, 3)     # returns 1  4  7 10
xd <- downsample(x, 3, 2)  # returns 3 6 9

x <- matrix(seq_len(12), 4, 3, byrow = TRUE)
xd <- downsample(x, 3)

Description

Compute the discrete sine transform of a signal.

Usage

dst(x, n = NROW(x))

Arguments

Details

The discrete sine transform (DST) is closely related to the discrete Fourier transform. but using a purely real matrix. It is equivalent to the imaginary parts of a DFT of roughly twice the length.

The DST has four standard variants. This function implements the DCT-I according to the definition in [1], which is the most common variant, and the original variant first proposed for image processing.

The 'Matlab' documentation for the DST warns that the use of the function is not recommended. They do not state the reason why, but it is likely that use of the discrete cosine transform (DCT)is preferred for image processing. Because cos(0) is 1, the first coefficient of the DCT (II) is the mean of the values being transformed. This makes the first coefficient of each 8x8 block represent the average tone of its constituent pixels, which is obviously a good start. Subsequent coefficients add increasing levels of detail, starting with sweeping gradients and continuing into increasingly fiddly patterns, and it just so happens that the first few coefficients capture most of the signal in photographic images. Sin(0) is 0, so the DSTs start with an offset of 0.5 or 1, and the first coefficient is a gentle mound rather than a flat plain. That is unlikely to suit ordinary images, and the result is that DSTs require more coefficients than DCTs to encode most blocks. This explanation was provided by Douglas Bagnall on Stackoverflow.

Value

Discrete sine transform, returned as a vector or matrix.

Author(s)

Paul Kienzle, [email protected].
Conversion to R by Geert van Boxtel, [email protected].

References

[1] https://en.wikipedia.org/wiki/Discrete_sine_transform

See Also

idst

Examples

x <- matrix(seq_len(100) + 50 * cos(seq_len(100) * 2 * pi / 40))
ct <- dct(x)
st <- dst(x)

Description

Compute the single-level discrete wavelet transform of a signal

Usage

dwt(x, wname = "d8", lo = NULL, hi = NULL)

wfilters(wname)

Arguments

Details

This function is only included because of compatibility with the 'Octave' 'signal' package. Specialized packages exist in R to perform the discrete wavelet transform, e.g., the wavelets package [1]. this function recognizes only a few wavelet names, namely those for which scale coefficients are available (Daubechies [2] and Coiflet [3]).

The wavelet and scaling coefficients are returned by the function wfilters, which returns the coefficients for reconstruction filters associated with the wavelet wname. Decomposition filters are the time reverse of the reconstruction filters (see examples).

Value

A list containing two numeric vectors:

a

approximation (average) coefficients, obtained from convolving x with the scaling (low-pass) filter lo, and then downsampled (keep the even-indexed elements).

d

detail (difference) coefficients, obtained from convolving x with the wavelet (high-pass) filter hi, and then downsampled (keep the even-indexed elements).

Note

The notations g and h are often used to denote low-pass (scaling) and high-pass (wavelet) coefficients, respectively, but inconsistently. Ref [4] uses it, as does the R wavelets package. 'Octave' uses the reverse notation. To avoid confusion, more neutral terms are used here.

There are two naming schemes for wavelet names in use. For instance for Daubechies wavelets (d), dN using the length or number of taps, and dbA referring to the number of vanishing moments. So d4 and db2 are the same wavelet transform. This function uses the formed (dN) notation; 'Matlab' uses the latter (dbA).

Author(s)

Lukas F. Reichlin.
Conversion to R by Geert van Boxtel, [email protected].

References

[1] https://CRAN.R-project.org/package=wavelets

[2] https://en.wikipedia.org/wiki/Daubechies_wavelet

[3] https://en.wikipedia.org/wiki/Coiflet

[4] https://en.wikipedia.org/wiki/Discrete_wavelet_transform

Examples

# get Coiflet 30 coefficients
wv <- wfilters('c30')
lo <- rev(wv$lo)
hi <- rev(wv$hi)

# general time-varying signal
time <- 1
fs <- 1000
x <- seq(0,time, length.out=time*fs)
y <- c(cos(2*pi*100*x)[1:300], cos(2*pi*50*x)[1:300],
       cos(2*pi*25*x)[1:200], cos(2*pi*10*x)[1:200])
op <- par(mfrow = c(3,1))
plot(x, y, type = "l", xlab = "Time", ylab = "Amplitude",
     main = "Original signal")
wt <- dwt(y, wname = NULL, lo, hi)

x2 <- seq(1, length(x) - length(hi) + 1, 2)
plot(x2, wt$a, type = "h", xlab = "Time", ylab = "",
    main = "Approximation coefficients")
plot(x2, wt$d, type = "h", xlab = "Time", ylab = "",
     main = "Detail coefficients")
par (op)

Description

Compute the transfer function coefficients of an elliptic filter.

Usage

ellip(n, ...)

## S3 method for class 'FilterSpecs'
ellip(n, Rp = n$Rp, Rs = n$Rs, w = n$Wc, type = n$type, plane = n$plane, ...)

## Default S3 method:
ellip(
  n,
  Rp,
  Rs,
  w,
  type = c("low", "high", "stop", "pass"),
  plane = c("z", "s"),
  output = c("Arma", "Zpg", "Sos"),
  ...
)

Arguments

Details

An elliptic filter is a filter with equalized ripple (equiripple) behavior in both the passband and the stopband. The amount of ripple in each band is independently adjustable, and no other filter of equal order can have a faster transition in gain between the passband and the stopband, for the given values of ripple.