3.2.4. Новый алгоритм

Существует несколько методов проектирования оптимальных одномерных нерекурсивных цифровых фильтров, которые являются итеративными по своей природе и позволяют найти наилучшую аппроксимацию требуемой частотной характеристики на интервале пользуясь при каждой итерации лишь очень небольшим числом точек. Таким свойством обладают второй алгоритм Ремеза [6], алгоритм Хофстеттера [8] и метод Паркса и Макклеллана [9]. К сожалению, ни один из этих методов нельзя непосредственно распространить на двумерный случай. Во всех трех методах в явном или неявном виде выдвигается требование о том, чтобы функции, используемые в аппроксимации, удовлетворяли условию Хаара; в случае проектирования двумерных фильтров это требование невыполнимо (см. п. 3.1.2). Следует также учесть, что алгоритм Хофстеттера и метод Паркса и Макклеллана основаны на интерполяционной формуле Лагранжа, которую нельзя обобщить на двумерный случай общего вида [27—29]. Однако если бы удалось разработать метод решения задачи двумерной аппроксимации, сходный с перечисленными методами, то это позволило бы значительно уменьшить число ограничений в виде равенства при использовании линейного программирования за счет увеличения числа промежуточных итераций, ведущих к нахождению действительно наилучшей аппроксимации требуемой частотной характеристики.

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

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

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

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

Чтобы начать указанную процедуру, необходимо иметь исходное решение по выбору сокращенного множества точек с ограничениями. Размещение исходных точек с ограничениями не является критичным. Метод определения местоположений точек совершенно не отличается от того, который описан в п. 3.2.3 и был использован для нахождения плотного множества точек с ограничениями. Точки с ограничениями располагаются в узлах прямоугольной решетки с шагом причем принимается, что в каждом из двух направлений импульсная характеристика простирается от до и что К — целочисленная константа. В данном случае константа К выбирается равной двум; однако в других случаях, когда аппроксимация определяется путем решения только одной задачи линейного программирования, выбираются большие значения этой константы. На границах полос точки с ограничениями размещаются точно так же, как описано выше и проиллюстрировано на фиг. 3.5 и 3.6. Как и прежде, точки решетки в переходной полосе не используются в качестве точек с ограничениями. Если не считать выбора исходных точек с ограничениями, все итерации производятся идентично; ниже приводится описание типичной итерации.

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

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

2. Вычисление «непрерывной» частотной характеристики в соответствующей части области частот. Практически частотная

характеристика вычисляется по 256X256 точкам решетки с использованием БПФ.

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

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

5. Построение различных графиков с выходными данными. Такие графики позволяют следить за ходом реализации алгоритма от итерации к итерации.

6. Получение данных для очередной задачи линейного программирования и решение этой задачи.

Первый этап алгоритма заключается в сохранении соответствующих данных решения последней задачи линейного программирования. Должны быть подготовлены следующие данные: коэффициенты импульсной характеристики, максимальное абсолютное значение ошибки и переменные, составляющие базис для двойственной задачи Решение, обеспечиваемое комплексом линейного программирования «Система математического программирования/360», содержит как оптимальный вектор программы для задачи линейного программирования, которая была решена (основная или двойственная задача), так и оптимальный вектор программы для задачи, двойственной к задаче линейного программирования, которая была решена. Если задача линейного программирования была сформулирована в основной форме, то компонентами оптимального вектора программы являются коэффициенты импульсной характеристики и максимальное абсолютное значение ошибки. Максимальное абсолютное значение ошибки также является оптимальным значением целевой функции. Базисные переменные для двойственной задачи являются ненулевыми компонентами оптимального вектора программы для задачи, двойственной к основной задаче. Если задача линейного программирования была сформулирована в двойственной форме, ненулевыми компонентами оптимального вектора программы являются базисные переменные для двойственной задачи. Коэффициенты импульсной характеристики и максимальное абсолютное значение ошибки являются отрицательными значениями компонентов оптимального вектора программы для задачи, двойственной к двойственной задаче. Максимальное абсолютное значение ошибки также

является отрицательным оптимальным значением целевой функции.

Коэффициенты импульсной характеристики требуется сохранить для того, чтобы можно было вычислить частотную характеристику. Максимальное абсолютное значение ошибки сохраняется при каждой итерации с той целью, чтобы его можно было сравнить со значением, полученным в результате выполнения предшествующей итерации. Если это значение ошибки уменьшается, выполнение алгоритма прекращается. Теоретически значение ошибки должно увеличиваться при переходе от данной итерации к следующей. Так почти всегда и происходит. Тем не менее из 250 (ориентировочно) выполненных итераций обнаружились три, давшие некоторое уменьшение этого значения. Однако следует указать, что во всех трех случаях алгоритм уже почти сходился, когда наступало указанное уменьшение абсолютного значения ошибки. Поэтому это явление было объяснено наличием накопленной вычислительной погрешности. Базисные переменные для двойственной задачи следует сохранить для того, чтобы соответствующие точки с ограничениями можно было включить в множество точек с ограничениями для следующей итерации. Одно из требований теоретического характера к новому множеству точек с ограничениями заключается в том, что система уравнений

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

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

Вторая операция алгоритма состоит в вычислении «непрерывной» частотной характеристики фильтра. С помощью БПФ вычисляются 256X256 отсчетов частотной характеристики. Решение об использовании размерности 256X256 массива является компромиссным. С одной стороны, массив должен быть достаточно плотным, чтобы можно было достаточно точно определить положение локальных максимумов и минимумов функции ошибки (четвертый этап алгоритма). С другой стороны, массив должен достаточно хорошо соответствовать имеющейся емкости памяти ЦВМ. Искомое двумерное БПФ вычисляется с использованием одномерного БПФ, как описывается в п. 3.1.1. Используется упрощенный вариант алгоритма БПФ, разработанного Бреннером [30]. Алгоритм был переработан для выполнения одномерного -точечного прямого преобразования с удвоенной точностью. До получения окончательного результата все вычисления БПФ производят в арифметическом устройстве, которое обеспечивает удвоенную точность. Окончательный результат, однако, записывается в память в виде массива размерности 256X256 с обеспечением «однократной» точности. Благодаря этому минимизируется влияние вычислительных погрешностей на точность определения локальных минимумов и максимумов функции ошибки (четвертый этап алгоритма).

Если принимается лишь допущение о действительности импульсной и частотной характеристик, то отсчеты частотной характеристики следует равномерно распределить в пределах области — В этом случае коэффициенты импульсной характеристики удовлетворяют условию симметрии Это условие позволяет ограничиться вычислением только (одномерных) преобразований столбцов. Остальные преобразований столбцов можно определить на основе симметрии. 256 отсчетов частотной характеристики в области можно определить, считая их первой половиной -точечного преобразования. Такое -точечное преобразование можно вычислить как преобразование двух -точечных последовательностей, если воспользоваться фундаментальным понятием алгоритма БПФ с прореживанием по времени. Поскольку такие последовательности являются действительными, они могут быть объединены, а их преобразования могут быть вычислены одновременно; после завершения вычислений результаты разделяются методом выделения четной и нечетной частей. После вычисления преобразований столбцов необходимо перейти к вычислению преобразований строк. Каждая строка является (в общем случае) комплексной. Однако, поскольку известно, что преобразование

каждой строки дает действительный результат, пары строк можно объединить и преобразовать совместно. Результаты разделяются после выполнения вычислений. Таким образом, в общей сложности требуется вычислить одномерных -точечных преобразований.

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

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

Если принято лишь одно допущение о действительности импульсной и частотной характеристик, то должны быть

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

для границы полосы пропускания и

для границы полосы задерживания. Эти отсчеты размещены равномерно по границам с угловым шагом

Поскольку гарантируется нечетность всегда приходится вычислять отсчет при углах и 180°.

Если принимается, что импульсная и частотная характеристики являются действительными и, кроме того, удовлетворяют условиям то частотная характеристика вычисляется только в пределах дуги окружности 45°. В этом случае полное число вычисляемых отсчетов частотной характеристики составляет

для границы полосы пропускания и

для границы полосы задерживания. Отсчеты размещены с угловым шагом

Первый отсчет всегда соответствует углу 0°, последний — углу 45°.

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

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

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

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

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

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

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

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

В случае когда частотная и импульсная характеристики являются действительными и кроме того удовлетворяют условиям правила для

нахождения локальных максимумов и минимумов в угловых точках, на границах и в пределах внутренней области частотной характеристики ничем не отличаются от только что описанных правил. Массив 256X256 отсчетов частотной характеристики заполняет область Три точки и считаются угловыми. Границами являются линии кроме того, могут существовать границы полос пропускания и задерживания. Внутренние точки занимают область треугольной формы

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

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

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

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

2. Контурные диаграммы частотной характеристики фильтра или функции ошибки.

3. Графики частотной характеристики вдоль границ полос.

4. Перспективные виды частотной характеристики или функции ошибки (или обеих функций).

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

Шестой этап алгоритма заключается в получении данных для задачи линейного программирования и решении этой задачи. Пример формулирования задачи проектирования фильтра в виде задачи линейного программирования приведен в п. 3.1.3. Чтобы задачу в такой общей формулировке воплотить в конкретной машинной программе, требуется располагать детальными сведениями, которые существенно зависят от особенностей используемого комплекса линейного программирования. Поэтому данная операция здесь подробно не рассматривается. Путь отыскания решения задачи линейного программирования в симплексном алгоритме описан в литературе [19, 20] и не имеет значения в контексте обсуждаемой темы.

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

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

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

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