invfreqz
Определение параметров дискретного фильтра по данным частотной характеристики.
| Библиотека |
|
Синтаксис
Аргументы
Входные аргументы
#
h —
частотная характеристика
вектор
Details
Частотная характеристика, заданная как вектор.
#
n,m —
желаемый порядок
скаляры
Details
Желаемые порядки многочленов в числителе и знаменателе, заданные как положительные целые числа.
#
iter —
количество итераций в алгоритме поиска
скаляр
Details
Количество итераций в алгоритме поиска, заданное как положительный вещественный скаляр. Аргумент iter указывает функции invfreqz завершить итерацию, когда решение сойдется, или после iter итераций, в зависимости от того, что произойдет раньше.
#
tol —
допуск
0.01 (по умолчанию) | скаляр
Details
Допуск, заданный как скаляр. Функция invfreqz определяет сходимость как наступление момента, когда норма (модифицированного) вектора градиента меньше, чем значение аргумента tol.
Для получения вектора весовых коэффициентов, состоящего из одних единиц, используйте
invfreqz(h, w, n, m, [], iter, tol)Выходные аргументы
#
b,a —
коэффициенты передаточной функции
векторы
Details
Коэффициенты передаточной функции, возвращаемые в виде векторов. Передаточная функция выражается через b и a следующим образом:
| Типы данных |
|
| Поддержка комплексных чисел |
Да |
Примеры
Устойчивая приближенная передаточная функция
Details
Преобразуем простую передаточную функцию в данные частотной характеристики, а затем обратно в исходные коэффициенты фильтра.
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]Изобразим нули и полюса функции, добавим единичную окружность для наглядности.
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!())
Переменные bb и aa эквивалентны b и a соответственно. Однако система неустойчива, поскольку имеет полюса вне единичной окружности. Используем итерационный алгоритм функции invfreqz, чтобы найти устойчивое приближение к системе.
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]Убедимся, что полюса находятся внутри единичной окружности.
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!())
Алгоритмы
По умолчанию функция invfreqz использует метод уравнения ошибок для определения наилучшей модели на основе данных. Этот метод находит b и a в
путем создания системы линейных уравнений и их решения с помощью оператора левого деления \. Здесь и — преобразования Фурье полиномов a и b соответственно на частоте , а — количество частотных точек (длина векторов h и w). Этот алгоритм основан на работе Леви [1].
Более эффективный (основанный на «ошибке выхода») алгоритм использует метод Гаусса — Ньютона с затуханием для итеративного поиска [2], при этом в качестве начальной оценки используется выход первого алгоритма. Это решает прямую задачу минимизации взвешенной суммы квадратов ошибок между фактической и желаемой точками частотной характеристики.