[医]常见问题
基于频率响应数据确定连续滤波器参数。
库::`工程师`
语法
函数调用
* [参数:ba]=invfreqs(___,"跟踪") -显示迭代进度的文本报告。
争论
例子:
将传递函数转换为频率响应
Details
我们将一个简单的传递函数转换为频率响应数据,然后返回到原始滤波器系数。
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]
变量 bb,bb 和 机管局 等同于 b 和 a 相应地。 但是,系统不稳定,因为 机管局 它有一个积极的实部极点。 考虑两极 bb,bb 和 机管局,为了清楚起见,我们添加了一个单位圆。
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!())
我们使用函数的迭代算法 *invfreqs* 以找到系统的稳定近似值。
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]
让我们通过绘制新极点的图形来确保系统稳定。 bbb,bbb 和 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!())