Цифровая обработка сигналов

Министерство образования и науки Российской Федерации

Федеральное государственное бюджетное образовательное учреждение

высшего профессионального образования

"Башкирский государственный университет" (БашГУ)

Факультет математики и информационных технологий

Кафедра информационных технологий

Специальность 010501 - Прикладная математика и информатика

Дипломная работа

Цифровая обработка сигналов

Научный руководитель:

д. ф.-м. н., доцент

В.А. Делев

УФА - 2013

Оглавление

Введение

. Линейные системы

.1 Дискретные и непрерывные сигналы

.1.1 Теорема Котельникова

.1.2 Наложение спектров (алиасинг)

.2 Свертка

. Дискретное преобразование Фурье

.1 Преобразования Фурье

.2 Комплексное ДПФ

.3 Двумерное ДПФ

. Восстановление расфокусированных и смазанных изображений

.1 Модель процесса искажения

.1.1 Искажающие функции

.2 Mодель шума

.3 Теорема о свертке

.4 Инверсная фильтрация

.5 Существующие подходы для деконволюции

.6 Сравнение изображений на практике

Заключение

Список литературы

Введение

сигнал преобразование фурье искажение

Цифрова́я обрабо́тка сигна́лов (ЦОС, DSP - англ. <#"730803.files/image001.gif"> будет записан правильно (по теореме Котельникова). Но так как помеха лежит выше частоты 22.05 кГц, то возникнет алиасинг, и помеха "отразится" в нижнюю часть спектра, на частоту около 5 кГц. Если мы теперь попробуем пропустить полученный цифровой сигнал через ЦАП и прослушать результат, то мы услышим на фоне музыки помеху на частоте 5кГц. Помеха "переместилась" из неслышимой ультразвуковой области в слышимую область.

Таким образом, мы видим, что алиасинг - нежелательное явление при дискретизации сигнала. Например, при оцифровке изображения алиасинг может при-вести к дефектам в изображении, таким как "блочные", "пикселизованные" границы или муар.

Как избежать алиасинга? Первый способ - использовать более высокую частоту дискретизации, чтобы весь спектр записываемого сигнала уместился ниже половины частоты дискретизации. Второй способ - искусственно ограничить спектр сигнала перед оцифровкой.

Существуют устройства, называемые фильтрами, которые позволяют изменять спектр сигнала. Например, фильтры низких частот (НЧ-фильтры, low-pass filters) пропускают без изменения все частоты ниже заданной, и удаляют из сиг-нала все частоты выше заданной. Эта граничная частота называется частотой среза (cutoff frequency) фильтра. Одно из важных применений НЧ-фильтров заключается в искусственном ограничении спектра сигнала перед оцифровкой. В этом случае фильтры называются анти-алиасинговыми, т.к. они предотвращают возникновение алиасинга при оцифровке сигнала. Частота среза анти-алиасинговых фильтров устанавливается равной половине частоты дискретизации.

Рассмотрим, что произойдет, если в примере с записью музыки и помехи применить анти-алиасинговый фильтр перед оцифровкой сигнала. Так как частота дискретизации составляет 44.1 кГц, то частота среза фильтра устанавливается на 22 кГц. Таким образом, фильтр будет пропускать без изменения все сигналы, спектр которых лежит ниже 22 кГц (музыку), и подавлять все сигналы, со спектром выше 22 кГц (в том числе - и помеху). После применения фильтра из сигнала исчезнет помеха, и спектр полученного сигнала будет лежать ниже 22 кГц. Когда этот сигнал будет подан на АЦП, алиасинга не возникнет, и по полученной цифровой записи можно будет правильно воссоздать исходную музыку (без помехи).

В реальные АЦП почти всегда встраивается анти-алиасинговый фильтр. Обычно эффект от искусственного ограничения спектра вполне приемлем, в то время как алиасинг - недопустим. Однако не всегда искусственное ограничение спектра так благотворно влияет на записываемый сигнал. Например, при оцифровке музыки на низкой частоте дискретизации 11 кГц приходится отфильтровывать из спектра музыки все частоты выше 5.5 кГц. В результате этого музыка теряет в качестве (хотя обычно такие потери лучше, чем алиасинг). При оцифровке изображений необходимо аккуратно проектировать анти-алиасинговый фильтр, чтобы изменение спектра изображения не повлекло видимых артефактов (таких как пульсации вблизи резких границ).

1.2 Свертка

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

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

 (формула свертки).

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

[n]=x[n]*h[n].

Функция h[n] называется ядром свертки (kernel) или импульсной характеристикой линейной системы.

Обычно все сигналы, обрабатываемые на компьютере, имеют конечную продолжительность (т.е. отличны от нуля лишь на конечном отрезке). Рассмотрим, что происходит с сигналом конечной продолжительности, когда его сворачивают с конечным ядром свертки. Пусть сигнал x[n] отличен от нуля только на отрезке от 0 до N-1 включительно ("имеет длину N"). Пусть ядро свертки h[n] отлично от нуля на отрезке от -m₁ до m₂ включительно, состоящем из M точек (M=m₁+m₂+1). Тогда при подстановке этих сигналов в формулу свертки, мы получим сигнал y[n], который отличен от нуля на отрезке от -m₁ до N-1+m₂ включительно. Таким образом, длина результирующего сигнала равна N+M-1, т.е. сумме длин исходного сигнала и ядра свертки минус один.

Итак, операция свертки расширяет сигнал на M-1 точку, где M - длина ядра свертки.

. Дискретное преобразование Фурье

.1 Преобразования Фурье

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

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

. Непериодический непрерывный сигнал можно разложить в интеграл Фурье.

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

. Непериодический дискретный сигнал можно разложить в интеграл Фурье.

. Периодический дискретный сигнал можно разложить в конечный ряд Фурье.

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

2.2 Комплексное ДПФ

До сих пор мы рассматривали ДПФ от действительных сигналов. Обобщим теперь ДПФ на случай комплексных сигналов. Пусть x[n], n=0,…,N-1 - исходный комплексный сигнал, состоящий из N комплексных чисел. Обозначим X[k], k=0,…N-1 - его комплексный спектр, также состоящий из N комплексных чисел. Тогда справедливы следующие формулы прямого и обратного преобразований Фурье:

,

.

Если по этим формулам разложить в спектр действительный сигнал, то первые N/2+1 комплексных коэффициентов спектра будут совпадать со спектром "обычного" действительного ДПФ, представленным в "комплексном" виде, а остальные коэффициенты будут их симметричным отражением относительно половины частоты дискретизации. Для косинусных коэффициентов отражение четное, а для синусных - нечетное.

.3 Двумерное ДПФ

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

 ,

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

Здесь N₁xN₂ - размер исходного сигнала, он же - размер спектра. k₁ и k₂ - это номера базисных функций (номера коэффициентов двумерного ДПФ, при которых эти функции находятся). Поскольку размер спектра равен размеру исходного сигнала, то k₁ = 0,…,N₁-1; k₂ = 0,…,N₂-1.₁ и n₂ - переменные-аргументы базисных функций. Поскольку область определения базисных функций совпадает с областью определения сигнала, то n₁ = 0,…,N₁-1; n₂ = 0,…,N₂-1.

,

.

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

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

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

Таким образом, эффективный алгоритм вычисления ДПФ изображения заключается в вычислении одномерных БПФ сначала от всех строк, а потом - от всех столбцов изображения.

3. Восстановление расфокусированных и смазанных изображений

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

Правое изображение получено из левого, без использования 48-битного формата (в этом случае будет 100% восстановление исходного изображения) - слева обычный PNG, размытый искусственно. В результате мы видим более качественное изображение.

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

Продемонстрируем это, используя небольшой пример для одномерного случая. Представим, что у нас есть ряд из пикселей со значениями:

x₁ | x₂ | x₃ | x₄… - Исходное изображение.

После искажения значение каждого пикселя суммируется со значением левого, т.е. x’_i = x_i + x_i-1. Надо еще поделить на 2, но опустим этот шаг для простоты. В результате имеем размытое изображение со значениями пикселей:

₁ + x₀ | x₂ + x₁ | x₃ + x₂ | x₄ + x₃… - Размытое изображение.

Теперь попробуем восстановить. Вычтем последовательно по цепочке значения по схеме: из второго пикселя - первый, из третьего - результат второго, из четвертого - результат третьего и так далее. Получаем:

x₁ + x₀ | x₂ - x₀ | x₃ + x₀ | x₄ - x₀… - Восстановленное изображение.

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

3.1 Модель процесса искажения

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

Сформулируем модель процесса искажения следующим образом:

g(x, y) = h(x, y) * f(x, y) + n(x, y) . (1)

Задача восстановления искаженного изображения заключается в нахождении наилучшего приближения f'(x,y) исходного изображения. Рассмотрим каждую составляющую более подробно. С f(x, y) и g(x, y) все достаточно понятно. Разберем функцию h(x, y): в процессе искажения каждый пиксель исходного изображения превращается в пятно для случая расфокусировки и в отрезок для случая простого смазывания. Либо же можно сказать наоборот, что каждый пиксель искаженного изображения "собирается" из пикселей некоторой окрестности исходного изображения. Все это друг на друга накладывается и в результате мы получаем искаженное изображение. То, по какой функции размазывается или собирается один пиксель и называется функцией искажения.

Другие синонимы - PSF (Point spread function, т.е. функция распределения точки), ядро искажающего оператора, kernel и другие. Размерность этой функции меньше размерности самого изображения. К примеру, в начальном рассмотрении примера, размерность функции была 2, т.к. каждый пиксель складывался из двух.

.1.1 Искажающие функции

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

в случае размытия по Гауссу функцией fspecial('gaussian', 30, 8).

в случае смазывания фунцией fspecial('motion', 40, 45).

Операция применения искажающей функции к другой функции (в данном случае к изображению) называется сверткой (convolution), т.е. некоторая область исходного изображения сворачивается в один пиксель искаженного изображения. Обозначается через оператор "*" (не путать с умножением). Математически для изображения f с размерами MxN и искажающей функции h c размерами mxn это записывается так:

g(x,y) = h(x,y)*f(x,y) = h(i,j)f(x+i,y+j). (2)

Где , . Операция, обратная свертке, называется деконволюцией (deconvolution), решение такой задачи нетривиально.

3.2 Mодель шума

Осталось рассмотреть последнее слагаемое, отвечающее за шум, n(x, y) в формуле (1). Причины шума в цифровых сенсорах могут быть разными, но основные - это тепловые колебания и темновые токи <#"730803.files/image014.gif"> исходного изображения:

 . (5)

Это называется инверсной фильтрацией, но на практике практически никогда не работает. Посмотрим на последнее слагаемое в формуле (5) - если функция H(u,v) принимает значения близкие к нулю или нулевые, то вклад этого слагаемого будет доминирующим. Это практически всегда встречается в реальных примерах - для объяснения этого вспомним, как выглядит спектр после преобразования Фурье.

Берем исходное изображение, преобразуем его в полутоновое и используя Matlab, получаем спектр:

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

noise_var = 0.0000001 noise_var = 0.000005

Хорошо видно, что добавление даже очень небольшого шума приводит к значительным помехам, что сильно ограничивает практическое применение метода.

3.5 Существующие подходы для деконволюции

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

 . (6) <#"730803.files/image022.gif"> . (7)

Где y - параметр регуляризации, а P(u, v) - Фурье-преобразование оператора Лапласа (матрицы 3 * 3).

Еще один интересный подход предложили независимо Ричардосн [1972] и Люси [1974]. Метод так и называется "метод Люси-Ричардсона". Его отличительная особенность в том, что он является нелинейным, в отличие от первых трех - что может дать лучший результат. Вторая особенность - метод является итерационным, соответственно возникают трудности с критерием останова итераций. Основная идея состоит в использовании метода максимального правдоподобия, для которого предполагается, что изображение подчиняется распределению Пуассона. Формулы для вычисления без использования преобразования Фурье - все делается в пространственной области:

 . (8)

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

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

3.6 Сравнение изображений на практике

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

Результаты:

Фильтр Винера

Фильтр Люси-Ричардсона

Регуляризация по Тихонову

Слепая деконволюция

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

Запускаем скрипт:

И получаем следующий результат:

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

Теперь рассмотрим пример с реальным смазыванием:

Результат:

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

Заключение

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

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

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

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

Так же в работе вводятся основные понятия ЦОС: линейные системы, свертка, дискретное преобразование Фурье.

Список литературы

1. Гонсалес Р., Вудс Р. Цифровая обработка изображений. - М.: Техносфера, 2005 г. - 1072 с.

. Глинченко А.С. Цифровая обработка сигналов. - Красноярск: Изд-во КГТУ, 2001 г. - 199 с.

. Апальков И.В., Хрящев В.В. Удаление шума из изображений на основе нелинейных алгоритмов с использованием ранговой статистики. - Ярославский государственный университет, 2007. - 105 с.

. Дьяконов В., Абраменкова И. MATLAB. Обработка сигналов и изображений. Специальный справочник. - СПб.: Питер, 2002, 608 с.

5. Сергиенко А.Б. Цифровая обработка сигналов. - СПб.: Питер <http://ru.wikipedia.org/wiki/%D0%9F%D0%B8%D1%82%D0%B5%D1%80_(%D0%B8%D0%B7%D0%B4%D0%B0%D1%82%D0%B5%D0%BB%D1%8C%D1%81%D1%82%D0%B2%D0%BE)>, 2007 г. - 751с.

. Оппенгейм А., Шафер Р. Цифровая обработка сигналов. - М.: Техносфера, 2007 г. - 856 с.

. Гонсалес Р., Вудс Р., Эддинс С. Цифровая обработка изображений в среде MATLAB. - М.: Техносфера, 2006 г. - 616 с.