gsignal: Signal Processing

Description

R implementation of the 'Octave' package 'signal', containing a variety of signal processing tools, such as signal generation and measurement, correlation and convolution, filtering, filter design, filter analysis and conversion, power spectrum analysis, system identification, decimation and sample rate change, and windowing.

Details

The package 'gsignal' is an implementation in R of the 'Octave' package 'signal'. It provides signal processing algorithms for use with R, include the creation of waveforms, FIR and IIR filter design, spectral analysis, Fourier and other transforms, window functions, and resampling and rate changing.

The 'Octave' package 'signal' release 1.4.1 (2019-02-08) was used to port the functions to R. Note that compatibility of function parameters and return values was not explicitly pursued.

Author(s)

Maintainer: Geert van Boxtel G.J.M.vanBoxtel@gmail.com (Maintainer)

Authors:

• Tom Short tshort@eprisolutions.com (Author of 'signal' package)

• Paul Kienzle (Majority of the original sources)

Other contributors:

• Ben Abbott [contributor]

• Juan Aguado [contributor]

• Muthiah Annamalai [contributor]

• Leonardo Araujo [contributor]

• William Asquith [contributor]

• David Bateman [contributor]

• David Billinghurst [contributor]

• Juan Pablo Carbajal [contributor]

• André Carezia [contributor]

• Vincent Cautaerts [contributor]

• Eric Chassande-Mottin [contributor]

• Luca Citi [contributor]

• Dave Cogdell [contributor]

• Carlo de Falco [contributor]

• Carne Draug [contributor]

• Pascal Dupuis [contributor]

• John W. Eaton [contributor]

• R.G.H Eschauzier [contributor]

• Andrew Fitting [contributor]

• Alan J. Greenberger [contributor]

• Mike Gross [contributor]

• Daniel Gunyan [contributor]

• Kai Habel [contributor]

• Kurt Hornik [contributor]

• Jake Janovetz [contributor]

• Alexander Klein [contributor]

• Peter V. Lanspeary [contributor]

• Bill Lash [contributor]

• Friedrich Leissh [contributor]

• Laurent S. Mazet [contributor]

• Mike Miller [contributor]

• Petr Mikulik [contributor]

• Paolo Neis [contributor]

• Georgios Ouzounis [contributor]

• Sylvain Pelissier [contributor]

• Francesco Potortì [contributor]

• Charles Praplan [contributor]

• Lukas F. Reichlin [contributor]

• Tony Richardson [contributor]

• Asbjorn Sabo [contributor]

• Thomas Sailer [contributor]

• Rolf Schirmacher [contributor]

• Rolf Schirmacher [contributor]

• Ivan Selesnick [contributor]

• Julius O. Smith III [contributor]

• Peter L. Soendergaard [contributor]

• Quentin Spencer [contributor]

• Doug Stewart [contributor]

• P. Sudeepam [contributor]

• Stefan van der Walt [contributor]

• Andreas Weber [contributor]

• P. Sudeepam [contributor]

• Andreas Weingessel [contributor]

References

https://wiki.octave.org/Signal_package,
https://octave.sourceforge.io/signal/

See Also

Useful links:

• https://github.com/gjmvanboxtel/gsignal

• Report bugs at https://github.com/gjmvanboxtel/gsignal/issues

Autoregressive moving average (ARMA) model

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, tshort@eprisolutions.com,
adapted by Geert van Boxtel, gjmvanboxtel@gmail.com.

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)

Filter specifications

Description

Filter specifications, including order, frequency cutoff, type, and possibly others.

Usage

FilterSpecs(n, Wc, type, ...)

Arguments

Value

A list of class 'FilterSpecs' with the following list elements (repeats of the input arguments):

n

filter order

Wc

cutoff frequency

type

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

...

other filter description characteristics, possibly including Rp for dB of pass band ripple or Rs for dB of stop band ripple, depending on filter type (Butterworth, Chebyshev, etc.).

Author(s)

Tom Short, tshort@eprisolutions.com,
renamed and adapted by Geert van Boxtel, gjmvanboxtel@gmail.com

See Also

filter, butter and buttord, cheby1 and cheb1ord, ellip and ellipord.

Examples

filt <- FilterSpecs(3, 0.1, "low")

Moving average (MA) model

Description

Create an MA model representing a filter or system model

Usage

Ma(b)

Arguments

Value

A list of class Ma with the polynomial coefficients

Author(s)

Tom Short, tshort@eprisolutions.com

See Also

See also Arma

Examples

f <- Ma(b = c(1, 2, 1) / 3)
freqz(f)
zplane(f)

Second-order sections

Description

Create or convert filter models to second-order sections form.

Usage

Sos(sos, g = 1)

as.Sos(x, ...)

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

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

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

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

Arguments

Details

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

Value

A list of class Sos with the following list elements:

sos

second-order section representation of the model, returned as an L x 6 matrix, one row for each section 1:L. Each row consists of an [B, A], pair, where B = c(b0, b1, b2), and A = c(1, a1, a2), the filer coefficients for each section. Each b0 entry must be nonzero for each section.

g

overall gain factor that scales any one of the B_i vectors. Default: 1

Author(s)

Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

See Also

Arma, Ma, Zpg

Examples

ba <- butter(3, 0.2)
sos <- as.Sos(ba)

Zero pole gain model

Description

Create an zero pole gain model of an ARMA filter, or convert other forms to a Zpg model.

Usage

Zpg(z, p, g)

as.Zpg(x, ...)

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

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

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

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

Arguments

Details

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

Value

A list of class Zpg with the following list elements:

z

complex vector of the zeros of the model

p

complex vector of the poles of the model

g

gain of the model

Author(s)

Tom Short, tshort@eprisolutions.com,
adapted by Geert van Boxtel, gjmvanboxtel@gmail.com.

See Also

See also Arma

Examples

## design notch filter at pi/4 radians = 0.5/4 = 0.125 * fs
w = pi/4
# 2 poles, 2 zeros
# zeroes at r = 1
r <- 1
z1 <- r * exp(1i * w)
z2 <- r * exp(1i * -w)
# poles at r = 0.9
r = 0.9
p1 <- r * exp(1i * w)
p2 <- r * exp(1i * -w)

zpg <- Zpg(c(z1, z2), c(p1, p2), 1)
zplane(zpg)
freqz(zpg)

## Sharper edges: increase distance between zeros and poles
r = 0.8
p1 <- r * exp(1i * w)
p2 <- r * exp(1i * -w)
zpg <- Zpg(c(z1, z2), c(p1, p2), 1)
zplane(zpg)
freqz(zpg)

Power spectrum of AR model

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, pvl@mecheng.adelaide.edu.au.
Conversion to R by Geert van Boxtel, gjmvanboxtel@gmail.com

Examples

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

Autoregressive model coefficients - Burg's method

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, pvl@mecheng.adelaide.edu.au.
Conversion to R by Geert van Boxtel, gjmvanboxtel@gmail.com.

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)

Autoregressive model coefficients - Yule-Walker method

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, pkienzle@users.sf.net,
Peter V. Lanspeary, pvl@mecheng.adelaide.edu.au.
Conversion to R by Geert van Boxtel, gjmvanboxtel@gmail.com.

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)

Modified Bartlett-Hann window

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, Andreas.Weingessel@ci.tuwien.ac.at. Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

See Also

bartlett, hann, hamming

Examples

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

Bartlett window

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, Andreas.Weingessel@ci.tuwien.ac.at. Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

See Also

triang

Examples

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

Bessel analog low-pass filter prototype

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, t.sailer@alumni.ethz.ch.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

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)

Bessel analog filter design

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, pkienzle@users.sf.net,
Doug Stewart, dastew@sympatico.ca,
Thomas Sailer, t.sailer@alumni.ethz.ch.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

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)

Bilinear transformation

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 pkienzle@users.sf.net. Conversion to R by Tom Short, adapted by Geert van Boxtel G.J.M.vanBoxtel@gmail.com.

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)

Permute input to bit-reversed order

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, G.J.M.vanBoxtel@gmail.com.

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)

Blackman window

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, Andreas.Weingessel@ci.tuwien.ac.at.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

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")

Blackman-Harris window

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, sylvain.pelissier@gmail.com.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

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")

Blackman-Nuttall window

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, muthiah.annamalai@uta.edu.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

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")

Bohman window

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, sylvain.pelissier@gmail.com.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

See Also

triang

Examples

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

Rectangular window

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, pkienzle@users.sf.net.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

See Also

triang

Examples

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

Buffer signal vector into matrix of data segments

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, adb014@gmail.com.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com

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)
}

Butterworth filter prototype

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, carandraug+dev@gmail.com.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

Examples

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

Butterworth filter design

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, pkienzle@users.sf.net,
Doug Stewart, dastew@sympatico.ca,
Alexander Klein, alexander.klein@math.uni-giessen.de,
John W. Eaton.
Conversion to R by Tom Short,
adapted by Geert van Boxtel G.J.M.vanBoxtel@gmail.com.

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)

Butterworth filter order and cutoff frequency

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, G.J.M.vanBoxtel@gmail.com.

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)

Complex cepstral analysis

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, G.J.M.vanBoxtel@gmail.com.

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")

Circular convolution

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, G.J.M.vanBoxtel@gmail.com.

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)

Chebyshev polynomials

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, acarezia@uol.com.br.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

Examples

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

Chebyshev Type I filter prototype

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, carandraug+dev@gmail.com.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

Examples

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

Chebyshev Type I filter order

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, G.J.M.vanBoxtel@gmail.com.

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)

Chebyshev Type II filter prototype

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, carandraug+dev@gmail.com.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

Examples

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

Chebyshev Type II filter order

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, G.J.M.vanBoxtel@gmail.com.

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)

Chebyshev window

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, acarezia@uol.com.br.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

Examples

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

Chebyshev Type I filter design

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, pkienzle@users.sf.net,
Doug Stewart, dastew@sympatico.ca.
Conversion to R Tom Short,
adapted by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

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")

Chebyshev Type II filter design

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, pkienzle@users.sf.net,
Doug Stewart, dastew@sympatico.ca.
Conversion to R Tom Short,
adapted by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

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")

Chirp signal

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, pkienzle@users.sf.net,
Mike Miller.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

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)

Constrained L2 bandpass FIR filter design

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, G.J.M.vanBoxtel@gmail.com.

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)

Cluster Segments

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, carbajal@ifi.uzh.ch.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

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)

Complex Morlet Wavelet

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, sylvain.pelissier@gmail.com.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

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))

Convolution and polynomial multiplication

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, arichard@stark.cc.oh.us, adapted by John W. Eaton.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

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')

2-D convolution

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, G.J.M.vanBoxtel@gmail.com.

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")

Convolution matrix

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 adb014@gmail.com.
Conversion to R by Geert van Boxtel G.J.M.vanBoxtel@gmail.com.

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)

Complex conjugate pairs

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, G.J.M.vanBoxtel@gmail.com.

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)))

Sort complex conjugate pairs and real

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, G.J.M.vanBoxtel@gmail.com.

See Also

cplxpair

Examples

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

Cross power spectral density

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, pvl@mecheng.adelaide.edu.au.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

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)")

Chirp Z-transform

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, G.J.M.vanBoxtel@gmail.com.

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)

Discrete Cosine Transform

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, pkienzle@users.sf.net.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

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)

2-D Discrete Cosine Transform

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, pkienzle@users.sf.net.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

See Also

idct2

Examples

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

Discrete Cosine Transform Matrix

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, pkienzle@users.sf.net.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

See Also

dct, dct2, idct, idct2

Examples

D <- dctmtx(8)

Decrease sample rate

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, pkienzle@users.sf.net.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

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")

Remove Polynomial Trend

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, Kurt.Hornik@wu-wien.ac.at.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

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))

Discrete Fourier Transform Matrix

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, adb014@gmail.com.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

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))

Permute input to digit-reversed order

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, G.J.M.vanBoxtel@gmail.com.

See Also

bitrevorder, fft, ifft

Examples

res <- digitrevorder(0:8, 3)

Dirichlet function

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, sylvain.pelissier@gmail.com.
Conversion to R by Geert van Boxtel G.J.M.vanBoxtel@gmail.com.

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)

Decrease sample rate

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, pkienzle@users.sf.net.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

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)

Discrete Sine Transform

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, pkienzle@users.sf.net.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

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)

1-D Discrete Wavelet Transform

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, G.J.M.vanBoxtel@gmail.com.

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)

Elliptic filter design

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.

As the ripple in the stopband approaches zero, the filter becomes a type I Chebyshev filter. As the ripple in the passband approaches zero, the filter becomes a type II Chebyshev filter and finally, as both ripple values approach zero, the filter becomes a Butterworth filter.

Because ellip is generic, it can be extended to accept other inputs, using ellipord 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)

Paulo Neis, p_neis@yahoo.com.br,
adapted by Doug Stewart, dastew@sympatico.ca.
Conversion to R Tom Short,
adapted by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

References

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

See Also

Arma, filter, butter, cheby1, ellipord

Examples

## compare the frequency responses of 5th-order Butterworth
## and elliptic filters.
bf <- butter(5, 0.1)
ef <- ellip(5, 3, 40, 0.1)
bfr <- freqz(bf)
efr <- freqz(ef)
plot(bfr$w, 20 * log10(abs(bfr$h)), type = "l", ylim = c(-80, 0),
     xlab = "Frequency (Rad)", ylab = c("dB"), lwd = 2,
     main = paste("Elliptic versus Butterworth filter",
     "low-pass -3 dB cutoff at 0.1 rad", sep = "\n"))
lines(efr$w, 20 * log10(abs(efr$h)), col = "red", lwd = 2)
legend ("topright", legend = c("Butterworh", "Elliptic"),
        lty = 1, lwd = 2, col = 1:2)

Low-pass analog elliptic filter

Description

Return the zeros, poles and gain of an analog elliptic low-pass filter prototype.

Usage

ellipap(n, Rp, Rs)

Arguments

Details

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

Value

list of class Zpg containing zeros, poles and gain of the filter.

Author(s)

Carne Draug, carandraug+dev@gmail.com. Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

Examples

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

Elliptic Filter Order

Description

Compute elliptic filter order and cutoff for the desired response characteristics.

Usage

ellipord(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, one of "low", "high", "stop", or "pass".

Rp

dB of passband ripple.

Rs

dB of stopband ripple.

Author(s)

Paulo Neis, p_neis@yahoo.com.br,
adapted by Charles Praplan.
Conversion to R by Tom Short,
adapted by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

See Also

buttord, cheb1ord, cheb2ord, ellip

Examples

fs <- 10000
spec <- ellipord(1000/(fs/2), 1200/(fs/2), 0.5, 29)
ef <- ellip(spec)
hf <- freqz(ef, fs = fs)
plot(c(0, 1000, 1000, 0, 0), c(0, 0, -0.5, -0.5, 0),
     type = "l", xlab = "Frequency (Hz)", ylab = "Attenuation (dB)",
     col = "red", ylim = c(-35,0), xlim = c(0,2000))
lines(c(5000, 1200, 1200, 5000, 5000), c(-1000, -1000, -29, -29, -1000),
      col = "red")
lines(hf$w, 20*log10(abs(hf$h)))

FFT-based convolution

Description

Convolve two vectors using the FFT for computation.

Usage

fftconv(x, y, n = NULL)

Arguments

Details

The computation uses the FFT by calling the function fftfilt. If the optional argument n is specified, an n-point overlap-add FFT is used.

Value

Convoluted signal, specified as a a vector of length equal to length (x) + length (y) - 1. If x and y are the coefficient vectors of two polynomials, the returned value is the coefficient vector of the product polynomial.

Author(s)

Kurt Hornik, Kurt.Hornik@wu-wien.ac.at,
adapted by John W. Eaton.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

See Also

conv, conv2

Examples

u <- rep(1L, 3)
v <- c(1, 1, 0, 0, 0, 1, 1)
w1 <- conv(u, v)              # time-domain convolution
w2 <- fftconv(u, v)           # frequency domain convolution
all.equal(w1, w2)             # same results

FFT-based FIR filtering

Description

FFT-based FIR filtering using the overlap-add method.

Usage

fftfilt(b, x, n = NULL)

## Default S3 method:
fftfilt(b, x, n = NULL)

## S3 method for class 'Ma'
fftfilt(b, x, n = NULL)

Arguments

Details

This function combines two important techniques to speed up filtering of long signals, the overlap-add method, and FFT convolution. The overlap-add method is used to break long signals into smaller segments for easier processing or preventing memory problems. FFT convolution uses the overlap-add method together with the Fast Fourier Transform, allowing signals to be convolved by multiplying their frequency spectra. For filter kernels longer than about 64 points, FFT convolution is faster than standard convolution, while producing exactly the same result.

The overlap-add technique works as follows. When an N length signal is convolved with a filter kernel of length M, the output signal is N + M - 1 samples long, i.e., the signal is expanded 'to the right'. The signal is then broken into k smaller segments, and the convolution of each segment with the f kernel will have a result of length N / k + M -1. The individual segments are then added together. The rightmost M - 1 samples overlap with the leftmost M - 1 samples of the next segment. The overlap-add method produces exactly the same output signal as direct convolution.

FFT convolution uses the principle that multiplication in the frequency domain corresponds to convolution in the time domain. The input signal is transformed into the frequency domain using the FFT, multiplied by the frequency response of the filter, and then transformed back into the time domain using the inverse FFT. With FFT convolution, the filter kernel can be made very long, with very little penalty in execution time.

Value

The filtered signal, returned as a vector or matrix with the same dimensions as x.

Author(s)

Kurt Hornik, Kurt.Hornik@wu-wien.ac.at,
adapted by John W. Eaton.
Conversion to R by Tom Short,
adapted by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

References

https://en.wikipedia.org/wiki/Overlap-add_method.

See Also

filter

Examples

t <- seq(0, 1, len = 10000)                          # 1 second sample
x <- sin(2* pi * t * 2.3) + 0.25 * rnorm(length(t))  # 2.3 Hz sinusoid+noise
filt <- rep(0.1, 10)                                 # filter kernel
y1 <- filter(filt, 1, x)                             # use normal convolution
y2 <- fftfilt(filt, x)                               # FFT convolution
plot(t, x, type = "l")
lines(t, y1, col = "red")
lines(t, y2, col = "blue")

## use 'filter' with different classes
t <- seq(0, 1, len = 10000)                          # 1 second sample
x <- sin(2* pi * t * 2.3) + 0.25 * rnorm(length(t))  # 2.3 Hz sinusoid+noise
ma <- Ma(rep(0.1, 10))                               # filter kernel
y1 <- filter(ma, x)                                  # convulution filter
y2 <- fftfilt(ma, x)                                 # FFT filter
all.equal(y1, y2)                                    # same result

Zero-frequency shift

Description

Perform a shift in order to move the frequency 0 to the center of the input.

Usage

fftshift(x, MARGIN = 2)

Arguments

Details

If x is a vector of N elements corresponding to N time samples spaced by dt, then fftshift(x) corresponds to frequencies f = c(-seq(ceiling((N-1)/2), 1, -1), 0, (1:floor((N-1)/2))) * df, where df = 1 / (N * dt). In other words, the left and right halves of x are swapped.

If x is a matrix, then fftshift operates on the rows or columns of x, according to the MARGIN argument, i.e. it swaps the the upper and lower halves of the matrix (MARGIN = 1), or the left and right halves of the matrix (MARGIN = 2). Specifying MARGIN = c(1, 2) swaps along both dimensions, i.e., swaps the first quadrant with the fourth, and the second with the third.

Value

vector or matrix with centered frequency.

Author(s)

Vincent Cautaerts, vincent@comf5.comm.eng.osaka-u.ac.jp,
adapted by John W. Eaton.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

See Also

ifftshift

Examples

Xeven <- 1:6
ev <- fftshift(Xeven)   # returns 4 5 6 1 2 3

Xodd <- 1:7
odd <- fftshift(Xodd)   # returns 5 6 7 1 2 3 4

fs <- 100                      # sampling frequency
t <- seq(0, 10 - 1/fs, 1/fs)   # time vector
S <- cos(2 * pi * 15 * t)
n <- length(S)
X <- fft(S)
f <- (0:(n - 1)) * (fs / n);   # frequency range
power <- abs(X)^2 / n          # power
plot(f, power, type="l")
Y <- fftshift(X)
fsh <- ((-n/2):(n/2-1)) * (fs / n)  # zero-centered frequency range
powersh <- abs(Y)^2 / n             # zero-centered power
plot(fsh, powersh, type = "l")

Fast Hartley Transform

Description

Compute the (inverse) Hartley transform of a signal using FFT

Usage

fht(x, n = NROW(x))

ifht(x, n = NROW(x))

Arguments

Details

The Hartley transform is an integral transform closely related to the Fourier transform, but which transforms real-valued functions to real-valued functions. Compared to the Fourier transform, the Hartley transform has the advantages of transforming real functions to real functions (as opposed to requiring complex numbers) and of being its own inverse [1].

This function implements the Hartley transform by calculating the difference between the real- and imaginary-valued parts of the Fourier-transformed signal [1]. The forward and inverse Hartley transforms are the same (except for a scale factor of 1/N for the inverse Hartley transform), but implemented using different functions.

Value

(inverse) Hartley transform, returned as a vector or matrix.

Author(s)

Muthiah Annamalai, muthiah.annamalai@uta.edu.\ Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

References

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

See Also

fft

Examples

# FHT of a 2.5 Hz signal with offset
fs <- 100
secs <- 10
freq <- 2.5
t <- seq(0, secs - 1 / fs, 1 / fs)
x <- 5 * t + 50 * cos(freq * 2 * pi * t)
X <- fht(x)
op <- par(mfrow = c(2, 1))
plot(t, x, type = "l", xlab = "", ylab = "", main = "Signal")
f <- seq(0, fs - (1 / fs), length.out = length(t))
to <- which(f >= 5)[1]
plot(f[1:to], X[1:to], type = "l", xlab = "", ylab = "",
     main = "Hartley Transform")
par(op)

Filter a signal

Description

Apply a 1-D digital filter compatible with 'Matlab' and 'Octave'.

Usage

filter(filt, ...)

## Default S3 method:
filter(filt, a, x, zi = NULL, ...)

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

## S3 method for class 'Ma'
filter(filt, x, ...)

## S3 method for class 'Sos'
filter(filt, x, ...)

## S3 method for class 'Zpg'
filter(filt, x, ...)

Arguments

Details

The filter is a direct form II transposed implementation of the standard linear time-invariant difference equation:

  N                  M
 SUM a(k+1)y(n-k) + SUM b(k+1)x(n-k) = 0;   1 <= n <= length(x)
 k=0                k=0

where N = length(a) - 1 and M = length(b) - 1.

The initial and final conditions for filter delays can be used to filter data in sections, especially if memory limitations are a consideration. See the examples.

Value

The filtered signal, of the same dimensions as the input signal. In case the zi input argument was specified, a list with two elements is returned containing the variables y, which represents the output signal, and zf, which contains the final state vector or matrix.

Author(s)

Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

See Also

filter_zi, sosfilt (preferred because it avoids numerical problems).

Examples

bf <- butter(3, 0.1)                                 # 10 Hz low-pass filter
t <- seq(0, 1, len = 100)                            # 1 second sample
x <- sin(2* pi * t * 2.3) + 0.25 * rnorm(length(t))  # 2.3 Hz sinusoid+noise
z <- filter(bf, x)                                   # apply filter
plot(t, x, type = "l")
lines(t, z, col = "red")

## specify initial conditions
## from Python scipy.signal.lfilter() documentation
t <- seq(-1, 1, length.out =  201)
x <- (sin(2 * pi * 0.75 * t * (1 - t) + 2.1)
      + 0.1 * sin(2 * pi * 1.25 * t + 1)
      + 0.18 * cos(2 * pi * 3.85 * t))
h <- butter(3, 0.05)
lab <- max(length(h$b), length(h$a)) - 1
zi <- filtic(h$b, h$a, rep(1, lab), rep(1, lab))
z1 <- filter(h, x)
z2 <- filter(h, x, zi * x[1])
plot(t, x, type = "l")
lines(t, z1, col = "red")
lines(t, z2$y, col = "green")
legend("bottomright", legend = c("Original signal",
        "Filtered without initial conditions",
        "Filtered with initial conditions"),
       lty = 1, col = c("black", "red", "green"))

Savitzky-Golay filtering

Description

Filter a signal with a Savitzky-Golay FIR filter.

Usage

## S3 method for class 'sgolayFilter'
filter(filt, x, ...)

sgolayfilt(x, p = 3, n = p + 3 - p%%2, m = 0, ts = 1)

Arguments

Details

Savitzky-Golay smoothing filters are typically used to "smooth out" a noisy signal whose frequency span (without noise) is large. They are also called digital smoothing polynomial filters or least-squares smoothing filters. Savitzky-Golay filters perform better in some applications than standard averaging FIR filters, which tend to filter high-frequency content along with the noise. Savitzky-Golay filters are more effective at preserving high frequency signal components but less successful at rejecting noise.

Savitzky-Golay filters are optimal in the sense that they minimize the least-squares error in fitting a polynomial to frames of noisy data.

Value

The filtered signal, of the same dimensions as the input signal.

Author(s)

Paul Kienzle, pkienzle@users.sf.net.
Conversion to R Tom Short,
adapted by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

See Also

sgolay

Examples

# Compare a 5 sample averager, an order-5 butterworth lowpass
# filter (cutoff 1/3) and sgolayfilt(x, 3, 5), the best cubic
# estimated from 5 points.
bf <- butter(5, 1/3)
x <- c(rep(0, 15), rep(10, 10), rep(0, 15))
sg <- sgolayfilt(x)
plot(sg, type="l", xlab = "", ylab = "")
lines(filtfilt(rep(1, 5) / 5, 1, x), col = "red") # averaging filter
lines(filtfilt(bf, x), col = "blue")              # butterworth
points(x, pch = "x")                              # original data
legend("topleft", c("sgolay (3,5)", "5 sample average", "order 5
Butterworth", "original data"), lty=c(1, 1, 1, NA),
pch = c(NA, NA, NA, "x"), col = c(1, "red", "blue", 1))

2-D digital filter

Description

Apply a 2-D digital filter to the data in x.

Usage

filter2(h, x, shape = c("same", "full", "valid"))

Arguments

Details

The filter2 function filters data by taking the 2-D convolution of the input x and the coefficient matrix h rotated 180 degrees. More specifically, filter2(h, x, shape) is equivalent to conv2(x, rot90(h, 2), shape).

Value

The filtered signal, returned as a matrix

Author(s)

Paul Kienzle. Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

See Also

conv2

Examples

op <- par(mfcol = c(1, 2))
x <- seq(-10, 10, length.out = 30)
y <- x
f <- function(x, y) { r <- sqrt(x^2+y^2); 10 * sin(r)/r }
z <- outer(x, y, f)
z[is.na(z)] <- 1
persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue")
title( main = "Original")

h <- matrix(c(1, -2, 1, -2, 3, -2, 1, -2, 1), 3, 3)
zf <-filter2(h, z, 'same')
persp(x, y, zf, theta = 30, phi = 30, expand = 0.5, col = "lightgreen")
title( main = "Filtered")
par(op)

Filter initial conditions

Description

Construct initial conditions for a filter

Usage

filter_zi(filt, ...)

## Default S3 method:
filter_zi(filt, a, ...)

## S3 method for class 'Arma'
filter_zi(filt, ...)

## S3 method for class 'Ma'
filter_zi(filt, ...)

## S3 method for class 'Sos'
filter_zi(filt, ...)

## S3 method for class 'Zpg'
filter_zi(filt, ...)

Arguments

Details

This function computes an initial state for the filter function that corresponds to the steady state of the step response. In other words, it finds the initial condition for which the response to an input of all ones is a constant. Therefore, the results returned by this function can also be obtained using the function filtic by setting x and y to all 1s (see the examples).

A typical use of this function is to set the initial state so that the output of the filter starts at the same value as the first element of the signal to be filtered.

Value

The initial state for the filter, returned as a vector.

Author(s)

Geert van Boxtel, G.J.M.vanBoxtel@gmail.com, converted to R from Python scipy.signal.lfilter_zi.

References

Gustafsson, F. (1996). Determining the initial states in forward-backward filtering. IEEE Transactions on Signal Processing, 44(4), 988 - 992.

See Also

filtic

Examples

## taken from Python scipy.signal.lfilter_zi documentation

h <- butter(5, 0.25)
zi <- filter_zi(h)
y <- filter(h, rep(1, 10), zi)
## output is all 1, as expected.
y2 <- filter(h, rep(1, 10))
## if the zi argument is not given, the output
## does not return the final conditions

x <- c(0.5, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0)
y <- filter(h, x, zi = zi*x[1])
## Note that the zi argument to filter was computed using
## filter_zi and scaled by x[1]. Then the output y has no
## transient until the input drops from 0.5 to 0.0.

## obtain the same results with filtic
lab <- max(length(h$b), length(h$a)) - 1
ic <- filtic(h, rep(1, lab), rep(1, lab))
all.equal(zi, ic)

Zero-phase digital filtering

Description

Forward and reverse filter the signal.

Usage

filtfilt(filt, ...)

## Default S3 method:
filtfilt(filt, a, x, ...)

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

## S3 method for class 'Ma'
filtfilt(filt, x, ...)

## S3 method for class 'Sos'
filtfilt(filt, x, ...)

## S3 method for class 'Zpg'
filtfilt(filt, x, ...)

Arguments

Details

Forward and reverse filtering the signal corrects for phase distortion introduced by a one-pass filter, though it does square the magnitude response in the process. That’s the theory at least. In practice the phase correction is not perfect, and magnitude response is distorted, particularly in the stop band.

Before filtering the input signal is extended with a reflected part of both ends of the signal. The length of this extension is 3 times the filter order. The Gustafsson [1] method is then used to specify the initial conditions used to further handle the edges of the signal.

Value

The filtered signal, normally of the same length of the input signal x, returned as a vector or matrix.

Author(s)

Paul Kienzle, pkienzle@users.sf.net,
Francesco Potortì, pot@gnu.org,
Luca Citi, lciti@essex.ac.uk.
Conversion to R and adapted by Geert van Boxtel G.J.M.vanBoxtel@gmail.com.

References

[1] Gustafsson, F. (1996). Determining the initial states in forward-backward filtering. IEEE Transactions on Signal Processing, 44(4), 988 - 992.

See Also

filter, filter_zi, Arma, Sos, Zpg

Examples

bf <- butter(3, 0.1)                                 # 10 Hz low-pass filter
t <- seq(0, 1, len = 100)                            # 1 second sample
x <- sin(2* pi * t * 2.3) + 0.25 * rnorm(length(t))  # 2.3 Hz sinusoid+noise
z <- filter(bf, x)                                   # apply filter
plot(t, x, type = "l")
lines(t, z, col = "red")
zz <- filtfilt(bf, x)
lines(t, zz, col="blue")
legend("bottomleft", legend = c("original", "filter", "filtfilt"), lty = 1,
 col = c("black", "red", "blue"))

Filter Initial Conditions

Description

Compute the initial conditions for a filter.

Usage

filtic(filt, ...)

## Default S3 method:
filtic(filt, a, y, x = 0, ...)

## S3 method for class 'Arma'
filtic(filt, y, x = 0, ...)

## S3 method for class 'Ma'
filtic(filt, y, x = 0, ...)

## S3 method for class 'Sos'
filtic(filt, y, x = 0, ...)

## S3 method for class 'Zpg'
filtic(filt, y, x = 0, ...)

Arguments

Details

This function computes the same values that would be obtained from the function filter given past inputs x and outputs y.

The vectors x and y contain the most recent inputs and outputs respectively, with the newest values first:

x = c(x(-1), x(-2), ... x(-nb)); nb = length(b)-1
y = c(y(-1), y(-2), ... y(-na)); na = length(a)-a

If length(x) < nb then it is zero padded. If length(y) < na then it is zero padded.

Value

Initial conditions for filter specified by filt, input vector x, and output vector y, returned as a vector.

Author(s)

David Billinghurst, David.Billinghurst@riotinto.com.
Adapted and converted to R by Geert van Boxtel G.J.M.vanBoxtel@gmail.com.

See Also

filter, sosfilt, filtfilt, filter_zi

Examples

## Simple low pass filter
b <- c(0.25, 0.25)
a <- c(1.0, -0.5)
ic <- filtic(b, a, 1, 1)

## Simple high pass filter
b <- c(0.25, -0.25)
a <- c(1.0, 0.5)
ic <- filtic(b, a, 0, 1)

## Example from Python scipy.signal.lfilter() documentation
t <- seq(-1, 1, length.out =  201)
x <- (sin(2 * pi * 0.75 * t * (1 - t) + 2.1)
      + 0.1 * sin(2 * pi * 1.25 * t + 1)
      + 0.18 * cos(2 * pi * 3.85 * t))
h <- butter(3, 0.05)
l <- max(length(h$b), length(h$a)) - 1
zi <- filtic(h, rep(1, l), rep(1, l))
z <- filter(h, x, zi * x[1])

Find local extrema

Description

Return peak values and their locations of the vector data.

Usage

findpeaks(
  data,
  MinPeakHeight = .Machine$double.eps,
  MinPeakDistance = 1,
  MinPeakWidth = 1,
  MaxPeakWidth = Inf,
  DoubleSided = FALSE
)

Arguments

Details

Peaks of a positive array of data are defined as local maxima. For double-sided data, they are maxima of the positive part and minima of the negative part. data is expected to be a one-dimensional vector.

Value

A list containing the following elements:

pks

The value of data at the peaks.

loc

The index indicating the position of the peaks.

parabol

A list containing the parabola fitted to each returned peak. The list has two fields, x and pp. The field pp contains the coefficients of the 2nd degree polynomial and x the extrema of the interval where it was fitted.

height

The estimated height of the returned peaks (in units of data).

baseline

The height at which the roots of the returned peaks were calculated (in units of data).

roots

The abscissa values (in index units) at which the parabola fitted to each of the returned peaks realizes its width as defined below.

Author(s)

Juan Pablo Carbajal, carbajal@ifi.uzh.ch.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

Examples

### demo 1
t <- 2 * pi * seq(0, 1,length = 1024)
y <- sin(3.14 * t) + 0.5 * cos(6.09 * t) +
     0.1 * sin(10.11 * t + 1 / 6) + 0.1 * sin(15.3 * t + 1 / 3)

data1 <- abs(y) # Positive values
peaks1 <- findpeaks(data1)

data2 <- y # Double-sided
peaks2 <- findpeaks(data2, DoubleSided = TRUE)
peaks3 <- findpeaks (data2, DoubleSided = TRUE, MinPeakHeight = 0.5)

op <- par(mfrow=c(1,2))
plot(t, data1, type="l", xlab="", ylab="")
points(t[peaks1$loc], peaks1$pks, col = "red", pch = 1)
plot(t, data2, type = "l", xlab = "", ylab = "")
points(t[peaks2$loc], peaks2$pks, col = "red", pch = 1)
points(t[peaks3$loc], peaks3$pks, col = "red", pch = 4)
legend ("topleft", "0: >2*sd, x: >0.5", bty = "n",
        text.col = "red")
par (op)
title("Finding the peaks of smooth data is not a big deal")

## demo 2
t <- 2 * pi * seq(0, 1, length = 1024)
y <- sin(3.14 * t) + 0.5 * cos(6.09 * t) + 0.1 *
     sin(10.11 * t + 1 / 6) + 0.1 * sin(15.3 * t + 1 / 3)
data <- abs(y + 0.1*rnorm(length(y),1))   # Positive values + noise
peaks1 <- findpeaks(data, MinPeakHeight=1)
dt <- t[2]-t[1]
peaks2 <- findpeaks(data, MinPeakHeight=1, MinPeakDistance=round(0.5/dt))
op <- par(mfrow=c(1,2))
plot(t, data, type="l", xlab="", ylab="")
points (t[peaks1$loc],peaks1$pks,col="red", pch=1)
plot(t, data, type="l", xlab="", ylab="")
points (t[peaks2$loc],peaks2$pks,col="red", pch=1)
par (op)
title(paste("Noisy data may need tuning of the parameters.\n",
            "In the 2nd example, MinPeakDistance is used\n",
            "as a smoother of the peaks"))

Window-based FIR filter design

Description

FIR filter coefficients for a filter with the given order and frequency cutoff.

Usage

fir1(
  n,
  w,
  type = c("low", "high", "stop", "pass", "DC-0", "DC-1"),
  window = hamming(n + 1),
  scale = TRUE
)

Arguments

Value

The FIR filter coefficients, a vector of length n + 1, of class Ma.

Author(s)

Paul Kienzle, pkienzle@users.sf.net, Conversion to R Tom Short,
adapted by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

References

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

See Also

Ma, filter, fftfilt, fir2

Examples

freqz(fir1(40, 0.3))
freqz(fir1(10, c(0.3, 0.5), "stop"))
freqz(fir1(10, c(0.3, 0.5), "pass"))

Frequency sampling-based FIR filter design

Description

Produce a FIR filter with arbitrary frequency response over frequency bands.

Usage

fir2(n, f, m, grid_n = 512, ramp_n = NULL, window = hamming(n + 1))

Arguments

Details

The function linearly interpolates the desired frequency response onto a dense grid and then uses the inverse Fourier transform and a Hamming window to obtain the filter coefficients.

Value

The FIR filter coefficients, a vector of length n + 1, of class Ma.

Author(s)

Paul Kienzle, pkienzle@users.sf.net.
Conversion to R Tom Short,
adapted by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

See Also

Ma, filter, fftfilt, fir1

Examples

f <- c(0, 0.3, 0.3, 0.6, 0.6, 1)
m <- c(0, 0, 1, 1/2, 0, 0)
fh <- freqz(fir2(100, f, m))
op <- par(mfrow = c(1, 2))
plot(f, m, type = "b", ylab = "magnitude", xlab = "Frequency")
lines(fh$w / pi, abs(fh$h), col = "blue")
# plot in dB:
plot(f, 20*log10(m+1e-5), type = "b", ylab = "dB", xlab = "Frequency")
lines(fh$w / pi, 20*log10(abs(fh$h)), col = "blue")
par(op)

Least-squares linear-phase FIR filter design

Description

Produce a linear phase filter such that the integral of the weighted mean squared error in the specified bands is minimized.

Usage

firls(n, f, a, w = rep(1L, length(a)/2))

Arguments

Value

The FIR filter coefficients, a vector of length n + 1, of class Ma.

Author(s)

Quentin Spencer, qspencer@ieee.org.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

See Also

Ma, filter, fftfilt, fir1

Examples

freqz(firls(255, c(0, 0.25, 0.3, 1), c(1, 1, 0, 0)))

Flat top window

Description

Return the filter coefficients of a flat top window.

Usage

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

Arguments

Details

The Flat Top window is defined by the function:

f(w) = 1 - 1.93 cos(2 \pi w) + 1.29 cos(4 \pi w) - 0.388 cos(6 \pi w) + 0.0322 cos(8 \pi w)

where w = i/(n-1) for i=0:n-1 for a symmetric window, or w = i/n for i=0:n-1 for a periodic window. The default is symmetric. The returned window is normalized to a peak of 1 at w = 0.5.

Flat top windows have very low passband ripple (< 0.01 dB) and are used primarily for calibration purposes. Their bandwidth is approximately 2.5 times wider than a Hann window.

Value

Flat top window, returned as a vector.

Author(s)

Paul Kienzle, pkienzle@users.sf.net.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

Examples

ft <- flattopwin(64)
plot (ft, type = "l", xlab = "Samples", ylab =" Amplitude")

Fractional shift

Description

Shift a signal by a (possibly fractional) number of samples.

Usage

fracshift(x, d, h = NULL)

Arguments

Details

The function calculates the initial index and end index of the sequences of 1’s in the rows of 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 matrices 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 1s and the second row is the end index of that sequence. If nr == 1 the output is a matrix with two rows.

Author(s)

Eric Chassande-Mottin, ecm@apc.univ-paris7.fr,
Juan Pablo Carbajal, carbajal@ifi.uzh.ch,
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

References

[1] A. V. Oppenheim, R. W. Schafer and J. R. Buck, Discrete-time signal processing, Signal processing series, Prentice-Hall, 1999.
[2] T.I. Laakso, V. Valimaki, M. Karjalainen and U.K. Laine Splitting the unit delay, IEEE Signal Processing Magazine, vol. 13, no. 1, pp 30–59 Jan 1996.

Examples

N = 1024
t <- seq(0, 1, length.out = N)
x <- exp(-t^2 / 2 / 0.25^2) * sin(2 * pi * 10 * t)
dt <- 0.25
d  <- dt / (t[2] - t[1])
y <- fracshift(x, d)
plot(t, x, type = "l", xlab = "Time", ylab = "Sigfnal")
lines (t, y, col = "red")
legend("topright", legend = c("original", "shifted"), lty = 1, col = 1:2)

Frequency response of analog filters

Description

Compute the s-plane frequency response of an IIR filter.

Usage

freqs(filt, ...)

## Default S3 method:
freqs(filt, a, w, ...)

## S3 method for class 'Arma'
freqs(filt, w, ...)

## S3 method for class 'Ma'
freqs(filt, w, ...)

## S3 method for class 'Sos'
freqs(filt, w, ...)

## S3 method for class 'Zpg'
freqs(filt, w, ...)

## S3 method for class 'freqs'
print(x, ...)

## S3 method for class 'freqs'
summary(object, ...)

## S3 method for class 'summary.freqs'
print(x, ...)

freqs_plot(x, ...)

Arguments

Details

The s-plane frequency response of the IIR filter B(s) / A(s) is computed as H = polyval(B, 1i * W) / polyval(A, 1i * W). If called with no output argument, a plot of magnitude and phase are displayed.

Value

For freqs, a list of class 'freqs' with items:

h

complex array of frequency responses at frequencies f.

w

array of frequencies.

Author(s)

Julius O. Smith III, jos@ccrma.stanford.edu.
Conversion to R by Geert van Boxtel gjmvanboxtel@gmail.com

Examples

b <- c(1, 2); a <- c(1, 1)
w <- seq(0, 4, length.out = 128)
freqs (b, a, w)

Frequency response of digital filter

Description

Compute the z-plane frequency response of an ARMA model or rational IIR filter.

Usage

freqz(filt, ...)

## Default S3 method:
freqz(
  filt,
  a = 1,
  n = 512,
  whole = ifelse((is.numeric(filt) && is.numeric(a)), FALSE, TRUE),
  fs = 2 * pi,
  ...
)

## S3 method for class 'Arma'
freqz(
  filt,
  n = 512,
  whole = ifelse((is.numeric(filt$b) && is.numeric(filt$a)), FALSE, TRUE),
  fs = 2 * pi,
  ...
)

## S3 method for class 'Ma'
freqz(
  filt,
  n = 512,
  whole = ifelse(is.numeric(filt), FALSE, TRUE),
  fs = 2 * pi,
  ...
)

## S3 method for class 'Sos'
freqz(filt, n = 512, whole = FALSE, fs = 2 * pi, ...)

## S3 method for class 'Zpg'
freqz(filt, n = 512, whole = FALSE, fs = 2 * pi, ...)

## S3 method for class 'freqz'
print(x, ...)

## S3 method for class 'freqz'
summary(object, ...)

## S3 method for class 'summary.freqz'
print(x, ...)

freqz_plot(w, h, ...)

Arguments

Details

The frequency response of a digital filter can be interpreted as the transfer function evaluated at z = e^{j\omega}.

The 'Matlab' and 'Octave' versions of freqz produce magnitude and phase plots. The freqz version in the 'signal' package produces separate plots of magnitude in the pass band (max - 3 dB to max) and stop (total) bands, as well as a phase plot. The current version produces slightly different plots. The magnitude plots are separate for stop and pass bands, but the pass band plot has an absolute lower limit of -3 dB instead of max - 3 dB. In addition a summary method was added that prints out the most important information about the frequency response of the filter.

Value

For freqz, a list of class 'freqz' with items:

h

complex array of frequency responses at frequencies f.

w

array of frequencies.

u

units of (angular) frequency; either rad/s or Hz.

Note

When results of freqz are printed, freqz_plot will be called to display frequency plots of magnitude and phase. As with lattice plots, automatic printing does not work inside loops and function calls, so explicit calls to print or plot are needed there.

Author(s)

John W. Eaton, Paul Kienzle, pkienzle@users.sf.net.
Port to R by Tom Short,
adapted by Geert van Boxtel, gjmvanboxtel@gmail.com

Examples

b <- c(1, 0, -1)
a <- c(1, 0, 0, 0, 0.25)
freqz(b, a)

hw <- freqz(b, a)
summary(hw)

Full width at half maximum

Description

Compute peak full-width at half maximum or at another level of peak maximum for a vector or matrix.

Usage

fwhm(
  x = seq_len(length(y)),
  y,
  ref = c("max", "zero", "middle", "min", "absolute"),
  level = 0.5
)

Arguments

Value

Full width at half maximum, returned as a vector with a length equal to the number of columns in y, or 1 in case of a vector.

Author(s)

Petr Mikulik.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

Examples

x <- seq(-pi, pi, 0.001)
y <- cos(x)
w <- fwhm(x, y)
m <- x[which.max(y)]
f <- m - w/2
t <- m + w/2
plot(x, y, type="l",
     panel.first = {
       usr <- par('usr')
       rect(f, usr[3], t, usr[4], col = rgb(0, 1, 0, 0.4), border = NA)
     })
abline(h = max(y) / 2, lty = 2, col = "gray")

Gaussian-modulated sinusoidal RF pulse

Description

Generate a Gaussian modulated sinusoidal pulse sampled at times t.

Usage

gauspuls(t, fc = 1000, bw = 0.5)

Arguments

Value

Inphase Gaussian-modulated sinusoidal pulse, returned as a vector of unit amplitude at the times indicated by the time vector t.

Author(s)

Sylvain Pelissier, Mike Miller.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

Examples

fs <- 11025    # arbitrary sample rate
t <- seq(-10, 10, 1/fs)
yi1 <- gauspuls(t, 0.1, 1)
yi2 <- gauspuls(t, 0.1, 2)
plot(t, yi1, type="l", xlab = "Time", ylab = "Amplitude")
lines(t, yi2, col = "red")

fs <- 11025  # arbitrary sample rate
f0 <- 100    # pulse train sample rate
x <- pulstran (seq(0, 4/f0, 1/fs), seq(0, 4/f0, 1/f0), "gauspuls")
plot (0:(length(x)-1) * 1000/fs, x, type="l",
      xlab = "Time (ms)", ylab = "Amplitude",
      main = "Gaussian pulse train at 10 ms intervals")

Gaussian convolution window

Description

Return a Gaussian convolution window of length n.

Usage

gaussian(n, a = 1)

Arguments

Details

The width of the window is inversely proportional to the parameter a. Use larger a for a narrower window. Use larger m for longer tails.

w = e^{(-(a*x)^{2}/2 )}

for x <- seq(-(n - 1) / 2, (n - 1) / 2, by = n).

Width a is measured in frequency units (sample rate/num samples). It should be f when multiplying in the time domain, but 1/f when multiplying in the frequency domain (for use in convolutions).

Value

Gaussian convolution window, returned as a vector.

Author(s)

Paul Kienzle, pkienzle@users.sf.net.
Conversion to R by Geert van Boxtel G.J.M.vanBoxtel@gmail.com.

Examples

g1 <- gaussian(128, 1)
g2 <- gaussian(128, 0.5)
plot (g1, type = "l", xlab = "Samples", ylab =" Amplitude", ylim = c(0, 1))
lines(g2, col = "red")

Gaussian window

Description

Return the filter coefficients of a Gaussian window of length n.

Usage

gausswin(n, a = 2.5)

Arguments

Details

The width of the window is inversely proportional to the parameter a. Use larger a for a narrower window. Use larger m for a smoother curve.

w = e^{(-(a*x)^{2}/2 )}

for x <- seq(-(n - 1) / n, (n - 1) / n, by = n).

The exact correspondence with the standard deviation of a Gaussian probability density function is \sigma = (n - 1) / (2a).

Value

Gaussian convolution window, returned as a vector.

Author(s)

Paul Kienzle, pkienzle@users.sf.net.
Conversion to R by Geert van Boxtel G.J.M.vanBoxtel@gmail.com.

Examples

g1 <- gausswin(64)
g2 <- gausswin(64, 5)
plot (g1, type = "l", xlab = "Samples", ylab =" Amplitude", ylim = c(0, 1))
lines(g2, col = "red")

Gaussian monopulse

Description

Returns samples of the unit-amplitude Gaussian monopulse.

Usage

gmonopuls(t, fc = 1000)

Arguments

Value

Samples of the Gaussian monopulse, returned as a vector of unit amplitude at the times indicated by the time vector t.

Author(s)

Sylvain Pelissier, sylvain.pelissier@gmail.com.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

Examples

fs <- 11025    # arbitrary sample rate
t <- seq(-10, 10, 1/fs)
y1 <- gmonopuls(t, 0.1)
y2 <- gmonopuls(t, 0.2)
plot(t, y1, type="l", xlab = "Time", ylab = "Amplitude")
lines(t, y2, col = "red")
legend("topright", legend = c("fc = 0.1", "fc = 0.2"),
       lty = 1, col = c(1, 2))

Group delay

Description

Compute the average delay of a filter (group delay).

Usage

grpdelay(filt, ...)

## S3 method for class 'grpdelay'
print(x, ...)

## S3 method for class 'grpdelay'
plot(
  x,
  xlab = if (x$HzFlag) "Frequency (Hz)" else "Frequency (rad/sample)",
  ylab = "Group delay (samples)",
  type = "l",
  ...
)

## Default S3 method:
grpdelay(filt, a = 1, n = 512, whole = FALSE, fs = NULL, ...)

## S3 method for class 'Arma'
grpdelay(filt, ...)

## S3 method for class 'Ma'
grpdelay(filt, ...)

## S3 method for class 'Sos'
grpdelay(filt, ...)

## S3 method for class 'Zpg'
grpdelay(filt, ...)

Arguments

Details

If the denominator of the computation becomes too small, the group delay is set to zero. (The group delay approaches infinity when there are poles or zeros very close to the unit circle in the z plane.)

Value

A list of class grpdelay with items:

gd

the group delay, in units of samples. It can be converted to seconds by multiplying by the sampling period (or dividing by the sampling rate fs).

w

frequencies at which the group delay was calculated.

ns

number of points at which the group delay was calculated.

Hzflag

TRUE for frequencies in Hz, FALSE for frequencies in radians.

Author(s)

Paul Kienzle, pkienzle@users.sf.net,
Julius O. Smith III, jos@ccrma.stanford.edu.
Conversion to R by Tom Short,
adapted by Geert van Boxtel, gjmvanboxtel@gmail.com

References

https://ccrma.stanford.edu/~jos/filters/Numerical_Computation_Group_Delay.html
https://en.wikipedia.org/wiki/Group_delay

Examples

# Two Zeros and Two Poles
b <- poly(c(1 / 0.9 * exp(1i * pi * 0.2), 0.9 * exp(1i * pi * 0.6)))
a <- poly(c(0.9 * exp(-1i * pi * 0.6), 1 / 0.9 * exp(-1i * pi * 0.2)))
gpd <- grpdelay(b, a, 512, whole = TRUE, fs = 1)
print(gpd)
plot(gpd)

Hamming window

Description

Return the filter coefficients of a Hamming window of length n.

Usage

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

Arguments

Details

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

Value

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

Author(s)

Andreas Weingessel, Andreas.Weingessel@ci.tuwien.ac.at.
Conversion to R by Geert van Boxtel G.J.M.vanBoxtel@gmail.com.

Examples

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

hs = hamming(64,'symmetric')
hp = hamming(63,'periodic')
plot (hs, type = "l", xlab = "Samples", ylab =" Amplitude")
lines(hp, col="red")

Hann window

Description

Return the filter coefficients of a Hann window of length n.

Usage

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

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

Arguments

Details

The Hann window is a member of the family of cosine sum windows. It was named after Julius von Hann, and is sometimes referred to as Hanning, presumably due to its linguistic and formulaic similarities to Hamming window.

Value

Hann window, returned as a vector.

Author(s)

Andreas Weingessel, Andreas.Weingessel@ci.tuwien.ac.at.
Conversion to R by Geert van Boxtel G.J.M.vanBoxtel@gmail.com.

Examples

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

hs = hann(64,'symmetric')
hp = hann(63,'periodic')
plot (hs, type = "l", xlab = "Samples", ylab =" Amplitude")
lines(hp, col="red")

Hilbert transform

Description

Computes the extension of a real valued signal to an analytic signal.

Usage

hilbert(x, n = ifelse(is.vector(x), length(x), nrow(x)))

Arguments

Details

The function returns returns a complex helical sequence, sometimes called the analytic signal, from a real data sequence. The analytic signal has a real part, which is the original data, and an imaginary part, which contains the Hilbert transform. The imaginary part is a version of the original real sequence with a 90 degrees phase shift. Sines are therefore transformed to cosines, and conversely, cosines are transformed to sines. The Hilbert-transformed series has the same amplitude and frequency content as the original sequence. The transform includes phase information that depends on the phase of the original.

Value

Analytic signal, of length n, returned as a complex vector or matrix, the real part of which contains the original signal, and the imaginary part of which contains the Hilbert transform of x.

Author(s)

Paul Kienzle, pkienzle@users.sf.net,
Peter L. Soendergaard.
Conversion to R by Geert van Boxtel, gjmvanboxtel@gmail.com

References

https://en.wikipedia.org/wiki/Hilbert_transform, https://en.wikipedia.org/wiki/Analytic_signal

Examples

## notice that the imaginary signal is phase-shifted 90 degrees
t <- seq(0, 10, length = 256)
z <- hilbert(sin(2 * pi * 0.5 * t))
plot(t, Re(z), type = "l", col="blue")
lines (t, Im(z), col = "red")
legend('topright', lty = 1, legend = c("Real", "Imag"),
       col = c("blue", "red"))

## the magnitude of the hilbert transform eliminates the carrier
t <- seq(0, 10, length = 1024)
x <- 5 * cos(0.2 * t) * sin(100 * t)
plot(t, x, type = "l", col = "green")
lines (t, abs(hilbert(x)), col = "blue")
legend('topright', lty = 1, legend = c("x", "|hilbert(x)|"),
        col = c("green", "blue"))

Inverse Discrete Cosine Transform

Description

Compute the inverse unitary discrete cosine transform of a signal.

Usage

idct(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.

Value

Inverse discrete cosine transform, returned as a vector or matrix.

Author(s)

Paul Kienzle, pkienzle@users.sf.net.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

See Also

dct

Examples

x <- 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)

Inverse 2-D Discrete Cosine Transform

Description

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

Usage

idct2(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, pkienzle@users.sf.net.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

See Also

dct2

Examples

A <- matrix(50 * runif(100), 10, 10)
B <- dct2(A)
B[which(B < 1)] <- 0
AA <- idct2(B)

Inverse Discrete Sine Transform

Description

Compute the inverse discrete sine transform of a signal.

Usage

idst(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.

Value

Inverse discrete sine transform, returned as a vector or matrix.

Author(s)

Paul Kienzle, pkienzle@users.sf.net.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

See Also

dst

Examples

x <- seq_len(100) + 50 * cos(seq_len(100) * 2 * pi / 40)
X <- dst(x)
xx <- idst(X)
all.equal(x, xx)

Inverse Fast Fourier Transform

Description

Compute the inverse Fast Fourier Transform compatible with 'Matlab' and 'Octave'.

Usage

ifft(x)

imvfft(x)

Arguments

Details

The 'fft' function in the 'stats' package can compute the inverse FFT by specifying inverse = TRUE. However, that function does not divide the result by length(x), nor does it return real values when appropriate. The present function does both, and is this compatible with 'Matlab' and 'Octave' (and differs from the 'ifft' function in the 'signal' package, which does not return real values).

Value

When x is a vector, the value computed and returned by ifft is the univariate inverse discrete Fourier transform of the sequence of values in x. Specifically, y <- ifft(x) is defined as stats::fft(x, inverse = TRUE) / length(x). The stats::fft function called with inverse = TRUE replaces exp(-2 * pi...) with exp(2 * pi) in the definition of the discrete Fourier transform (see fft).

When x contains an array, ifft computes and returns the normalized inverse multivariate (spatial) transform. By contrast, imvfft takes a real or complex matrix as argument, and returns a similar shaped matrix, but with each column replaced by its normalized inverse discrete Fourier transform. This is useful for analyzing vector-valued series.

Author(s)

Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

See Also

fft

Examples

res <- ifft(stats::fft(1:5))
res <- ifft(stats::fft(c(1+5i, 2+3i, 3+2i, 4+6i, 5+2i)))
res <- imvfft(stats::mvfft(matrix(1:20, 4, 5)))

Inverse zero-frequency shift

Description

Rearranges a zero-frequency-shifted Fourier transform back to the original.

Usage

ifftshift(x, MARGIN = 2)

Arguments

Details

Undo the action of the fftshift function. For even length x, fftshift is its own inverse, but not for odd length input.

Value

back-transformed vector or matrix.

Author(s)

Vincent Cautaerts, vincent@comf5.comm.eng.osaka-u.ac.jp,
adapted by John W. Eaton.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

See Also

fftshift

Examples

Xeven <- 1:6
res <- fftshift(fftshift(Xeven))

Xodd <- 1:7
res <- fftshift(fftshift(Xodd))
res <- ifftshift(fftshift(Xodd))

Fast Walsh-Hadamard Transform

Description

Compute the (inverse) Fast Walsh-Hadamard transform of a signal.

Usage

ifwht(x, n = NROW(x), ordering = c("sequency", "hadamard", "dyadic"))

fwht(x, n = NROW(x), ordering = c("sequency", "hadamard", "dyadic"))

Arguments

Value

(Inverse) Fast Walsh Hadamard transform, returned as a vector or matrix.

Author(s)

Mike Miller.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

References

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

https://en.wikipedia.org/wiki/Fast_Walsh-Hadamard_transform

Examples

x <- c(19, -1, 11, -9, -7, 13, -15, 5)
X <- fwht(x)
all.equal(x, ifwht(X))

IIR lowpass filter to IIR multiband

Description

Transform an IIR lowpass filter prototype to an IIR multiband filter.

Usage

iirlp2mb(b, ...)

## S3 method for class 'Arma'
iirlp2mb(b, Wo, Wt, type, ...)

## S3 method for class 'Zpg'
iirlp2mb(b, Wo, Wt, type, ...)

## S3 method for class 'Sos'
iirlp2mb(b, Wo, Wt, type, ...)

## Default S3 method:
iirlp2mb(b, a, Wo, Wt, type = c("pass", "stop"), ...)

Arguments

Details

The utility of a prototype filter comes from the property that all other filters can be derived from it by applying a scaling factor to the components of the prototype. The filter design need thus only be carried out once in full, with other filters being obtained by simply applying a scaling factor. Especially useful is the ability to transform from one bandform to another. In this case, the transform is more than a simple scale factor. Bandform here is meant to indicate the category of passband that the filter possesses. The usual bandforms are lowpass, highpass, bandpass and bandstop, but others are possible. In particular, it is possible for a filter to have multiple passbands. In fact, in some treatments, the bandstop filter is considered to be a type of multiple passband filter having two passbands. Most commonly, the prototype filter is expressed as a lowpass filter, but other techniques are possible[1].

Filters with multiple passbands may be obtained by applying the general transformation described in [2].

Because iirlp2mb is generic, it can be extended to accept other inputs.

Value

List of class Arma numerator and denominator polynomials of the resulting filter.

Author(s)

Alan J. Greenberger, alanjg@ptd.net.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

References

[1] https://en.wikipedia.org/wiki/Prototype_filter
[2] https://en.wikipedia.org/wiki/Prototype_filter#Lowpass_to_multi-band

Examples

## Design a prototype real IIR lowpass elliptic filter with a gain of about
## –3 dB at 0.5pi rad/sample.
el <- ellip(3, 0.1, 30, 0.409)
## Create a real multiband filter with two passbands.
mb1 <- iirlp2mb(el, 0.5, c(.2, .4, .6, .8), 'pass')
## Create a real multiband filter with two stopbands.
mb2 <- iirlp2mb(el, 0.5, c(.2, .4, .6, .8), 'stop')
## Compare the magnitude responses of the filters.
hfl <- freqz(el)
hf1 <- freqz(mb1)
hf2 <- freqz(mb2)
plot(hfl$w, 20 * log10(abs(hfl$h)), type = "l",
    xlab = "Normalized frequency (* pi rad/sample)",
    ylab = "Magnitude (dB)")
lines(hf1$w, 20 * log10(abs(hf1$h)), col="red")
lines(hf2$w, 20 * log10(abs(hf2$h)), col="blue")
legend('bottomleft',
       legend = c('Prototype', 'Two passbands', 'Two Stopbands'),
       col=c("black", "red", "blue"), lty = 1)

Impulse invariance method for A/D filter conversion

Description

Convert analog filter with coefficients b and a to digital, conserving impulse response.

Usage

impinvar(b, ...)

## S3 method for class 'Arma'
impinvar(b, ...)

## Default S3 method:
impinvar(b, a, fs = 1, tol = 1e-04, ...)

Arguments

Details

Because impinvar is generic, it can also accept input of class Arma.

Value

A list of class Arma containing numerator and denominator polynomial filter coefficients of the A/D converted filter.

Author(s)

Tony Richardson, arichard@stark.cc.oh.us,
Ben Abbott, bpabbott@mac.com,
adapted by John W. Eaton.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com

See Also

invimpinvar

Examples

f <- 2
fs <- 10
but <- butter(6, 2 * pi * f, 'low', 's')
zbut <- impinvar(but, fs)
freqz(zbut, n = 1024, fs = fs)

Impulse response of digital filter

Description

Compute the z-plane impulse response of an ARMA model or rational IIR filter. A plot of the impulse and step responses is generated.

Usage

impz(filt, ...)

## S3 method for class 'impz'
print(x, ...)

## S3 method for class 'Arma'
impz(filt, ...)

## S3 method for class 'Ma'
impz(filt, ...)

## S3 method for class 'Sos'
impz(filt, ...)

## S3 method for class 'Zpg'
impz(filt, ...)

## Default S3 method:
impz(filt, a = 1, n = NULL, fs = 1, ...)

Arguments

Value

For impz, a list of class "impz" with items:

x

impulse response signal.

t

time.

Note

When results of impz are printed, plot will be called to display a plot of the impulse response against frequency. As with lattice plots, automatic printing does not work inside loops and function calls, so explicit calls to print or plot are needed there.

Author(s)

Paul Kienzle, pkienzle@users.sf.net.
Conversion to R by Tom Short;
adapted by Geert van Boxtel, gjmvanboxtel@gmail.com

Examples

## elliptic low-pass filter
elp <- ellip(4, 0.5, 20, 0.4)
impz(elp)

xt <- impz(elp)

Interpolation

Description

Increase sample rate by integer factor.

Usage

interp(x, q, n = 4, Wc = 0.5)

Arguments

Value

interpolated signal, returned as a vector.

Author(s)

Paul Kienzle, pkienzle@users.sf.net.
Conversion to R by Geert van Boxtel G.J.M.vanBoxtel@gmail.com.

See Also

decimate, resample

Examples

# Generate a signal
t <- seq(0, 2, 0.01)
x <- chirp(t, 2, .5, 10,'quadratic') + sin(2 * pi * t * 0.4)
w <- seq(1, 121, 4)
plot(t[w] * 1000, x[w], type = "h", xlab = "", ylab = "")
points(t[w] * 1000, x[w])
abline (h = 0)
y <- interp(x[seq(1, length(x), 4)], 4, 4, 1)
lines(t[1:121] * 1000, y[1:121], type = "l", col = "red")
points(t[1:121] * 1000, y[1:121], col = "red", pch = '+')
legend("topleft", legend = c("original", "interpolated"),
  lty = 1, pch = c(1, 3), col = c(1, 2))

Inverse Frequency Response

Description

Identify filter parameters from frequency response data.

Usage

invfreq(
  h,
  w,
  nb,
  na,
  wt = rep(1, length(w)),
  plane = c("z", "s"),
  method = c("ols", "tls", "qr"),
  norm = TRUE
)

invfreqs(
  h,
  w,
  nb,
  na,
  wt = rep(1, length(w)),
  method = c("ols", "tls", "qr"),
  norm = TRUE
)

invfreqz(
  h,
  w,
  nb,
  na,
  wt = rep(1, length(w)),
  method = c("ols", "tls", "qr"),
  norm = TRUE
)

Arguments

Details

Given a desired (one-sided, complex) spectrum h(w) at equally spaced angular frequencies w = (2 \pi k) / N, k = 0, ... N-1, this function finds the filter B(z)/A(z) or B(s)/A(s) with nb zeroes and na poles. Optionally, the fit-errors can be weighted with respect to frequency according to the weights wt.

Value

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

b

moving average (MA) polynomial coefficients

a

autoregressive (AR) polynomial coefficients

Author(s)

Julius O. Smith III, Rolf Schirmacher, Andrew Fitting, Pascal Dupuis.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

References

https://ccrma.stanford.edu/~jos/filters/FFT_Based_Equation_Error_Method.html

Examples

order <- 6  # order of test filter
fc <- 1/2   # sampling rate / 4
n <- 128    # frequency grid size
ba <- butter(order, fc)
hw <- freqz(ba, n)
BA = invfreq(hw$h, hw$w, order, order)
HW = freqz(BA, n)
plot(hw$w, abs(hw$h), type = "l", xlab = "Frequency (rad/sample)",
  ylab = "Magnitude")
lines(HW$w, abs(HW$h), col = "red")
legend("topright", legend = c("Original", "Measured"), lty  = 1, col = 1:2)
err <- norm(hw$h - HW$h, type = "2")
title(paste('L2 norm of frequency response error =', err))

Inverse impulse invariance method

Description

Convert digital filter with coefficients b and a to analog, conserving impulse response.

Usage

invimpinvar(b, ...)

## S3 method for class 'Arma'
invimpinvar(b, ...)

## Default S3 method:
invimpinvar(b, a, fs = 1, tol = 1e-04, ...)

Arguments

Details

Because invimpinvar is generic, it can also accept input of class Arma.

Value

A list of class Arma containing numerator and denominator polynomial filter coefficients of the A/D converted filter.

Author(s)

R.G.H. Eschauzier, reschauzier@yahoo.com,
Carne Draug, carandraug+dev@gmail.com.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com

References

Thomas J. Cavicchi (1996) Impulse invariance and multiple-order poles. IEEE transactions on signal processing, Vol 40 (9): 2344–2347.

See Also

impinvar

Examples

f <- 2
fs <- 10
but <- butter(6, 2 * pi * f, 'low', 's')
zbut <- impinvar(but, fs)
sbut <- invimpinvar(zbut, fs)
all.equal(but, sbut, tolerance = 1e-7)

Kaiser window

Description

Return the filter coefficients of a kaiser window of length n.

Usage

kaiser(n, beta = 0.5)

Arguments

Details

The Kaiser, or Kaiser-Bessel, window is a simple approximation of the DPSS window using Bessel functions, discovered by James Kaiser.

        besselI(0, Beta * sqrt(1-(2*x/m)^2))
k(x) =  -------------------------------------,  -m/2 <= x <= m/2
        besselO(0, Beta)

The variable parameter \beta determines the trade-off between main lobe width and side lobe levels of the spectral leakage pattern. Increasing \beta widens the main lobe and decreases the amplitude of the side lobes (i.e., increases the attenuation).

Value

Kaiser window, returned as a vector.

Author(s)

Kurt Hornik, Kurt.Hornik@ci.tuwien.ac.at,
Paul Kienzle, pkienzle@users.sf.net.
Conversion to R by Geert van Boxtel G.J.M.vanBoxtel@gmail.com.

Examples

k <- kaiser(200, 2.5)
plot (k, type = "l", xlab = "Samples", ylab =" Amplitude")

Kaiser filter order and cutoff frequency

Description

Return the parameters needed to produce a FIR filter of the desired specification from a Kaiser window.

Usage

kaiserord(f, m, dev, fs = 2)

Arguments

Details

Given a set of specifications in the frequency domain, kaiserord estimates the minimum FIR filter order that will approximately meet the specifications. kaiserord converts the given filter specifications into passband and stopband ripples and converts cutoff frequencies into the form needed for windowed FIR filter design.

kaiserord uses empirically derived formulas for estimating the orders of lowpass filters, as well as differentiators and Hilbert transformers. Estimates for multiband filters (such as band-pass filters) are derived from the low-pass design formulas.

The design formulas that underlie the Kaiser window and its application to FIR filter design are

\beta =

0.1102(\alpha - 8.7), \alpha > 50

0.5842(\alpha -21)^{0.4} + 0.07886(\alpha - 21), 21 \le \alpha \le 50

0, \alpha < 21

where \alpha = -20log_{10}(\delta) is the stopband attenuation expressed in decibels, n=(\alpha - 8) / 2.285(\Delta\omega), where n is the filter order and \Delta\omega is the width of the smallest transition region.

Value

A list of class FilterSpecs with the following list elements:

n

filter order

Wc

cutoff frequency

type

filter type, one of "low", "high", "stop", "pass", "DC-0", or "DC-1".

beta

shape parameter

Author(s)

Paul Kienzle.
Conversion to R by Tom Short,
adapted by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

See Also

hamming, kaiser

Examples

fs <- 11025
op <- par(mfrow = c(2, 2), mar = c(3, 3, 1, 1))
for (i in 1:4) {
  if (i == 1) {
    bands <- c(1200, 1500)
    mag <- c(1, 0)
    dev <- c(0.1, 0.1)
  }
  if (i == 2) {
    bands <- c(1000, 1500)
    mag <- c(0, 1)
    dev <- c(0.1, 0.1)
  }
  if (i == 3) {
    bands <- c(1000, 1200, 3000, 3500)
    mag <- c(0, 1, 0)
    dev <- 0.1
  }
  if (i == 4) {
    bands <- 100 * c(10, 13, 15, 20, 30, 33, 35, 40)
    mag <- c(1, 0, 1, 0, 1)
    dev <- 0.05
  }
  kaisprm <- kaiserord(bands, mag, dev, fs)
  d <- max(1, trunc(kaisprm$n / 10))
  if (mag[length(mag)] == 1 && (d %% 2) == 1) {
     d <- d + 1
  }
  f1 <- freqz(fir1(kaisprm$n, kaisprm$Wc, kaisprm$type,
                   kaiser(kaisprm$n + 1, kaisprm$beta),
                   scale = FALSE),
              fs = fs)
  f2 <- freqz(fir1(kaisprm$n - d, kaisprm$Wc, kaisprm$type,
                   kaiser(kaisprm$n - d + 1, kaisprm$beta),
                   scale = FALSE),
              fs = fs)
  plot(f1$w, abs(f1$h), col = "blue", type = "l",  xlab = "", ylab = "")
  lines(f2$w, abs(f2$h), col = "red")
  legend("right", paste("order", c(kaisprm$n-d, kaisprm$n)),
         col = c("red", "blue"), lty = 1, bty = "n")
  b <- c(0, bands, fs/2)
  for (i in seq(2, length(b), by=2)) {
    hi <- mag[i/2] + dev[1]
    lo <- max(mag[i/2] - dev[1], 0)
    lines(c(b[i-1], b[i], b[i], b[i-1], b[i-1]), c(hi, hi, lo, lo, hi))
  }
}
par(op)

Durbin-Levinson Recursion

Description

Use the Durbin-Levinson algorithm to compute the coefficients of an autoregressive linear process.

Usage

levinson(acf, p = NROW(acf))

Arguments

Details

levinson uses the Durbin-Levinson algorithm to solve:

toeplitz(acf(1:p)) * x = -acf(2:p+1)

The solution c(1, x) is the denominator of an all pole filter approximation to the signal x which generated the autocorrelation function acf.

From ref [2]: Levinson recursion or Levinson–Durbin recursion is a procedure in linear algebra to recursively calculate the solution to an equation involving a Toeplitz matrix. Other methods to process data include Schur decomposition and Cholesky decomposition. In comparison to these, Levinson recursion (particularly split Levinson recursion) tends to be faster computationally, but more sensitive to computational inaccuracies like round-off errors.

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.

Author(s)

Paul Kienzle, pkienzle@users.sf.net,
Peter V. Lanspeary, pvl@mecheng.adelaide.edu.au.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

References

[1] Steven M. Kay and Stanley Lawrence Marple Jr. (1981). Spectrum analysis – a modern perspective. Proceedings of the IEEE, Vol 69, 1380-1419.
[2] https://en.wikipedia.org/wiki/Levinson_recursion

Examples

## Estimate the coefficients of an autoregressive process given by
## x(n) = 0.1x(n-1) - 0.8x(n-2) - 0.27x(n-3) + w(n).
a <- c(1, 0.1, -0.8, -0.27)
v <- 0.4
w <- sqrt(v) * rnorm(15000)
x <- filter(1, a, w)
xc <- xcorr(x, scale = 'biased')
acf <- xc$R[-which(xc$lags < 0)]
lev <- levinson(acf, length(a) - 1)

Marcum Q function

Description

Compute the generalized Marcum Q function

Usage

marcumq(a, b, m = 1)

Arguments

Details

The code for this function was taken from the help file of the cdfkmu function in the lmomco package, based on a suggestion of Daniel Wollschlaeger.

Value

Marcum Q function.

Author(s)

William Asquith, william.asquith@ttu.edu.

References

https://cran.r-project.org/package=lmomco

Examples

mq <- marcumq(12.4, 12.5)

1-D median filtering

Description

Apply a running median of odd span to the input x

Usage

medfilt1(x, n = 3, MARGIN = 2, na.omit = FALSE, ...)

Arguments

Details

This function computes a running median over the input x, using the runmed function. Because of that, it works a little differently than the 'Matlab' or 'Octave' versions (i.e., it does not produce exactly the same values).

missing values

The 'Mablab' and 'Octave' functions have a 'nanflag' option that allows to include or remove missing values. If inclusion is specifies, then the function returns a signal so that the median of any segment containing NAs is also NA. Because the 'runmed' function does not include an na.omit option, implementing this functionality would lead to a considerable speed loss. Instead, a na.omit parameter was implemented that allows either omitting NAs or interpolating them with a spline function.

endpoint filtering

Instead of the 'zeropad' and 'truncate' options to the 'padding' argument in the 'Matlab' and 'Octave' functions, the present version uses the standard endrule parameter of the 'runmed' function, with options keep, constant, or median.

Value

Filtered signal, returned as a numeric vector, matrix, or array, of the same size as x.

Author(s)

Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

See Also

runmed, splinefun

Examples

## noise suppression
fs <- 100
t <- seq(0, 1, 1/fs)
x <- sin(2 * pi * t * 3) + 0.25 * sin(2 * pi * t * 40)
plot(t, x, type = "l", xlab = "", ylab = "")
y <- medfilt1(x, 11)
lines (t, y, col = "red")
legend("topright", c("Original", "Filtered"), lty = 1, col = 1:2)

Mexicat Hat

Description

Generate a Mexican Hat (Ricker) wavelet sampled on a regular grid.

Usage

mexihat(lb = -5, ub = 5, n = 1000)

Arguments

Details

The Mexican Hat or Ricker wavelet is the negative normalized second derivative of a Gaussian function, i.e., up to scale and normalization, the second Hermite function. It is a special case of the family of continuous wavelets (wavelets used in a continuous wavelet transform) known as Hermitian wavelets. The Ricker wavelet is frequently employed to model seismic data, and as a broad spectrum source term in computational electrodynamics. It is usually only referred to as the Mexican hat wavelet in the Americas, due to taking the shape of a sombrero when used as a 2D image processing kernel. It is also known as the Marr wavelet (source: Wikipedia)

Value

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

Author(s)

Sylvain Pelissier, sylvain.pelissier@gmail.com.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

Examples

mh <- mexihat(-5, 5, 1000)
plot(mh$x, mh$psi, type="l", main = "Mexican Hat Wavelet",
     xlab = "", ylab = "")

Meyer wavelet auxiliary function

Description

Compute the Meyer wavelet auxiliary function.

Usage

meyeraux(x)

Arguments

Details

The code y = meyeraux(x) returns values of the auxiliary function used for Meyer wavelet generation evaluated at the elements of x. The input x is a vector or matrix of real values. The function is

y = 35x^{4} - 84x^{5} + 70x^{6} - 20x^{7}.

x and y have the same dimensions. The range of meyeraux is the closed interval c(0, 1).

Value

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

Author(s)

Sylvain Pelissier, sylvain.pelissier@gmail.com.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

Examples

x <- seq(0, 1, length.out = 100)
y <- meyeraux(x)
plot(x, y, type="l", main = "Meyer wavelet auxiliary function",
     xlab = "", ylab = "")

Morlet Wavelet

Description

Compute the Morlet wavelet on a regular grid.

Usage

morlet(lb = -4, ub = 4, n = 1000)

Arguments

Details

The code m <- morlet(lb, ub, n) returns values of the Morlet wavelet on an n-point regular grid in the interval c(lb, ub).

The Morlet waveform is defined as

\psi(x) = e^{-x^{2}/2} cos (5x)

Value

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

Author(s)

Sylvain Pelissier, sylvain.pelissier@gmail.com.
Conversion to R by Geert van Boxtel G.J.M.vanBoxtel@gmail.com.

Examples

m <- morlet(-4, 4, 1000)
plot(m$x, m$psi, type="l", main = "Morlet Wavelet", xlab = "", ylab = "")

Moving Root Mean Square

Description

Compute the moving root mean square (RMS) of the input signal.

Usage

movingrms(x, width = 0.1, rc = 0.001, fs = 1)

Arguments

Details

The signal is convoluted against a sigmoid window of width w and risetime rc. The units of these parameters are relative to the value of the sampling frequency given in fs.

Value

A list containing 2 variables:

rmsx

Output signal with the same dimensions as x

w

Window, returned as a vector

Author(s)

Juan Pablo Carbajal, carbajal@ifi.uzh.ch.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

See Also

sigmoid_train

Examples

       
N <- 128
fs <- 5
t <- seq(0, 1, length.out = N)
x <- sin(2 * pi * fs * t) + runif(N)
y <- movingrms(x, 5)

Multiplicity of poles

Description

Identify unique poles and their associated multiplicity.

Usage

mpoles(p, tol = 0.001, reorder = TRUE, index.return = FALSE)

Arguments

Value

If index.return = TRUE, a list consisting of two vectors:

m

vector specifying the multiplicity of the poles

n

index

If index.return = FALSE, only m is returned (as a vector).

Author(s)

Ben Abbott, bpabbott@mac.com.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com

See Also

poly, residue

Examples

p <- c(2, 3, 1, 1, 2)
ret <- mpoles(p, index = TRUE)

Magnitude-squared coherence

Description

Compute the magnitude-squared coherence estimates of input signals.

Usage

mscohere(
  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")
)

cohere(
  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

mscohere estimates the magnitude-squared coherence 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.

coh

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 mscohere (and its deprecated alias cohere) is a wrapper for the function pwelch, which is more complete and more flexible.

Author(s)

Peter V. Lanspeary, pvl@mecheng.adelaide.edu.au.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

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 <- mscohere(cbind(s1, s2), fs = fs)
plot(rv$freq, rv$coh, type="l", xlab = "Frequency", ylab = "Coherence")

ncauer analog filter design

Description

Compute the transfer function coefficients of a Cauer analog filter.

Usage

ncauer(Rp, Rs, n)

Arguments

Details

Cauer filters have equal maximum ripple in the passband and the stopband. The Cauer filter has a faster transition from the passband to the stopband than any other class of network synthesis filter. The term Cauer filter can be used interchangeably with elliptical filter, but the general case of elliptical filters can have unequal ripples in the passband and stopband. An elliptical filter in the limit of zero ripple in the passband is identical to a Chebyshev Type 2 filter. An elliptical filter in the limit of zero ripple in the stopband is identical to a Chebyshev Type 1 filter. An elliptical filter in the limit of zero ripple in both passbands is identical to a Butterworth filter. The filter is named after Wilhelm Cauer and the transfer function is based on elliptic rational functions.Cauer-type filters use generalized continued fractions.[1]

Value

A list of class Zpg with the following list elements:

zero

complex vector of the zeros of the model

pole

complex vector of the poles of the model

gain

gain of the model

Author(s)

Paulo Neis, p_neis@yahoo.com.br.
Conversion to R Tom Short,
adapted by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

References

[1] https://en.wikipedia.org/wiki/Network_synthesis_filters#Cauer_filter

See Also

Zpg, filter, ellip

Examples

zpg <- ncauer(1, 40, 5)
freqz(zpg)
zplane(zpg)

Nuttall-defined minimum 4-term Blackman-Harris window

Description

Return the filter coefficients of a Blackman-Harris window defined by Nuttall of length n.

Usage

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

Arguments

Details

The window is minimum in the sense that its maximum sidelobes are minimized. The coefficients for this window differ from the Blackman-Harris window coefficients computed with blackmanharris and produce slightly lower sidelobes.

Value

Nuttall-defined Blackman-Harris window, returned as a vector.

Author(s)

Sylvain Pelissier, sylvain.pelissier@gmail.com.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

See Also

blackman, blackmanharris

Examples

n <- nuttallwin(64)
plot (n, type = "l", xlab = "Samples", ylab =" Amplitude")

Pad data

Description

Pre- or postpad the data object x with the value c until it is of length l.

Usage

pad(x, l, c = 0, MARGIN = 2, direction = c("both", "pre", "post"))

prepad(x, l, c = 0, MARGIN = 2)

postpad(x, l, c = 0, MARGIN = 2)

Arguments

Value

Padded data, returned as a vector or matrix.

Author(s)

Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

Examples

v <- 1:24
res <- postpad(v, 30)
res <- postpad(v, 20)
res <- prepad(v, 30)
res <- prepad(v, 20)

m <- matrix(1:24, 4, 6)
res <- postpad(m, 8, 100)
res <- postpad(m, 8, 100, MARGIN = 1)
res <- prepad(m, 8, 100)
res <- prepad(m, 8, 100, MARGIN = 1)

res <- postpad(m, 2)
res <- postpad(m, 2, MARGIN = 1)
res <- prepad(m, 2)
res <- prepad(m, 2, MARGIN = 1)

Parzen (de la Vallée Poussin) window

Description

Return the filter coefficients of a Parzen window of length n.

Usage

parzenwin(n)

Arguments

Details

Parzen windows are piecewise-cubic approximations of Gaussian windows.

Value

Parzen window, returned as a vector.

Author(s)

Sylvain Pelissier, sylvain.pelissier@gmail.com.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

Examples

p <- parzenwin(64)
g <- gausswin(64)
plot (p, type = "l", xlab = "Samples", ylab =" Amplitude", ylim = c(0, 1))
lines(g, col = "red")

Autoregressive PSD estimate - Burg's method

Description

Calculate Burg maximum-entropy power spectral density.

Usage

pburg(
  x,
  p,
  criterion = NULL,
  freq = 256,
  fs = 1,
  range = NULL,
  method = if (length(freq) == 1 && bitwAnd(freq, freq - 1) == 0) "fft" else "poly"
)

Arguments

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.

Note

This function is a wrapper for arburg and ar_psd.

Author(s)

Peter V. Lanspeary, pvl@mecheng.adelaide.edu.au.
Conversion to R by Geert van Boxtel, gjmvanboxtel@gmail.com

See Also

ar_psd, arburg

Examples

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

Maximum-to-minimum difference

Description

Compute the maximum-to-minimum difference of the input data x.

Usage

peak2peak(x, MARGIN = 2)

Arguments

Details

The input x can be a vector, a matrix or an array. If the input is a vector, a single value is returned representing the maximum-to-minimum difference of the vector. If the input is a matrix or an array, a vector or an array of values is returned representing the maximum-to-minimum differences of the dimensions of x indicated by the MARGIN argument.

Support for complex valued input is provided. In this case, the function peak2peak identifies the maximum and minimum in complex magnitude, and then subtracts the complex number with the minimum modulus from the complex number with the maximum modulus.

Value

Vector or array of values containing the maximum-to-minimum differences of the specified MARGIN of x.

Author(s)

Georgios Ouzounis, ouzounis_georgios@hotmail.com.
Conversion to R by Geert van Boxtel G.J.M.vanBoxtel@gmail.com.

Examples

## numeric vector
x <- c(1:5)
pp <- peak2peak(x)

## numeric matrix
x <- matrix(c(1,2,3, 100, 150, 200, 1000, 1500, 2000), 3, 3)
pp <- peak2peak(x)
pp <- peak2peak(x, 1)

## numeric array
x <- array(c(1, 1.5, 2, 100, 150, 200, 1000, 1500, 2000,
             10000, 15000, 20000), c(2,3,2))
pp <- peak2peak(x, 1)
pp <- peak2peak(x, 2)
pp <- peak2peak(x, 3)

## complex input
x <- c(1+1i, 2+3i, 3+5i, 4+7i, 5+9i)
pp <- peak2peak(x)

Peak-magnitude-to-RMS ratio

Description

Compute the ratio of the largest absolute value to the root-mean-square (RMS) value of the object x.

Usage

peak2rms(x, MARGIN = 2)

Arguments

Details

The input x can be a vector, a matrix or an array. If the input is a vector, a single value is returned representing the peak-magnitude-to-RMS ratio of the vector. If the input is a matrix or an array, a vector or an array of values is returned representing the peak-magnitude-to-RMS ratios of the dimensions of x indicated by the MARGIN argument.

Support for complex valued input is provided.

Value

Vector or array of values containing the peak-magnitude-to-RMS ratios of the specified MARGIN of x.

Author(s)

Andreas Weber, octave@tech-chat.de.
Conversion to R by Geert van Boxtel G.J.M.vanBoxtel@gmail.com.

Examples

## numeric vector
x <- c(1:5)
p <- peak2rms(x)

## numeric matrix
x <- matrix(c(1,2,3, 100, 150, 200, 1000, 1500, 2000), 3, 3)
p <- peak2rms(x)
p <- peak2rms(x, 1)

## numeric array
x <- array(c(1, 1.5, 2, 100, 150, 200, 1000, 1500, 2000,
           10000, 15000, 20000), c(2,3,2))
p <- peak2rms(x, 1)
p <- peak2rms(x, 2)
p <- peak2rms(x, 3)

## complex input
x <- c(1+1i, 2+3i, 3+5i, 4+7i, 5+9i)
p <- peak2rms(x)

Pei-Tseng notch filter

Description

Compute the transfer function coefficients of an IIR narrow-band notch filter.

Usage

pei_tseng_notch(w, bw)

Arguments

Details

The filter construction is based on an all-pass which performs a reversal of phase at the filter frequencies. Thus, the mean of the phase-distorted and the original signal has the respective frequencies removed.

Value

List of class Arma with list elements:

b

moving average (MA) polynomial coefficients

a

autoregressive (AR) polynomial coefficients

Author(s)

Alexander Klein, alexander.klein@math.uni-giessen.de.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

References

Pei, Soo-Chang, and Tseng, Chien-Cheng "IIR Multiple Notch Filter Design Based on Allpass Filter"; 1996 IEEE Tencon, doi: doi:10.1109/TENCON.1996.608814

See Also

Arma, filter

Examples

## 50 Hz notch filter
fs <- 256
nyq <- fs / 2
notch <- pei_tseng_notch(50 / nyq, 2 / nyq)
freqz(notch, fs = fs)

Polynomial with specified roots

Description

Compute the coefficients of a polynomial when the roots are given, or the characteristic polynomial of a matrix.

Usage

poly(x)

Arguments

Details

If a vector is passed as an argument, then poly(x) is a vector of the coefficients of the polynomial whose roots are the elements of x.

If an N x N square matrix is given, poly(x) is the row vector of the coefficients of det (z * diag (N) - x), which is the characteristic polynomial of x.

Value

A vector of the coefficients of the polynomial in order from highest to lowest polynomial power.

Author(s)

Kurt Hornik.
Conversion to R by Tom Short,
adapted by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com

See Also

roots

Examples

p <- poly(c(1, -1))
p <- poly(pracma::roots(1:3))
p <- poly(matrix(1:9, 3, 3))

Reduce polynomial

Description

Reduce a polynomial coefficient vector to a minimum number of terms by stripping off any leading zeros.

Usage

polyreduce(pc)

Arguments

Value

Vector of reduced polynomial coefficients.

Author(s)

Tony Richardson, arichard@stark.cc.oh.us,
adapted by John W. Eaton.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com

Examples

p <- polyreduce(c(0, 0, 1, 2, 3))

Stabilize polynomial

Description

Stabilize the polynomial transfer function by replacing all roots outside the unit circle with their reflection inside the unit circle.

Usage

polystab(a)

Arguments

Value

Vector of stabilized polynomial coefficients.

Author(s)

Paul Kienzle, pkienzle@users.sf.net.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com

Examples

unstable <- c(-0.5, 1)
zplane(unstable, 1)
stable <- polystab(unstable)
zplane(stable, 1)

Power - decibel conversion

Description

Convert power to decibel and decibel to power.

Usage

pow2db(x)

db2pow(x)

Arguments

Value

Converted data, same type and dimensions as x.

Author(s)

P. Sudeepam
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

Examples

db <- pow2db(c(0, 10, 100))
pow <- db2pow(c(-10, 0, 10))

Primitive

Description

Calculate the indefinitive integral of a function.

Usage

primitive(FUN, t, C = 0)

Arguments

Details

This function is a fancy way of calculating the cumulative sum.

Value

Vector of integrated function values.

Author(s)

Juan Pablo Carbajal.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com

See Also

cumsum

Examples

f <- function(t) sin(2 * pi * 3 * t)
t <- c(0, sort(runif(100)))
F <- primitive (f, t, 0)
t_true <- seq(0, 1, length.out = 1e3)
F_true <- (1 - cos(2 * pi * 3 * t_true)) / (2 * pi * 3)
plot (t, F, xlab = "", ylab = "")
lines (t_true, F_true, col = "red")
legend("topright", legend = c("Numerical primitive", "True primitive"),
  lty = c(0, 1), pch = c(1, NA), col = 1:2)

Pulse train

Description

Generate a train of pulses based on samples of a continuous function.

Usage

pulstran(
  t,
  d,
  func,
  fs = 1,
  method = c("linear", "nearest", "cubic", "spline"),
  ...
)

Arguments

Details

Generate the signal y <- sum(func(t + d, ...)) for each d. If d is a matrix of two columns, the first column is the delay d and the second column is the amplitude a, and y <- sum(a * func(t + d)) for each d, a. Clearly, func must be a function which accepts a vector of times. Any extra arguments needed for the function must be tagged on the end.

If instead of a function name you supply a pulse shape sampled at frequency fs (default 1 Hz), an interpolated version of the pulse is added at each delay d. The interpolation stays within the the time range of the delayed pulse. The interpolation method defaults to linear, but it can be any interpolation method accepted by the function interp1

Value

Pulse train generated by the function, returned as a vector.

Author(s)

Sylvain Pelissier, sylvain.pelissier@gmail.com.
Conversion to R by Geert van Boxtel G.J.M.vanBoxtel@gmail.com.

Examples

## periodic rectangular pulse
t <- seq(0, 60, 1/1e3)
d <- cbind(seq(0, 60, 2), sin(2 * pi * 0.05 * seq(0, 60, 2)))
y <- pulstran(t, d, 'rectpuls')
plot(t, y, type = "l", xlab = "Time (s)", ylab = "Waveform",
     main = "Periodic rectangular pulse")

## assymetric sawtooth waveform
fs <- 1e3
t <- seq(0, 1, 1/fs)
d <- seq(0, 1, 1/3)
x <- tripuls(t, 0.2, -1)
y <- pulstran(t, d, x, fs)
plot(t, y, type = "l", xlab = "Time (s)", ylab = "Waveform",
     main = "Asymmetric sawtooth waveform")

## Periodic Gaussian waveform
fs <- 1e7
tc <- 0.00025
t <- seq(-tc, tc, 1/fs)
x <- gauspuls(t, 10e3, 0.5)
plot(t, x, type="l", xlab = "Time (s)", ylab = "Waveform",
     main = "Gaussian pulse")
ts <- seq(0, 0.025, 1/50e3)
d <- cbind(seq(0, 0.025, 1/1e3), sin(2 * pi * 0.1 * (0:25)))
y <- pulstran(ts, d, x, fs)
plot(ts, y, type = "l", xlab = "Time (s)", ylab = "Waveform",
     main = "Gaussian pulse train")

# Custom pulse trains
fnx <- function(x, fn) sin(2 * pi * fn * x) * exp(-fn * abs(x))
ffs <- 1000
tp <- seq(0, 1, 1/ffs)
pp <- fnx(tp, 30)
plot(tp, pp, type = "l",xlab = 'Time (s)', ylab = 'Waveform',
     main = "Custom pulse")
fs <- 2e3
t <- seq(0, 1.2, 1/fs)
d <- seq(0, 1, 1/3)
dd <- cbind(d, 4^-d)
z <- pulstran(t, dd, pp, ffs)
plot(t, z, type = "l", xlab = "Time (s)", ylab = "Waveform",
     main = "Custom pulse train")

Welch’s power spectral density estimate

Description

Compute power spectral density (PSD) using Welch's method.

Usage

pwelch(
  x,
  window = nextpow2(sqrt(NROW(x))),
  overlap = 0.5,
  nfft = if (isScalar(window)) window else length(window),
  fs = 1,
  detrend = c("long-mean", "short-mean", "long-linear", "short-linear", "none"),
  range = if (is.numeric(x)) "half" else "whole"
)

## S3 method for class 'pwelch'
plot(
  x,
  xlab = NULL,
  ylab = NULL,
  main = NULL,
  plot.type = c("spectrum", "cross-spectrum", "phase", "coherence", "transfer"),
  yscale = c("linear", "log", "dB"),
  ...
)

## S3 method for class 'pwelch'
print(
  x,
  plot.type = c("spectrum", "cross-spectrum", "phase", "coherence", "transfer"),
  yscale = c("linear", "log", "dB"),
  xlab = NULL,
  ylab = NULL,
  main = NULL,
  ...
)

Arguments

Details

The Welch method [1] reduces the variance of the periodogram estimate to the PSD by splitting the signal into (usually) overlapping segments and windowing each segment, for instance by a Hamming window. The periodogram is then computed for each segment, and the squared magnitude is computed, which is then averaged for all segments. See also [2].

The spectral density is the mean of the modified periodograms, scaled so that area under the spectrum is the same as the mean square of the data. This equivalence is supposed to be exact, but in practice there is a mismatch of up to 0.5 the data.

In case of multivariate signals, Cross-spectral density, phase, and coherence are also returned. The input data can be demeaned or detrended, overall or for each segment separately.

Value

An object of class "pwelch", which is 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.

spec

Vector (for univariate series) or matrix (for multivariate series) of estimates of the spectral density at frequencies corresponding to freq.

cross

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

phase

NULL for univariate series. For multivariate series, a matrix containing the cross-spectrum phase between different series. The format is the same as cross.

coh

NULL for univariate series. For multivariate series, a matrix containing the squared coherence between different series. The format is the same as cross.

trans

NULL for univariate series. For multivariate series, a matrix containing estimates of the transfer function between different series. The format is the same as cross.

x_len

The length of the input series.

seg_len

The length of each segment making up the averages.

psd_len

The number of frequencies. See freq

nseries

The number of series

series

The name of the series

snames

For multivariate input, the names of the individual series

window

The window used to compute the modified periodogram

fs

The sampling frequency

detrend

Character string specifying detrending option

Note

Unlike the 'Octave' function 'pwelch', the current implementation does not compute confidence intervals because they can be inaccurate in case of overlapping segments.

Author(s)

Peter V. Lanspeary pvl@mecheng.adelaide.edu.au.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

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.

[2] https://en.wikipedia.org/wiki/Welch%27s_method

Examples

fs <- 256
secs <- 10
freq <- 30
ampl <- 1
t <- seq(0, secs, length.out = fs * secs)

x <- ampl * cos(freq * 2 * pi * t) + runif(length(t))
Pxx <- pwelch(x, fs = fs)              # no plot
pwelch(x, fs = fs)                     # plot

# 90 degrees phase shift with with respect to x
y <- ampl * sin(freq * 2 * pi * t) + runif(length(t))
Pxy <- pwelch(cbind(x, y), fs = fs)
plot(Pxy, yscale = "dB")
plot(Pxy, plot.type = "phase")
# note the phase shift around 30 Hz is pi/2
plot(Pxy, plot.type = "coherence")

# Transfer function estimate example
fs <- 1000                 # Sampling frequency
t <- (0:fs) / fs           # One second worth of samples
A <- c(1, 2)               # Sinusoid amplitudes
f <- c(150, 140)           # Sinusoid frequencies
xn <- A[1] * sin(2 * pi * f[1] * t) +
      A[2] * sin(2 * pi * f[2] * t) +  0.1 * runif(length(t))
h <- Ma(rep(1L, 10) / 10)      # Moving average filter
yn <- filter(h, xn)
atfm <- freqz(h, fs = fs)
etfm <- pwelch(cbind(xn, yn), fs = fs)
op <- par(mfrow = c(2, 1))
xl <- "Frequency (Hz)"; yl <- "Magnitude"
plot(atfm$w, abs(atfm$h), type = "l", main = "Actual", xlab = xl, ylab = yl)
plot(etfm$freq, abs(etfm$trans), type = "l", main = "Estimated",
     xlab = xl, ylab = yl)
par(op)

Autoregressive PSD estimate - Yule-Walker method

Description

Calculate Yule-Walker autoregressive power spectral density.

Usage

pyulear(
  x,
  p,
  freq = 256,
  fs = 1,
  range = NULL,
  method = if (length(freq) == 1 && bitwAnd(freq, freq - 1) == 0) "fft" else "poly"
)

Arguments

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.

Note

This function is a wrapper for arburg and ar_psd.

Author(s)

Peter V. Lanspeary, pvl@mecheng.adelaide.edu.au.
Conversion to R by Geert van Boxtel, gjmvanboxtel@gmail.com

See Also

ar_psd, arburg

Examples

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

Kaiser FIR filter design

Description

Compute FIR filter for use with a quasi-perfect reconstruction polyphase-network filter bank.

Usage

qp_kaiser(nb, at, linear = FALSE)

Arguments

Value

The FIR filter coefficients, of class Ma.

Author(s)

André Carezia, andre@carezia.eng.br.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

See Also

Ma, filter, fftfilt, fir2

Examples

freqz(qp_kaiser(1, 20))
freqz(qp_kaiser(1, 40))

Real cepstrum

Description

Return the real cepstrum and minimum-phase reconstruction of a signal

Usage

rceps(x, minphase = FALSE)

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 rceps implements the real cepstrum by computing the inverse of the log-transformed FFT while discarding phase, i.e.,

rceps(x) <- ifft(log(Mag(fft(x))))

The real cepstrum is related to the power spectrum by the relation pceps = 4 * rceps^2.

The function rceps() can also return a minimum-phase reconstruction of the original signal. The concept of minimum phase originates from filtering theory, and denotes a filter transfer function with all of its poles and zeroes in the Z-transform domain lie inside the unit circle on the complex plane. Such a transfer function represents a stable filter.

A minimum-phase signal is a signal that has its energy concentrated near the front of the signal (near time 0). Such signals have many applications, e.g. in seismology and speech analysis.

Value

If minphase equals FALSE, the real cepstrum is returned as a vector. If minphase equals TRUE, a list is returned containing two vectors; y containing the real cepstrum, and ym containing the minimum-phase reconstructed signal

Author(s)

Paul Kienzle, pkienzle@users.sf.net,
Mike Miller.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

References

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

See Also

cceps

Examples

## Simulate a speech signal with a 70 Hz glottal wave
f0 <- 70; fs = 10000           # 100 Hz fundamental, 10 kHz sampling rate
a <- Re(poly(0.985 * exp(1i * pi * c(0.1, -0.1, 0.3, -0.3))))
s <- 0.05 * runif(1024)
s[floor(seq(1, length(s), fs / f0))] <- 1
x <- filter(1, a, s)

## compute real cepstrum and min-phase of x
cep <- rceps(x, TRUE)
hx <- freqz(x, fs = fs)
hxm <- freqz (cep$ym, fs = fs)
len <- 1000 * trunc(min(length(x), length(cep$ym)) / 1000)
time <- 0:(len-1) * 1000 / fs

op <- par(mfcol = c(2, 2))
plot(time, x[1:len], type = "l", ylim = c(-10, 10),
  xlab = "Time (ms)", ylab = "Amplitude",
  main = "Original and reconstructed signals")
lines(time, cep$ym[1:len], col = "red")
legend("topright", legend = c("original", "reconstructed"),
  lty = 1, col = c(1, 2))

plot(time, cep$y[1:len], type = "l",
  xlab = "Quefrency (ms)", ylab = "Amplitude",
  main = "Real cepstrum")

plot (hx$w, log(abs(hx$h)), type = "l",
  xlab = "Frequency (Hz)", ylab = "Magnitude",
  main = "Magnitudes are identical")
lines(hxm$w, log(abs(hxm$h)), col = "red")
legend("topright", legend = c("original", "reconstructed"),
  lty = 1, col = c(1, 2))

phx <- unwrap(Arg(hx$h))
phym <- unwrap(Arg(hxm$h))
range <- c(round(min(phx, phym)), round(max(phx, phym)))
plot (hx$w, phx, type = "l", ylim = range,
  xlab = "Frequency (Hz)", ylab = "Phase",
  main = "Unwrapped phase")
lines(hxm$w, phym, col = "red")
legend("bottomright", legend = c("original", "reconstructed"),
  lty = 1, col = c(1, 2))
par(op)

## confirm the magnitude spectrum is identical in the signal
## and the reconstruction and that there are peaks in the
## cepstrum at 14 ms intervals corresponding to an F0 of 70 Hz.

Rectangular pulse

Description

Return samples of the unit-amplitude rectangular pulse at the times indicated by t.

Usage

rectpuls(t, w = 1)

Arguments

Details

y <- rectpuls(t) returns a continuous, aperiodic, unit-height rectangular pulse at the sample times indicated in array t, centered about t = 0.

y <- rectpuls(t, w) generates a rectangular pulse over the interval from -w/2 to w/2, sampled at times t. This is useful with the function pulstran for generating a series of pulses.

Value

Rectangular pulse of unit amplitude, returned as a vector.

Author(s)

Paul Kienzle, Mike Miller.
Conversion to R by Geert van Boxtel, G.J.M.vanBoxtel@gmail.com.

See Also

pulstran