gsignal (version 0.3-1)

invfreq: 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

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

Value

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.

References

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

Examples

# NOT RUN {
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))

# }