4.14. Применение методов оптимизации для расчета БИХ-фильтров

Перейдем к описанию последнего класса методов расчета БИХ-фильтров, называемых методами оптимизации. Отличительная черта этих методов заключается в том, что система уравнений, составленная относительно коэффициентов фильтра, не может быть решена в явной форме.

Поэтому для нахождения коэффициентов приходится использовать расчетные формулы математических методов оптимизации, минимизирущих, согласно выбранному критерию, некоторую ошибку. С помощью последовательных приближений можно в конечном счете свести ошибку к минимуму; можно также задать определенное число выполняемых итераций и после выполнения их считать расчет законченным. В данном разделе будет описано применение нескольких методов оптимизации для расчета фильтров.

1. Минимизация среднеквадратической ошибки

Пусть -преобразование импульсной характеристики БИХ) фильтра имеет вид

т. е. фильтр строится из последовательно соединенных блоков второго порядка. Обозначим заданную амплитудную характеристику фильтра через . Пусть , где , — дискретный ряд необязательно равномерно расположенных частот, на которых вычисляются отклонения получаемой и заданной характеристик фильтра. Тогда, согласно Штейглицу, квадрат суммарной ошибки на всех частотах, рассматриваемый в функции параметров фильтра, можно представить следующим образом:

(4.143)

где -мерный вектор искомых коэффициентов

(4.144)

Минимизация квадрата ошибки (4.143) сводится к нахождению оптимального значения вектора (обозначим его через ), для которого

(4.145)

Задачу минимизации можно решить, применив методы нелинейной оптимизации, скажем алгоритм Флетчера — Пауэла, при использовании которого предполагается, что градиент минимизируемой функции известен. Прежде чем перейти к рассмотрению способа расчета градиента, целесообразно исключить из вычислений коэффициент усиления А фильтра, поскольку он может быть рассчитан аналитически. Введя

(4.146)

и

(4.147)

получим

(4.148)

Оптимальное значение А (равное А*) можно найти, продифференцировав правую часть (4.148) по А и приравняв производную к нулю, что дает

(4.149)

Теперь задача сводится к минимизации функции ошибки

(4.150)

Градиент функции по равен

(4.151)

Второй член в правой части (4.151) равен нулю, так как значение А минимизирует Q. Итак, формулу (4.151) можно записать следующим образом:

Так как

(4.153)

то

(4.154)

Формулу (4.154) можно использовать для вычислений. Итак, все вычисления, необходимые для расчета фильтров с использованием алгоритмов оптимизации типа алгоритма Флетчера—Пауэла, оказываются вполне выполнимыми.

При использовании методов оптимизации учитывается поведение только амплитудной характеристики, поэтому некоторые нолюсы или нули после завершения итераций могут оказаться за пределами единичного круга. В этом случае можно прежде всего заменить полюс с полярными координатами , оказавшийся вне единичного круга, на полюс с координатами находящийся внутри единичного круга.

Амплитудная характеристика фильтра при такой замене остается неизменной, так как полюс заменяется на его зеркальное отображение. Однако после того, как все полюсы оказываются внутри единичного круга, появляется возможность с помощью дополнительного анализа еще больше оптимизировать квадрат ошибки. Такая ситуация возникает достаточно часто, и в этих случаях оптимизация должна производиться двумя этапами:

1. Использование программы оптимизации для минимизации ошибки без каких-либо ограничений на расположение нулей и полюсов.

2. После завершения итераций инвертирование всех полюсов и нулей, оказавшихся вне единичного круга. После этого продолжение программы оптимизации для нахождения нового минимума ошибки.

На фиг. 4.38 в качестве примера приведены кривые для широкополосного дифференциатора, рассчитанного Штейглицем с помощью описанного метода.

Фиг. 4.38. Ошибка аппроксимации амплитудной и фазовой характеристик дифференциатора, рассчитанного с использованием критерия минимума среднеквадратической ошибки (по Штейглицу).

При расчете использовалось

На фиг. 4.38 приведены кривые ошибок аппроксимации амплитудной и фазовой характеристик дифференциатора.

2. Минимизация Lp-ошибки

Дечки показал, что от рассмотренного выше критерия минимума среднеквадратической ошибки можно перейти к критериям ошибки более высокого порядка. Более того, ошибку аппроксимации для характеристики групповой задержки фильтра можно определить так же, как для амплитудной характеристики.

Выразим z-преобразование импульсной характеристики фильтра через z-преобразования К последовательно включенных блоков 2-го порядка, представив полюсы и нули в полярных координатах:

(4.155)

Искомый вектор неизвестных параметров определим следующим образом:

(4.156)

Амплитудная характеристика фильтра будет равна

(4.157)

а характеристика групповой задержки этого же фильтра будет описываться формулой

Задачу расчета БИХ-фильтра по заданной амплитудной характеристике или характеристике групповой задержки можно рассматривать как задачу минимизации ошибок -аппроксимации, определяемых следующими формулами:

(4.159)

(4.160)

Эти формулы представляют ошибки аппроксимации амплитудной характеристики и характеристики групповой задержки соответственно в функции вектора параметров . При и (для всех ) минимизация Lp-ошибки будет идентична минимизации по критерию минимума среднеквадратической ошибки, рассмотренной в предыдущем разделе. Можно показать, что случай будет соответствовать критерию Чебышева (т. е. минимаксному критерию).

Фиг. 4.39. Ошибка аппроксимации амплитудной характеристики дифференциатора, рассчитанного с использованием критерия ошибки 4-го порядка (по Дечки).

Фиг. 4.40. Многополосный фильтр, рассчитанный с использованием методов оптимизации (по Дечки).

Фиг. 4.41. Выравнивание характеристики групповой задержки фильтра нижних частот с использованием выравнивающего фильтра, рассчитанного методами оптимизации (по Дечки).

Итак, задача расчета коэффициентов фильтра с использованием Lp-критерия сводится к задаче минимизации ошибок или путем подбора вектора . Можно показать, что если и весовая функция положительна, то ошибка имеет локальный минимум. Это дает возможность для нахождения вектора параметров , минимизирующего соответствующую ошибку, использовать алгоритмы минимизации без ограничений типа алгоритма Флетчера — Пауэла.

На фиг. 4.39-4.41 приведены примеры использования критерия минимума Lp-ошибки, взятые из работы Дечки. На фиг. 4.39 представлена ошибка аппроксимации одно каскадного широкополосного дифференциатора, при расчете которого было взято . В этом примере минимизировалась ошибка аппроксимации амплитудной характеристики, причем для любой из возможных частот величина ошибки не превышает 1%. На фиг. 4.40 показана амплитудная характеристика рассчитанного этим методом фильтра, который был получен из фильтра 10-го порядка с двумя полосами пропускания и тремя полосами непропускания (К = 5).

В полосе ненропускания величина ошибки составляет приблизительно 0,1 (что обеспечивает ослабление на 20 дБ). Последний пример приведен на фиг. 4.41, где представлены характеристики групповой задержки исходного эллиптического фильтра и эллиптического фильтра, полученного после выравнивания его групповой задержки. Выравнивающая цепь состояла из включенных последовательно с исходным эллиптическим фильтром всепропускающих цепей, не оказывающих влияния на амплитудную характеристику фильтра. Порядок выравнивающей всепропускающей цепи был равен , индекс ошибки . Как видно из фиг. 4.41, после выравнивания пульсации групповой задержки фильтра стали равновеликими.

3. Оптимизация в w-плоскости с использованием всепропускающих цепей

Весьма простая методика оптимизации, предложенная Дечки, может быть использована в случае, когда рассчитываемый БИХ-фильтр имеет равновеликие пульсации в полосе пропускания и обеспечивает аппроксимацию с равновеликими пульсациями произвольной характеристики в полосе непропускания или в полосе пропускания. Рассмотрим эту методику. Для этого запишем квадрат амплитудной характеристики фильтра в виде

(4.161)

где — рациональная передаточная функция, подобная, но не идентичная рациональной функции Чебышева, использовавшейся при расчете эллиптических фильтров. Вместо того чтобы сразу найти функцию , целесообразно сначала перенести решение задачи аппроксимации из z-плоскости в некоторую новую плоскость (назовем ее w-плоскостью), такую, чтобы полоса пропускания (или непропускания) фильтра отображалась на всю мнимую ось в w-плоскости. В этом случае оказывается возможным достаточно просто записывать передаточные функции всепропускающих цепей в -плоскости; используя эти всепропускающие функции, можно получить такую функцию которая будет осциллировать между 0 и 1 при изменении w вдоль мнимой оси от 0 до Таким образом, функция будет иметь в w-плоскости равновеликие пульсации. В зависимости от характера отображения эти равновеликие пульсации можно отобразить либо в полосу пропускания фильтра в z-плоскости, либо в полосу непропускания. При этом поведение характеристики фильтра в другой полосе будет полностью определяться еще не найденными значениями коэффициентов всепропускающего фильтра.

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

Чтобы определить передаточную функцию, модуль которой постоянен вдоль всей мнимой оси, рассмотрим всепропускающую функцию вида

где либо действительные, либо образуют комплексно сопряженные пары. Так как , то можно записать в виде

(4.163)

где

(4.164)

причем

Введем действительную функцию следующим образом:

(4.166)

Аналитически продолжая (4.166), получим

(4.167)

Таким образом, найдена искомая функция имеющая в -плоскости равновеликие пульсации модуля вдоль оси независимо от значений коэффициентов всепропускающего фильтра. Рассмотрим теперь случай отображения из z-плоскости в плоскость, соответствующий равновеликим пульсациям в полосе пропускания при произвольной характеристике в полосе непропускания. Аналогично можно было бы рассмотреть и другой случай отображения, соответствующий равновеликим пульсациям в полосе непропускания при произвольной характеристике в полосе пропускания, но так как этот случай обычно представляет меньший интерес, то ниже он не рассматривается.

Для отображения полосы пропускания фильтра из z-плоскости на всю мнимую ось в w-плоскости воспользуемся следующим преобразованием:

Обратное преобразование из w-плоскости в z-плоскость можно найти, решив (4.168) относительно z:

(4.169)

где

С помощью преобразования (4.168) дуга единичной окружности из z-плоскости, соответствующая частотам (т. е. полосе пропускания фильтра), отображается на всю мнимую ось в w-плоскости (предполагается, что фильтр имеет единственную полосу пропускания и две или более полосы непропускания). Полосы пропускания фильтра из z-плоскости отображаются следующим образом:

(4.171)

Здесь — действительная часть w. Для фильтра нижних частот во всех приведенных выше формулах

Преобразование z-плоскости в w-плоскость иллюстрируется на фиг. 4.42, где в каждой из плоскостей изображены область пропускания (с равновеликими пульсациями и области непропускания.

Теперь остается лишь привести методику расчета коэффициентов входящих в формулу (4.167), которые дали бы возможность аппроксимировать произвольную амплитудную характеристику в полосе непропускания, а также методику получения передаточной функции , или, что то же самое, функции , по функции Рассмотрим сначала вторую из этих двух задач как более простую.

Если нули функции расположены на единичной окружности (как это обычно имеет место при аппроксимации с равновеликими пульсациями), то корни в w-плоскости будут действительными и иметь четную кратность, так как в этом случае оба комплексно сопряженных нуля будут отображаться в одну и ту же точку на действительной оси в w-плоскости.

Фиг. 4.42. Отображение из z-плоскости в w-плоскость (по Дечки).

Поэтому выражение (4.167) можно представить следующим образом:

где — преобразованные нули фильтра. Введя вспомогательные многочлены и , равные

и

(4.174)

запишем формулу (4.161) в виде

из которой после разложения на множители получим

(4.176)

Корни функции располагаются в правой полуплоскости w, что гарантирует устойчивость искомого фильтра с передаточной функцией Н (z). [Кроме того, при расчетах удобнее находить корни функции в -плоскости и преобразовывать их затем обратно в z-плоскость, так как в w-плоскости они обычно легче разделяются, чем корни Н(z) в z-плоскости.]

Наконец приведем алгоритм расчета коэффициентов , обеспечивающих заданную форму амплитудной характеристики в полосе непропускания. Простой способ аппроксимации получается при использовании для описания амплитудной характеристики в полосе непропускания функции ее логарифма:

откуда

При этом для расчета значений , таких, чтобы величина ошибки аппроксимации функции а, задаваемой формулой (4.178), была минимаксной для всей полосы непропускания, можно использовать достаточно простые методы (например, алгоритм Ремеза).

Таким образом, выше было показано, что БИХ-фильтры с равновеликими пульсациями в полосе пропускания (или в полосе непропускания) и произвольной характеристикой в полосах непропускания (или в полосах пропускания) можно рассчитывать путем перенесения решения задачи аппроксимации из z-плоскости в (-плоскость, такую, чтобы полоса пропускания фильтра в z-плоскости отображалась на всю мнимую ось в w-плоскости. В этой новой w-плоскости синтезируется всепропускающая функция, модуль которой постоянен на всей мнимой оси. Затем из этой всепропус-кающей функции с помощью простой подстановки формируется передаточная функция, модуль которой имеет вдоль мнимой оси в w-плоскости равновеликие пульсации независимо от значений коэффициентов всепропускающей функции. Наконец, с помощью простой методики рассчитываются оптимальные значения коэффициентов всепропускающего фильтра, которые используются для аппроксимации в z-плоскости требуемой частотной характеристики в полосе непропускания. Проиллюстрируем применение этого метода на двух примерах, взятых из работы Дечки.

Фиг. 4.43. Амплитудная характеристика режекторного фильтра, рассчитанного в w-плоскости (по Дечки).

Фиг. 4.44. Амплитудная характеристика фильтра нижних частот, рассчитанного в w-плоскости (по Дечки).

На фиг. 4.43 и 4.44 изображены амплитудные характеристики (в логарифмическом масштабе) двух фильтров с равновеликими пульсациями в полосе пропускания и произвольными характеристиками в полосе непропускания. В примере, представленном на фиг. 4.43, заданный уровень пульсаций в полосе пропускания составлял 1 дБ, а характеристика в полосе непропускания должна была удовлетворять следующим условиям:

Порядок рассчитанного фильтра оказался равным восьми, причем этот фильтр удовлетворяет заданным характеристикам с точностью до 0,7 дБ. Второй пример, приведенный на фиг. 4.44, соответствует фильтру нижних частот с линейным (в логарифмическом масштабе) увеличением ослабления в полосе непропускания.

4. Расчет БИХ-фильтров методами линейного программирования

Методы линейного программирования могут быть использованы для расчета БИХ-фильтров, обеспечивающих аппроксимацию с равновеликими пульсациями заданной амплитудной характеристики. Если передаточная функция цифрового фильтра равна

то можно представить в виде

(4.180)

где

(4.1816)

Поэтому квадрат амплитудной характеристики фильтра [т. е. значения (4.180) на единичной окружности] равен отношению тригонометрических полиномов:

(4.182)

Обе функции линейно зависят от коэффициентов Рассмотрим, каким образом можно использовать методы линейного программирования для нахождения таких значений коэффициентов которые обеспечили бы аппроксимацию заданного квадрата амплитудной характеристики функцией причем максимум ошибки аппроксимации был бы минимизирован (т. е. чтобы аппроксимация имела равновеликие пульсации).

Итак, задача аппроксимации заданной функции квадрата амплитудной характеристики сводится к нахождению таких коэффициентов фильтра, при которых

(4.183)

Здесь — функция допуска для ошибки аппроксимации, позволяющая учитывать неодинаковый вес ошибок аппроксимации на различных частотах. Функции обычно известны (или, как будет показано в приведенном ниже примере, зависят от некоторого параметра), поэтому неравенство (4.183) можно представить с помощью следующей системы неравенств, линейных относительно неизвестных

или

(4.185)

Эти, а также следующие дополнительные линейные неравенства:

(4.186)

(4.187)

полностью определяют задачу аппроксимации. Для решения системы линейных неравенств (4.185)-(4.187) из левой части каждого из них вычитается вспомогательная переменная w, которая затем минимизируется. Если значение w оказывается равным нулю, то это означает, что решение задачи аппроксимации существует, причем значения коэффициентов можно получить обычным методом линейного программирования. Если же оказывается, что w О, то это означает, что решения задачи аппроксимации не существует, поэтому для ее решения нужно изменить либо либо .

Рассмотрим в качестве примера расчет фильтра нижних частот с равновеликими пульсациями и коэффициентом передачи, равным 1 в полосе пропускания и 0 в полосе непропускания. Допустим, что максимальная ошибка аппроксимации равна в полосе непропускания и (величина постоянной К выбирается разработчиком) в полосе пропускания. Заданная амплитудная характеристика такого фильтра нижних частот изображена на фиг. 4.45, а. Величина неизвестна, причем в процессе расчета она должна быть минимизирована (результирующий фильтр в данном случае будет, очевидно, эллиптическим, но здесь он используется только для иллюстрации метода). Заданная функция квадрата амплитудной характеристики фильтра, равная квадрату функции, изображенной на фиг. 4.45, а, приведена на фиг. 4.45, б.

Фиг. 4.45. Исходные характеристики фильтра нижних частот, используемые при расчете его методом линейного программирования.

Ее можно рассматривать как аппроксимацию с равновеликими пульсациями функции , изображенной на фиг. 4.45, в; функция амплитуды ошибки аппроксимации представлена на фиг. 4.45, г. [Читатель может убедиться в том, что сумма дает верхнюю границу квадрата амплитудной характеристики, а разность — ее нижнюю границу.] Решив неравенства (4.185)-(4.187) при заданных и , можно найти оценку . Для фильтров нижних частот величина ограничена по определению следующими пределами:

Это позволяет получить начальную оценку и тем самым задать функции на фиг. 4.45. Метод линейного программирования дает возможность для выбранного значения определить, имеет ли заданная система неравенств какое-либо решение.

Фиг. 4.46. Амплитудная характеристика фильтра нижних частот, рассчитанного методом линейного программирования.

Если решения не существует, величину следует увеличивать до до тех пор, пока не будет получено решение. Если же система неравенств имеет решение, начальное значение следует заменить на минимальное значение , для которого решение еще существует. Поднимая методом последовательных приближений нижнюю границу (для которой решения не существует) и опуская верхнюю границу, можно с любой заданной точностью (по крайней мере теоретически) найти минимальное значение .

Хотя при использовании рассматриваемого метода и встречаются трудности, связанные с чувствительностью коэффициентов фильтра к выбору функции квадрата его амплитудной характеристики, тем не менее он часто и с успехом применялся для расчета цифровых фильтров. Так, на фиг. 4.46 представлена амплитудная характеристика (в логарифмическом масштабе) фильтра нижних частот, а на фиг. 4.47 — функция ошибки аппроксимации широкополосного дифференциатора, рассчитанных описанным методом.

Фиг. 4.47. Ошибка аппроксимации амплитудной характеристики дифференциатора, рассчитанного методом линейного программирования.

Фильтр нижних частот имеет шестой порядок граничные частоты полосы пропускания и полосы непропускания составляют 0,20 и, 0,25 соответственно, , окончательное значение максимума ошибки аппроксимации .

Порядок дифференциатора равен четырем, амплитуда пульсаций в диапазоне частот составляет 0,00000763 (поведение характеристики в диапазоне не задавалось). Из фиг. 4.47 видно, что только пульсации относительной ошибки аппроксимации амплитудной характеристики дифференциатора являются равновеликими. Две прямые линии являются граничными для функции абсолютной ошибки. Наибольшее несовпадение с граничными линиями является результатом погрешности алгоритма линейного программирования при таких малых значениях .