Цифровое моделирование двумерных линейных систем Текст научной статьи по специальности «Математика»

В.В. Сергеев, А.В. Усачев

ЦИФРОВОЕ МОДЕЛИРОВАНИЕ ДВУМЕРНЫХ ЛИНЕЙНЫХ СИСТЕМ

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

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

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

Предлагаемая методика моделирования применяется к двумерным линейным системам с постоянными параметрами (ЛПП-системам).

Постановка задачи

Сигнал д(х,у) на выходе Лпп-системы равен двумерной апериодической свертке входного сигнала £(х,у) и импульсной характеристике системы Ых,у):

00

д(х,у) = // £(Е,п)Ь(х-Е, у-Т1)<ЗЕ<1г1г (1)

— 00

где х,у - пространственные аргументы сигнала.

в частотной области равенству (1) соответствует выражение

С(“х' V = Р(®Х' V ‘ Н(“х' V' где в, Р, н - Фурье-образы функций д, 1, Ь соответственно, со , а> - пространст-

венные частоты.

Пусть двумерная система задана своей частотной характеристикой, обладающей

свойством радиальной симметрии:

И(«х' “у* = Нр(р)' р = /сох + “у '

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

шагом Т:

£(шТ, пТ) при (т,п) Є [О,М-1]х [О,N-1],

£ »

тп 1 0 в противном случае.

Здесь Мг N — размеры информативной части массива входных отсчетов.

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

Г(<ох, шу) = 0 при |шх1 >2 или 1шу| г £ .

При этом, в силу теоремы Котельникова, дискретизация с шагом Т не сопровожда-

т

ется потерей информации. Тогда непрерывная свертка (1) без потери точности может быть заменена на дискретную |1|

М-1 N-1 = Т2ш=0 п=0

ГДе дк1 = 9(кТ' 1Т)' Ьрд = Ь<РТ' 9Т)•

(2)

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

(к,1) е [к,, ка] * [1,, 1а],

где к1, ка, 11, 1а - границы интересующей нас области выходных отсчетов.

Для вычисления свертки (2) известен быстрый метод, основанный на выполнении трех дискретных преобразований Фурье (ДПФ) - двух прямых и обратного [2].

В рамках рассматриваемой задачи метод ДПФ приобретает ряд особенностей.

С одной стороны, описание системы непосредственно в частотной области позволяет сэкономить одно ДПФ. С другой стороны, оказывается неизвестным "размах” дискретизированной импульсной характеристики, т.е. количество существенно отличных от нуля отсчетов импульсной характеристики, взятых с шагом Т, на положительной

полуоси.

Оценка размаха импульсной характеристики

Для оценки размаха Р импульсной характеристики Ь(х,у) предлагается следующий простой способ, использующий одномерное преобразование Фурье вместо двумерного. Этот способ заключается в определении размаха интегральной функции

00

(х) = / Ых,у)<}у

* — 00

одного аргумента, Фурье-образ которой есть функция Н (р), четным образом про-

г**

долженная на область р < 0. Очевидно (рис. 1), что размах функции (х) равен

размаху импульсной характеристики Ых,у).

Поэтому параметр Р можно определить, например, как наименьшее положительное

число, удовлетворяющее неравенству

Рис. 1. К определению размаха импульсной

характеристики; И = РТ

РТ

/ |Ь (х)12йх

о *

1---------------------< є,

/ |І1 (х) I 2с1х

О Ъ

»

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

Определение размеОа ДПФ для вычисления свертки

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

Эту процедуру удобно описать применительно к одномерному случаю.

Необходимо вычислить

М-1

дк - * УЪ-т' к Є [кі, ка] ,

т=и

причем I при тЄ [0,М-1] и її = Опри Іт| > Р

т т

(3)

Как известно, по методу ДПФ вычисляется циклическая свертка, в которой все участвующие последовательности являются периодическими с периодом 1-размером (длиной) ДПФ. Необходимо, чтобы результаты вычислений циклической свертки совпадали с результатами апериодической свертки (3) для выходного "окна" несмотря на возможные эффекты наложения при переходе от последовательностей конечной длины к периодическим.

На рис. 2 схематически показана взаимосвязь отсчетов входной и выходной последовательностей в апериодической свертке.

Как следует из рис. 2а в общем случае нам нужна информация о значениях £

во входном "окне":

го

т Є [пн , ша] , где пн = к! - Р, та = к2 + Р

к«+Р

ка-Р

I£т I • • • I 1 1 ] ’ ’ * 1^а |

| ••• |9к,Г • Ф • 9ка 1 ''*

• •

р+1-

— ша-т1+1

Р+1

I = Р+М-кп

I = Р+ка+1

Рис. 2. Определение минимально-необходимого размера а) общий случай; I = та-пи+1; б) т-,<0, т3>М-1, Р<М; I * Р + шах{М - к<|, ка + 1}; в) пн <0, та > М - 1, Р £ М; Р = шах{М - к! - 1, ка), I = 2Р + 1

ДПФ

Следовательно, при вычислении можно не принимать во внимание значения вх

ной последовательности £ вне указанного интервала и взять за размер цикличе

ш

ской свертки

I = Ша - пн +1 = ка " + 2Р + 1.

для элементов выходного

исключены

Дополнительное сокращение длины ДПФ возможно в случае, если выходное "окно" или размер импульсной характеристики таковы, что входное "окно" покрывает интервал ненулевых значений последовательности £ , т.е.

ІТІ1 < О, ГП.2 > М-1.

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

(5)

ш1 =0, I = Р + шах{М - к-|, к2+1} г та = I - 1.

Если размах импульсной характеристики настолько велик, что Р > М, то для вычисления каждого требуемого отсчета выходной последовательности нужны все ненулевые отсчеты входной, и в этом случае из (3) следует, что для вычисления

лишь

ными

симметричное

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

т.е. принять (рис. 2в)

Р = шах {М - кі - 1, к2} , I = 2Р + 1, ПН = шіп{0, к.,}, ш2 = ПЧ +1-1. (6)

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

При определении размеров циклической свертки и входного "окна" следует

формулам

(4), (5) или (6) значение I нужно затем увеличить до ближайшего значения, до-

пускающего построение быстрого алгоритма преобразования.

Секционирование свертки

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

Этого можно достичь двумя методами, сравнимыми по сложности вычислений [3].

В первом из них - методе перекрытия с суммированием - входная последовательность ± разбивается на смежные блоки длиной и.

Апериодическая свертка каждого из этих блоков с последовательностью Ь длиной

(V ) ^

2Р + 1 даст последовательность д^ на выходе длиной и + 2Р, которые, перекры-ваясь с соседними блоками в 2Р отсчетах, дают на выходе искомые отсчеты последовательности д^ (рис. За).

В методе перекрытия с накоплением разбиваем входную последовательность на перекрывающиеся блоки, а на выходе отсечением ненужных отсчетов получаем смежные блоки, из которых формируется искомая последовательность (рис. Зб) .

Второй метод несколько эффективнее по сложности вычислений и, следовательно, является более предпочтительным.

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

£^-1) (V)

и и

/

(V-!)

£ ^+1)

Рис, 3. Схема расположения блоков секционированной свертки а) для метода перекрытия с суммированием? б) для метода

перекрытия с накоплением

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

Применение преобразования Хартли

В случае, когда сигнал и частотная характеристика системы вещественны, еще примерно двукратное сокращение вычислительной сложности моделирования и объема обрабатываемых данных обеспечивается использованием вещественного ДПФ в форме

Хартли [4].

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

бы одной из сворачиваемых последовательностей вещественный. Эффективность его применения объясняется тем, что одномерное преобразование Хартли может быть вычислено через одномерное ДПФ половинного размера. Преимущество преобразования Хартли перед алгоритмом совмещенного ДПФ [5] при вычислении свертки состоит в том, что в первом случае промежуточные данные (так же как входные и выходные) вещественны, тогда как преобразование Фурье вещественной функции в общем случае комплексно.

Общий алгоритм моделирования

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

сигнала представляется следующим образом.

1. Определение размеров дискретизированной импульсной характеристики системы с помощью одномерного обратного ДПФ (или обратного ДПФ в форме Хартли, если частотная характеристика вещественна).

2. Расчет параметров циклической свертки, исходя из размеров информативной (ненулевой) части входных отсчетов, размеров импульсной характеристики и ограничений на области входных отсчетов по формулам (4), (5), (6). Определение эффективного размера преобразования (Фурье или Хартли) или (для случая, представ ленного на рис. 2а) оптимального размера блока секционированной свертки.

3. Выполнение прямого преобразования (Фурье или Хартли) для двумерного мае сива входных данных (или - для случая секционированной свертки - для каждого из блоков).

4. Умножение полученного дискретного спектра на отсчеты частотной характеристики .

5. Выполнение обратного преобразования (Фурье или Хартли).

6. Получение искомых отсчетов выходного сигнала. Для случая секционированной свертки - предварительная стыковка блоков перекрытия.

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

Литература

1. П р э т т У. Цифровая обработка изображений. Т. 1. М.: Мир,

1982.

2. Рабинер Л., Голд Б. Теория и применение цифровой обработки сигналов. М.: Мир, 19 78.

3. Нуссбаумер Г. Быстрое преобразование Фурье и алгоритмы вычисления сверток. М.: Радио и связь, 1984.

4. С о р о к о Л.М., Стриж Т.А. Спектральные преобразования на ЭВМ. Дубна: ОИЯИ, 1972.

5. Ярославский Л.п. Введение в цифровую обработку изображений. М.: Наука, 1979.