invfreq function

Inverse Frequency Response

Identify filter parameters from frequency response data.

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

• h: Frequency response, specified as a vector

• w: Angular frequencies at which h is computed, specified as a vector

• nb, na: Desired order of the numerator and denominator polynomials, specified as positive integers.

• wt: Weighting factors, specified as a vector of the same length as w. Default: rep(1, length(w))

• plane: "z" (default) for discrete-time spectra; "s" for continuous-time spectra

• method: minimization method used to solve the normal equations, one of:

  • "ols": ordinary least squares (default)
  • "tls": total least squares
  • "qr": QR decomposition

• norm: logical indicating whether frequencies must be normalized to avoid matrices with rank deficiency. Default: TRUE

Returns

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

• b: moving average (MA) polynomial coefficients
• a: autoregressive (AR) polynomial coefficients

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.

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

References

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

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 .