invfreqs
Determination of continuous filter parameters based on frequency response data.
| Library |
|
Syntax
Function call
-
b,a = invfreqs(___,"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 invfreqs complete the iteration when the algorithm converges to a solution, 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 invfreqs 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
invfreqs(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
Converting a transfer function to a frequency response
Details
We convert a simple transfer function into frequency response data, and then back to the original filter coefficients.
import EngeeDSP.Functions: invfreqs
function freqs(b, a, n::Int)
w = exp10.(range(-1, stop=1, length=n))
H = complex(zeros(n))
for i in 1:n
s = im * w[i]
H[i] = polyval(b, s) / polyval(a, s)
end
return H, w
end
function polyval(coeffs, x)
val = 0.0 + 0.0im
for i in 1:length(coeffs)
val += coeffs[i] * x^(length(coeffs)-i)
end
return val
end
h, w = freqs(b, a, 64)
bb, aa = invfreqs(h, w, 4, 5)
println("bb = ", bb, "\naa = ", aa)bb = [0.9999999999990669 1.9999999990772468 2.999999999727699 2.0000000001897797 2.999999998506294]
aa = [1.0 1.9999999990749304 2.9999999998375264 1.9999999999506228 1.000000001133613 3.9999999974508746]Variables bb and aa are equivalent to b and a accordingly. However, the system is unstable because aa it has poles with a positive real part. Consider the poles bb and aa, we 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!())
We use the iterative algorithm of the function invfreqs to find a stable approximation to the system.
bbb, aaa = invfreqs(h, w, 4, 5, [], 30)
println("bbb = ", bbb, "\naaa = ", aaa)bbb = [0.785752197354587 3.361463179107808 3.3188282441237305 3.763898679714697 -0.1609352003083818]
aaa = [1.0 4.552000944361605 9.45821054631851 8.730229236531377 7.354415003366297 0.0005058837699203762]Let’s make sure that the system is stable by plotting a graph of new poles. bbb and aaa.
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 invfreqs 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.