invfreqs

Определение параметров непрерывного фильтра по данным частотной характеристики.

Библиотека

EngeeDSP

Синтаксис

Вызов функции

b,a = invfreqs(h,w,n,m) — возвращает векторы вещественных коэффициентов числителя и знаменателя b и a передаточной функции h.

b,a = invfreqs(h,w,n,m,wt) — взвешивает ошибки аппроксимации в зависимости от частоты с помощью аргумента wt.

b,a = invfreqs(___,iter) — предоставляет алгоритм, гарантирующий устойчивость результирующей линейной системы путем поиска наилучшей аппроксимации с использованием численной итерационной схемы. Данный синтаксис может включать любую комбинацию входных аргументов из предыдущих вариантов синтаксиса.

b,a = invfreqs(___,tol) — использует аргумент tol для определения сходимости итерационного алгоритма.

b,a = invfreqs(___,"trace") — выводит текстовый отчет о ходе выполнения итерации.

b,a = invfreqs(h,w,"complex",n,m,___) — создает комплексный фильтр. В этом случае симметрия не соблюдается, а частота задается в радианах в диапазоне от до .

Аргументы

Входные аргументы

# h — частотная характеристика
вектор

Details

Частотная характеристика, заданная как вектор.

# w — угловые частоты
вектор

Details

Угловые частоты, на которых вычисляется аргумент h, заданные как вектор.

# n,m — желаемый порядок
скаляры

Details

Желаемые порядки многочленов в числителе и знаменателе, заданные как положительные целые числа.

Типы данных

Float32, Float64

# wt — весовые коэффициенты
вектор

Details

Весовые коэффициенты, заданные как вектор. Вектор весовых коэффициентов wt имеет ту же длину, что и вектор w.

Типы данных

Float32, Float64

# iter — количество итераций в алгоритме поиска
скаляр

Details

Количество итераций в алгоритме поиска, заданное как положительный вещественный скаляр. Аргумент iter указывает функции invfreqs завершить итерацию, когда алгоритм сойдется к решению, или после iter итераций, в зависимости от того, что произойдет раньше.

# tol — допуск
0.01 (по умолчанию) | скаляр

Details

Допуск, заданный как скаляр. Функция invfreqs определяет сходимость как наступление момента, когда норма (модифицированного) вектора градиента меньше, чем значение аргумента tol.

Для получения вектора весовых коэффициентов, состоящего из одних единиц, используйте

invfreqs(h, w, n, m, [], iter, tol)

Выходные аргументы

# b,a — коэффициенты передаточной функции
векторы

Details

Коэффициенты передаточной функции, возвращаемые в виде векторов. Передаточная функция выражается через b и a следующим образом:

Типы данных	`Float64`, `Float32`
Поддержка комплексных чисел	Да

Примеры

Преобразование передаточной функции в частотную характеристику

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 и aa эквивалентны b и a соответственно. Однако система неустойчива, поскольку aa имеет полюса с положительной вещественной частью. Рассмотрим полюса bb и aa, добавим единичную окружность для наглядности.

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 1

Используем итерационный алгоритм функции 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 и 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!())

invfreqs 2

Советы

При построении моделей более высокого порядка с использованием высоких частот важно масштабировать частоты, разделив их, например, на половину самой высокой частоты, присутствующей в w, чтобы получить хорошо обусловленные значения a и b. Это соответствует изменению масштаба времени.

Алгоритмы

По умолчанию функция invfreqs использует метод уравнения ошибок для определения наилучшей модели на основе данных. Этот метод находит b и a в

путем создания системы линейных уравнений и их решения с помощью оператора левого деления \. Здесь и — преобразования Фурье полиномов a и b соответственно на частоте , а — количество частотных точек (длина векторов h и w). Этот алгоритм основан на работе Леви [1].

Более эффективный (основанный на «ошибке выхода») алгоритм использует метод Гаусса — Ньютона с затуханием для итеративного поиска [2], при этом в качестве начальной оценки используется выход первого алгоритма. Это решает прямую задачу минимизации взвешенной суммы квадратов ошибок между фактической и желаемой точками частотной характеристики.

Литература

Levi, E. C. «Complex-Curve Fitting.» IRE Transactions on Automatic Control. Vol. AC-4, 1959, pp. 37–44.
Dennis, J. E., Jr., and R. B. Schnabel. Numerical Methods for Unconstrained Optimization and Nonlinear Equations. Englewood Cliffs, NJ: Prentice-Hall, 1983.