 |
Алгоритм деконволюции в частотной области
|
|
|
|
Существуют два принципиально различающихся способа устранения искажений, вносимых реальным регистрирующим прибором в наблюдаемый сигнал. Первый состоит в том, что при помощи искусных технических решений совершенствуют конструкцию реального прибора, однако по мере улучшения характеристик неизбежно возрастают сложность и стоимость аппаратуры. К тому же иногда характеристики сигнала (например, близкая по частоте помеха) не позволяют улучшать передающие свойства прибора.
Другим способом устранения искажений является апостериорная обработка сигнала с целью уменьшения влияния реальных характеристик прибора (сведение их к "идеальному" виду). Вопрос о редукции к идеальному прибору, не искажающему полезный входной сигнал был поставлен еще Рэлеем в 1871 г. в связи с некоторыми задачами спектроскопии и состоял в следующем: можно ли после опыта так обработать функцию, наблюдаемую на выходе некоторого реального прибора, чтобы полностью восстановить сигнал, поступивший на его вход. В такой постановке задача относится к классу некорректно поставленных, и в общем случае не может быть решена, однако на практике при соблюдении некоторых требований (ограниченность спектра рассматриваемого сигнала по частоте, отсутствие в рассматриваемом диапазоне частот нулей в передаточной функции прибора и т.д.) она может быть решена с достаточной точностью и иметь практическое значение.
Различные методы подобного восстановления сигнала хорошо развиты в области обработки изображений (см., например, работу [3]), однако в других областях в настоящее время не получили широкого распространения из-за необходимости цифровой обработки сигнала и большого объема требуемых вычислений (в частности, нам неизвестны медицинские приборы, реализующие такую коррекцию). Тем не менее существующий уровень компьютерной техники позволяет реализовать рассматриваемый алгоритм за приемлемое время (так обработка 1024 отсчетов на IBM PC/AT 80386 без сопроцессора занимает около 5 секунд и может быть ускорена за счет ряда специальных приемов). К тому же, как показывает анализ литературы, большинство исследователей уделяют внимание коррекции только амплитудных, но не фазовых искажений сигналов. В дальнейшем будет показана важность коррекции именно фазовых характеристик сигнала.
|
|
|
Как правило, реальные приборы с достаточной степенью точности являются линейными, и сигнал f(t) на их выходе описывается интегральным уравнением:
(1)
где h(t) - импульсная переходная характеристика прибора, s(t) - сигнал, поступающий на вход прибора. Это уравнение является уравнением свертки, и может быть записано как
(2)
где F(w), H(w), S(w) - фурье-образы функций f(t), h(t) и s(t), т.е. соответственно спектр отклика, передаточная функция системы и спектр входного сигнала. Из уравнения (2) следует, что по крайней мере формально, можно записать
(3)
откуда обратным преобразованием Фурье находим:
(4)
Формула (4) справедлива, если f(t) L2(-,), h(t) L1(-,), F(w)/H(w) L2(-,), и сохраняет силу, если H(w) обращается в нуль на множестве нулевой меры.
На практике применение формулы (4) сопряжено с рядом сложностей:
1. В действительности H(w) отлична от нуля лишь в конечном интервале частот. Радиотехнические приборы разрабатываются таким образом, что H(w) 0 вне частотного диапазона сигнала, поэтому естественной становится коррекция в диапазоне частот, ограниченном частотными характеристиками сигнала и зануление всех частот, где спектральная плотность мощности шума превышает сигнал. Если нуль частотной характеристики попадает в полосу частот сигнала, то некоторые способы разрешения данной проблемы приведены, например, в работе [2].
|
|
|
2. Присутствие в сигнале аддитивной помехи. Выходной сигнал прибора f(t) известен всегда приближенно, и содержит в себе случайную помеху v(t):
f(t) = fт(t) + v(t),
тогда согласно выражению (3) спектр восстановленного сигнала:
S(w) = Sт(w) + V(w)/H(w),
где V(w) - спектр помехи. Спектральная плотность случайной ошибки восстановления (дисперсия решения задачи) равна [3]:
(5)
где Rv(w)= |V(w)|2 - спектральная плотность мощности помехи. Реальная помеха содержит компоненту белого шума, и ее спектральная мощность стремится к конечному пределу при w. Поэтому интеграл (5) может не сходиться либо быть недопустимо большим. В этом случае можно рекомендовать ограничение частотного диапазона коррекции с подавлением ненужных частот спектра.
3. Погрешности, связанные с квантованием уровня сигнала, округлением при машинном представлении действительных чисел, и т.д. Оценка ошибок восстановления, обусловленных данными погрешностями рассматривается в литературе (например, [3]) и должна в каждом случае оцениваться отдельно. В наших задачах сколько-нибудь заметного пpоявления этих факторов замечено не было.
4.Граничные эффекты. Наиболее неприятным фактом при коррекции искажений сигналов являются краевые эффекты. На практике сигнал известен не на всей временной оси, а в диапазоне (-Т, Т) (рис. 1,а). Поэтому для формулы фильтрации в пространственной области (1) возникает вопрос о доопределении сигнала за границами рассматриваемого интервала. В простейшем случае считают, что за границами (-Т, Т) сигнал равен своему математическому ожиданию (рис.1,б). Также используется различная экстраполяция (напр. считают, что сигнал вне (-Т, Т) равен своему значению на границах интервала).
Рис 1. Варианты доопределения сигналов за границами выборки
При цифровой фильтрации (особенно при использовании аппарата преобразований Фурье), часто исходный сигнал считают периодически продолженным с периодом, равным длине исходной последовательности. Так как начало и окончание записи сигнала расположены в произвольные (независимые от сигнала) моменты времени, то длина исследуемого интервала не равна целому числу периодов сигнала, что приводит к несовпадению амплитудного значения (и значения производных амплитуды) сигнала в начальной и конечной точках и, как следствие, к появлению в спектре высокочастотных гармоник, а после фильтрации - к возникновению осцилляций на краях последовательности (эффект Гиббса). Некоторые исследователи рекомендуют в этом случае использование линейной аппроксимации между начальной и конечной точками (рис. 1,в, длина дополнения должна быть не менее удвоенной длительности импульсной характеристики фильтра) или четного продолжения исходной последовательности (рис. 1,г). Проведенные численные эксперименты показали, что хотя эти способы уменьшают высокочастотные выбросы спектра, однако сильно искажают его низкочастотную область. Более удачные результаты при использовании дискретного преобразования Фурье (ДПФ) получаются при использовании так называемых "весовых окон". Обычно они применяются при оценке спектра по малой выборке: предварительно сигнал умножают на быстро спадающую функцию, например, на кривую Гаусса:
|
|
|
(6)
где 2T - длительность выборки, | t | < T, а - параметр (наилучшие результаты получаются при 2.5 < а < 3), затем проводят Фурье-преобразование для оценки спектра. Другие весовые окна рассматриваются, например, в работе [4]. Их применимость в нашем случае определяют два фактора: скорость спада боковых лепестков и малые вычислительные погрешностями. Окно Ханна дает лучшие результаты в смысле отсутствия боковых лепестков, но на границах его значения близки к нулю, поэтому после фильтрации часть информации на концах интервала будет потеряна.
Исходный сигнал можно рассматривать как бесконечный, умноженный на прямоугольное весовое окно, что в частотной области эквивалентно свертке с Фурье-образом этого окна, функцией sinc(x) = sin(x)/x. В результате истинный спектр как бы "размывается" по частоте (боковые лепестки имеют уровень -13.3 дБ), что приводит
|
|
|
к неограниченности полученного спектра (даже если в качестве исходного сигнала взят sin(t));
к тому, что гармоника вносит ненулевой вклад во весь спектр сигнала (т.е. их невозможно разделить при фильтрации).
Предварительное умножение на функцию Гаусса приводит к тому, что свертка производится с ее Фурье-образом, а не с sinc(x), таким образом, за счет некоторого уширения самой спектральной полосы, сильно подавляются ее боковые лепестки (до уровня -45..-55 дБ в зависимости от параметра а), и спектр более точно аппроксимирует спектр реального бесконечного сигнала. На рис.2,а - представлен вычисленный спектр синусоидального сигнала с использованием прямоугольного окна (сплошная линия) и окна Гаусса (теоретически, спектр должен был бы иметь вид дельта-функции). После фильтрации, разумеется, надо вернуть сигнал к первоначальному виду, т.е. разделить на функцию Гаусса. Здесь может возникнуть единственная тонкость, связанная с округлением чисел при машинной обработке.
Рис. 2 Спектр мощности выборки синусоидального сигнала (сплошная линия) и спектр того же сигнала с использованием окна Гаусса (пунктирная линия) (а); видно расширение спектральной полосы и подавление боковых лепестков; исходная выборка (сплошная линия), она же после фильтрации (пунктирная), она же после фильтрации с применением окна Гаусса (штрих-пунктирная) (б)
Использование прямоугольного окна при фильтрации вносит заметные искажения на протяжении примерно трех периодов сигнала с каждой стороны, в то время как окно Гаусса позволяет уменьшить влияние краев выборки до 0.5 периода На рис.2,б представлен синусоидальный сигнал (сплошная линия), он же после проведения узкополосной фильтрации с применением окна Гаусса (штрих-пунктирная линия) и без применения окна Гаусса (штриховая линия). Вторым моментом, связанным с граничными эффектами, является то, что импульсная передаточная характеристика корректирующего фильтра должна быть конечной во времени, и ее длина должна быть много меньше 2Т. На практике это означает невозможность реализации прямоугольной АЧХ фильтра. Поскольку именно такую характеристику считают "идеальной", можно рекомендовать следующий алгоритм построения фильтра. В частотной области определяются амплитудно-фазочастотные характеристики требуемого фильтра, затем при помощи обратного фурье-преобразования вычисляется его импульсная характеристика, которая обрезается до нужной длины (чем короче - тем меньше влияние границ выборки, но и тем ниже порядок получающегося фильтра), и, для повышения гладкости АФЧХ, умножается на функцию Гаусса (6). После этого, вычисляя прямое фурье-преобразование, получим сглаженные АЧХ и ФЧХ нашего фильтра. В работе [5] приведен интересный способ фильтрации, когда помимо восстановления сигнала производится усиление его гармоник пропорционально соотношению сигнал/шум для каждой частоты:
|
|
|
где H(w) и H(w)* - передаточная функция прибора, и ее комплексно сопряженная величина, Rs(w) и Rv(w) - спектральные плотности мощности входного сигнала и шума соответственно.
Этот способ предполагает знание спектральных характеристик шума, и, вообще говоря, искажает спектр входного сигнала.
Вычислительные алгоритмы. Поскольку при обработке сигналов приходится искать компромисс между точностью и скоростью работы математических алгоритмов, то и мы в своем изложении будем уделять основное внимание именно этим двум факторам.
Коррекция исходного сигнала есть не что иное как фильтрация с передаточной функцией 1/H(w). Синтез подобного фильтра аналоговыми методами в общем случае невозможен, то единственным общим способом его реализации является апостериорная цифровая обработка.
Цифровой сигнал представляет собой значения реального аналогового сигнала, взятые через равные промежутки времени. Необходимую частоту дискретизации (время между отсчетами) определяют из теоремы Котельникова (она должна быть больше 2Fмакс, где Fмакс - максимальная присутствующая в спектре сигнала частота). На практике используют частоту, превышающую частоту Котельникова (Найквиста) в 2-5 раз. Фильтрация в цифровой области имеет свои особенности, преимущества и недостатки, и сходна с аналоговой фильтрацией. Простейшим и часто используемым методом цифровой фильтрации является прямое вычисление свертки входного сигнала с импульсной передаточной характеристикой (ИПХ) фильтра по формуле (для нерекурсивного фильтра):
которая представляет собой дискретный аналог формулы (1), здесь Nипх - количество отсчетов импульсной передаточной характеристики. Для вычислений по этой формуле требуется произвести Nc Nипх умножений и столько же сложений (Nc - число отсчетов сигнала), что при длинной ИПХ занимает больше времени, чем реализация такого же фильтра в частотной области.
Согласно формуле (2) свертка сигналов во временной области аналогична перемножению их спектров. Аналогом преобразования Фурье для дискретных сигналов является дискретное преобразование Фурье (ДПФ):
(7)
Его главное отличие от непрерывного преобразования Фурье - цикличность, т.е. и сигнал, и его спектр являются периодическими с периодом N. Наиболее употребительные свойства ДПФ:
A*N-S -- ДПФ комплексно-сопряженного сигнала а;
AS = A*N-S -- ДПФ действительного сигнала;
AS = AN-S -- ДПФ четно-симметричного сигнала;
AS = -AN-S -- ДПФ нечетно-симметричного сигала;
При непосредственных вычислениях по формуле (7) требуется произвести N2 умножений и сложений комплексных чисел, поэтому на практике пользуются быстрым преобразованием Фурье (БПФ), в котором за счет некоторого дополнительного расхода машинной памяти и использования рекурсии число операций уменьшено до N logN, что при больших выборках дает существенный выигрыш во времени. Существуют различные алгоритмы БПФ, рассматриваемые в специальной литературе (например, [3]), суть которых в том, что исходная последовательность делится на несколько меньших частей, для каждой из которых вычисляется ДПФ (время вычисления которого меньше, т.к. N12+N22 < (N1+N2)2), а затем по найденным спектрам вычисляется спектр исходной последовательности.
Кроме алгоритма БПФ, двукратного выигрыша по времени выполнения можно достичь, применяя совмещенные алгоритмы ДПФ. Для последовательностей вещественных чисел коэффициенты Фурье с номерами k и N-k являются комплексно-сопряженными числами (свойство 2). Поэтому можно либо преобразовывать две последовательности одновременно, либо разбить исходную последовательность на две, одновременно преобразовать их, а потом пересчитать результаты на всю последовательность в целом. Используя первый способ, из последовательностей Ak и Bk создают последовательность Ck, такую, что: Ck = Ak + i Bk и вычисляют ее ДПФ (Fr). ДПФ исходных последовательностей находят по соотношениям:
учитывая, что AN-r = Ar, BN-r = Br.
При втором способе последовательность A длины 2N разбивается на две подпоследовательности: четных элементов (Аr_even) и нечетных (Аr_odd). Далее вычисления ведутся как и в первом способе, а затем по формуле
находят ДПФ исходной последовательности. Если в исходной последовательности или в ее ДПФ заведомо много нулей (например, некоторые гармоники в процессе фильтрации должны быть обнулены), то дополнительной экономии времени можно достичь просто их не вычисляя.
Поскольку обычно реальные сигналы вводятся в ЭВМ через АЦП, т.е. являются целыми числами, то естественно для них вести обработку в целочисленной области (которая в 3-10 раз по скорости превышает обработку действительных чисел). Однако в этом случае встает вопрос о точности проводимых преобразований, что ограничивает применимость этих методов для коррекции сигналов. Здесь не приведены численные эксперименты, а показанные ниже методы указаны для полноты обзора, поскольку в литературе встречаются редко.
Наиболее просто использование целых чисел реализуется в методе прямой свертки. Найденные коэффициенты просто домножаются на некоторое число (чтобы получить необходимую точность), а в процессе свертки полученное отфильтрованное значение делится на то же число.
В работах [3, 6] приводится метод квантованного БПФ. Суть его в том, что значения синуса и косинуса квантуются, т.е. принимают дискретные значения, например -1, 0, 1 при трех уровнях и -1, -1/2, 0, 1/2, 1 при пяти уровнях (вообще говоря, это разложение сигнала по другим базисным функциям). В этом случае умножения либо вообще отсутствуют, либо заменяются сдвигами. Для фильтрации изображений этот метод дает вполне приемлемые результаты.
Существует еще ряд базисов (функции Уолша, Хаара, и т.д.), разложение по которым привлекает внимание тем, что операции можно вести с целыми числами, однако мы не исследовали их применимость для цифровой коррекции сигналов.
Обобщение. После проведенных исследований мы остановились на следующей схеме вычислений.
Предварительно измеряем или вычисляем передаточную функцию прибора H(w). Вычисляем передаточная функция восстанавливающего фильтра 1/H(w) (если предполагается прямая свертка во временной области, то затем вычисляем импульсную характеристику (обратное ДПФ от передаточной характеристики), с учетом приведенных выше замечаний.
Исходный сигнал умножаем на функцию окна.
Проводим фильтрацию-восстановление (во временной или в частотной области).
Отфильтрованный сигнал делим на функцию окна.
Вопрос о том, в какой области проводить фильтрацию-восстановление - во временной или в частотной, решается исходя из требуемой точности и скорости обработки. В частотной области можно достичь наивысшей точности при восстановлении и наиболее крутых фронтов фильтров при приемлемом времени вычислений (как уже упоминалось - примерно 5 с на IBM PC/AT для двух последовательностей в 1024 отсчета при одновременной обработке), однако неудобство состоит в необходимости сначала ввести всю последовательность, а затем ее обработать, к тому же при использовании БПФ ограничения накладываются на число точек в последовательности и появляются довольно жесткие требования к объему машинной памяти. Во временной области три последних недостатка отсутствуют, но сильно возрастает время вычислений, которое можно сократить за счет выбора короткой ИПХ (ориентировочно меньше 0.1N, где N - длина исходной последовательности) и проведения операций с целыми числами, однако при этом сильно проигрываем в точности.
|
|
|
