invfreqz
Determination of discrete filter parameters based on frequency response data.
| Library |
|
Syntax
Function call
-
b,a = invfreqz(___,"trace")— displays a text report on the progress of the iteration.
Arguments
Input arguments
# n,m — desired order
+
scalars
Details
The desired orders of the polynomials in the numerator and denominator, given as positive integers.
# iter — the number of iterations in the search algorithm
+
scalar
Details
The number of iterations in the search algorithm, given as a positive real scalar. Argument iter specifies functions invfreqz finish the iteration when the solution comes together, or after iter iterations, depending on what happens first.
# tol — tolerance
+
0.01 (by default) | scalar
Details
The tolerance is set as a scalar. Function invfreqz defines convergence as the occurrence of the moment when the norm of the (modified) gradient vector is less than the value of the argument tol.
To get a vector of weighting coefficients consisting of one unit, use
invfreqz(h, w, n, m, [], iter, tol)Output arguments
# b,a are the coefficients of the transfer function
+
vectors
Details
The coefficients of the transfer function returned as vectors. The transfer function is expressed in terms of b and a as follows:
| Типы данных |
|
| Support for complex numbers |
Yes |
Examples
Stable approximate transfer function
Details
We convert a simple transfer function into frequency response data, and then back to the original filter coefficients.
import EngeeDSP.Functions: freqz, invfreqz
a = [1 2 3 2 1 4]
b = [1 2 3 2 3]
h, w = freqz(b, a, 64)
bb, aa = invfreqz(h, w, 4, 5)
println("bb = ", bb, "\naa = ", aa)bb = [1.000000000000525 1.9999999999978302 2.9999999999980442 1.999999999998664 2.9999999999957274]
aa = [1.0 1.9999999999999936 2.999999999994893 1.9999999999994074 1.0000000000005878 3.9999999999943077]Let’s draw the zeros and poles of the function, and add a unit circle for clarity.
import Pkg
Pkg.add("ControlSystems")
using ControlSystems
H = tf(vec(bb), vec(aa))
z, p = tzeros(H), poles(H)
pzmap(H)
annotate!([(real(zi)+0.1, imag(zi), text("Zero", 8)) for zi in z])
annotate!([(real(pi)+0.1, imag(pi), text("Pole", 8)) for pi in p])
θ = range(0, 2π, length=200)
plot!(cos.(θ), sin.(θ),
ls=:dash, color=:black,
label="Unit Circle", aspect_ratio=:equal,
xlabel="Real Part", ylabel="Imaginary Part")
display(plot!())
Variables bb and aa are equivalent to b and a accordingly. However, the system is unstable because it has poles outside the unit circle. We use the iterative algorithm of the function invfreqz to find a stable approximation to the system.
bbb, aaa = invfreqz(h, w, 4, 5, [], 30)
println("bbb = ", bbb, "\naaa = ", aaa)bbb = [0.24268035528110668 0.27880314196742667 0.006865008205357481 0.09712543692553093 0.19804400564173713]
aaa = [1.0 -0.8943888916069611 0.6953713986167904 0.9997072238456083 -0.8932896310322123 0.6949656991019518]Let’s make sure that the poles are inside the unit circle.
H = tf(vec(bbb), vec(aaa))
z, p = tzeros(H), poles(H)
pzmap(H)
θ = range(0, 2π, length=200)
plot!(cos.(θ), sin.(θ),
ls=:dash, color=:black,
label="Unit Circle", aspect_ratio=:equal,
xlabel="Real Part", ylabel="Imaginary Part")
display(plot!())
Algorithms
The default function is invfreqz uses the error equation method to determine the best model based on the data. This method finds b and a in
by creating a system of linear equations and solving them using the left division operator \. Here and — Fourier transforms of polynomials a and b accordingly, on the frequency , and — number of frequency points (length of vectors h and w). This algorithm is based on the work of Levi [1].
A more efficient (based on the "exit error") algorithm uses the Gauss method. — Newton with attenuation for iterative search [2], while the output of the first algorithm is used as an initial estimate. This solves the direct problem of minimizing the weighted sum of squared errors between the actual and desired frequency response points.