CC BY
187
УДК 621.396
АНАЛИЗ ИМПУЛЬСНОГО ОТКЛИКА СИСТЕМЫ НА ОСНОВЕ РАЗЛОЖЕНИЯ ПО ЭМПИРИЧЕСКИМ МОДАМ
Ю.А. КРАСНИТСКИЙ
Предложен метод фильтрации неминимально-фазовых импульсных сигналов с применением разложения по эмпирическим модам (РЭМ), применение которого позволяет представить логарифмический комплексный спектр минимально-фазовой части сигнала в виде суммы небольшого числа комплексных характеристических функций (ХФ), что позволяет найти компоненты эмпирической сверточной модели исследуемого сигнала. Разработан ряд программ на МаЙаЬ, которые могут быть использованы при анализе одиночных импульсных сигналов различной физической природы.
Ключевые слова: нестационарный сигнал, минимально-фазовая аппроксимация, логарифмический спектр, преобразование Гильберта, комплексный кепстр.
1. Введение
При обработке нестационарных сигналов широко применяется их разложение в различных базисах. В дополнение к хорошо известным математическим методам в последние годы предложены новые, среди которых особое место занимает разложение по эмпирическим модам (РЭМ, БМБ) [1]. С его помощью действительно-значный сигнал может быть представлен в виде суммы относительно небольшого числа характеристических функций (ХФ), которые не принадлежат к какому-то заранее навязываемому базису: индивидуальный базис извлекается непосредственно из исследуемого сигнала. Однако при работе с комплексно-значными данными возникают некоторые проблемы, поскольку их действительная и мнимая части должны одновременно удовлетворять определенным условиям. Ярким примером могут служить комплексные спектры, обработка которых требует особого внимания.
Интерес к обработке нестационарных импульсных сигналов (или видеоимпульсов) наблюдается в таких областях, как сверхширокополосная радиолокация, разведочная геофизика, гидролокация, исследование природных источников электромагнитного излучения и др. Для их описания чаще всего вводится линейная модель, которая предполагает, что сигнал образован линейной комбинацией
¿0'^) = ¿шрь(*®)4р(*®) (1)
нескольких импульсов, пришедших в точку наблюдения по несовпадающим путям. Такой сигнал, как правило, принадлежит к классу неминимально-фазовых сигналов, а его частотный спектр может быть представлен произведением минимально-фазовой 8трь('а) и неминимально-фазовой 8ар('а) компонент. Многие алгоритмы цифровой обработки используют предположение, что основная информация о сигнале содержится в первом сомножителе спектра (1).
2. Нуль-полюсная модель сигнала и его минимально-фазовый аналог
Любой действительный сигнал можно рассматривать как импульсный отклик (реакцию) некоторой системы, возбуждаемой ¿-импульсом Дирака. 2-преобразование
м м
X Ъ(к ) 7 - П (1 - 2к2 _1)
Н (7) = -*=%------- = Ъ(0) --------- (2)
1 + X а(к ) 7 - к П (1 - Рк2 _1)
к=1 к=1
позволяет описать соответствующую системную функцию в виде отношения двух полиномов по г'1, где корни полинома числителя суть нули функции Н(2), а корни рк полинома знаменателя называются полюсами этой функции.
Существует много алгоритмов расчета коэффициентов а(к) и b(k) в (2) с применением методов моделирования нулей и полюсов (zero-pole (Z-P) modeling) сигнала. Среди других следует отметить фильтры Стиглица-Мак Брайда и Прони [2]. Первый минимизирует среднеквадратическую ошибку
E(z) = X(z) -H(z), (3)
где X(z) - это Z-преобразование исследуемого сигнала x(t). В методе Прони ищутся коэффициенты a(k) так называемой полюсной модели
H (z) =----^--------- (4)
1 + £ a(k ) z ~к
к=1
путем минимизации энергетической ошибки по норме L2 во временной области. Эта ошибка равна
2
ep = £ x(l) - Xp(l)| , (5)
l=0
где xp - последовательность, восстановленная по модели (4), L - число отсчетов сигнала x(t).
Оба фильтра, рассмотренные выше, легко реализуются на базе существующих версий Matlab. Выбор фильтра основывается на предполагаемой структуре модели, поскольку ее внутренние блоки могут быть соединены последовательно, параллельно или более сложным образом. Для сигнала, сформированного линейной комбинацией нескольких задержанных импульсов, пришедших в точку наблюдения по несовпадающим путям, метод Прони представляется более адекватным.
Перед Z-P моделированием следует задать порядки M числителя и N знаменателя системной функции (2). Это достаточно сложная проблема, не имеющая стандартных путей решения. Один из способов состоит в вычислении ошибки (5) при постепенном увеличении величины N с сохранением неизменной разности между N и M. Этот процесс следует прервать при достижении заранее установленного уровня относительной ошибки. Для сохранения устойчивости модели рекомендуется выбирать N < (0.3 - 0.4)L. Величину M обычно задают равной N-1 или меньше.
Рассмотренная процедура снижает уровень шума в восстанавливаемом сигнале xp по сравнению с исходным сигналом x благодаря ограничению числа членов в полиномах заданными значениями N и M. Выбором числа отсчетов в xp можно корректировать влияние прямоугольного окна, неизбежно применяемого к исходным данным перед началом цифровой обработки.
Примером одиночных импульсных сигналов могут служить атмосферики -электромагнитные импульсы, генерируемые молниевыми разрядами. Научная проблема, связанная с их изучением, состоит в решении обратной задачи, состоящей, в частности, в оценке геометрии и электрофизических свойств канала молнии, исходя из структурных особенностей излученного им атмосферика. На рис. 1 a показан сигнал, зарегистрированный на расстоянии около 30 км от разряда и выбранный в качестве пробного при последующем изложении. Здесь же приведен результат восстановления этого сигнала с помощью моделирования по методу Прони (4) при высоком (N=74) порядке фильтра.
На рис. 1 б сравниваются энергии ошибок (3) и (5). Видно, что для рассматриваемого сигнала они одинаковы при порядке фильтра Стиглица-Мак Брайда, примерно в два раза меньшем, чем при моделировании по алгоритму Прони. Тем не менее, как уже упоминалось, выбор метода должен учитывать предполагаемую внутреннюю структуру модели.
Замена сигнала его минимально-фазовым аналогом (МФА, MPA) перед процедурой основной обработки применяется довольно часто. Существует несколько способов оценки МФА, например, на основе классического метода [3] вычисления действительного кепстра. В общем случае при Z-P моделировании некоторых сигналов на картах их нулей и полюсов будут присутствовать нули, расположенные вне единичной окружности, что и служит признаком неминимально-фазовой природы исследуемого сигнала.
Рис. 1. Анализируемый атмосферик - а; ошибка нуль-полюсного моделирования в зависимости от порядка модели - б
На рис. 2 показан (слева внизу) тот же сигнал, что и ранее, но при меньшем значении порядка 2-Р модели. На том же рисунке (слева вверху) представлена его карта нулей и полюсов.
Рис. 2. Минимально-фазовое преобразование
Поскольку некоторые нули расположены вне единичной окружности, рассматриваемый сигнал принадлежит к неминимально-фазовым. Его можно преобразовать в минимальнофазовый, зеркально отобразив упомянутые “лишние” нули внутрь окружности. Полученная при этом преобразовании новая карта нулей и полюсов изображена на рис. 2 справа вверху; под ней показан искомый МФА. Следует отметить, что задержка первого пика в МФА оказывается существенно меньшей, чем у исходного сигнала.
Естественно, что любая операция, совершаемая над сигналом, приводит к потере некоторой доли содержащейся в нем информации. Мерой близости сигналов х и у соответственно до и после обработки может служить коэффициент корреляции
Ь
X ( Х1 - тх )( У1 - ту )
Сху ^ (6)
'Х Ь Ь
X(х1- тх)2 X(У1- ту)2 _/=1 1=1
где тх и ту - средние значения этих сигналов.
Обозначив исходный сигнал (рис. 1) как и, его нуль-полюсную модель как V, а МФА как w, и применив затем соотношение (6), реализуемое функцией еоггеоеГ из пакета МайаЬ, найдем: Сиу = согтсоеДи^) = 0.99943, Сип = 0.74623, Ст = 0.74659.
Уменьшение коэффициента корреляции для МФА можно объяснить очевидным сдвигом этого сигнала на временной оси (рис. 2, справа внизу) относительно исходного.
3. Разложение по эмпирическим модам во временной области
Разложение по эмпирическим модам (РЭМ), предложенное Н. Хуангом и др. в 1998 году, может применяться практически к любым сигналам, в том числе, нелинейным, нестационарным и т.д., независимо от их предполагаемой внутренней структуры [1]. Оно основано на адаптивной процедуре, позволяющей представить сигнал суммой небольшого числа характеристических функций (ХФ), не принадлежащих к какому-то заранее навязываемому базису. ХФ извлекаются непосредственно из обрабатываемого сигнала, формируя адаптированный к нему естественный набор базисных функций. Процесс РЭМ, называемый «просеиванием», отображает цифровой сигнал х[/] в виде
N
х[/] = X сп [/] + г[/], (7)
П = 1
где символами сП[/], п = 1,.. .0 обозначены ХФ, а г[/] - остаточный член разложения или тренд.
Характеристические функции должны удовлетворять двум основным предположениям:
а) число экстремумов и переходов через нуль у этой функции должно быть таким же или отличаться от обрабатываемого сигнала не более, чем на единицу, на каждом этапе просеивания;
б) должна сохраняться локальная симметрия относительно нулевого среднего значения.
Алгоритмы просеивания рассмотрены в многих литературных источниках [4-6].
Примеров сравнительного анализа РЭМ для сигналов и их минимально-фазовых аналогов в доступных источниках не обнаружено. Ниже кратко изложены результаты, полученные на материале, описанном в разделе 2.
На рис. 3 показаны ХФ, обозначенные как с1 . с5, для исходного сигнала (левая колонка) и его минимально-фазового аналога.
Степень близости одноименных ХФ может быть оценена по величине соответствующего коэффициента корреляции (6) и иллюстрируется табл. 1.
Таблица 1
С1-1 С2-2 Сз_з С4-4 С5-5
-0.49183 -0,17519 -0.87089 0.9595 0.994
Как видно, корреляция растет по мере того, как ХФ становится все более гладкой функцией времени.
РЭМ по определению предполагает аддитивную модель сигнала. Поэтому, используя аналогию с вейвлет-анализом, можно утверждать, что наиболее медленные ХФ в (7) соответствуют так называемым аппроксимирующим членам в разложении по вейвлетам, которые определяют основные энергетические соотношения в сигнале. С этой точки зрения сигнал и его МФА оказываются взаимозаменяемыми.
Рис. 3. Характеристические функции для исходного сигнала (левая колонка) и его минимально-фазовый аналог
4. Разложение по эмпирическим модам в спектральной области
Аддитивная модель (7) сигнала во временной области не описывает структуру каскадных линейных систем, выходной сигнал которых образован сверткой импульсных характеристик (ИХ) каскадов. Поэтому РЭМ (7) не позволяет оценить сверточные компоненты, т.е. найти ИХ отдельных звеньев системы. Чтобы решить эту задачу, необходимо выполнить гомоморфное преобразование исходного сигнала u(t), отобразив его из временной области в область логарифмических спектров
SL(iw) = Log[F[w(t)]] = Log(S(iw)) = log S(iw) + ij(w). (8)
Здесь F и Log - символы преобразования Фурье и комплексного логарифмирования соответственно; ф(о>) - фаза комплексного спектра. На эту область сверточные компоненты сигнала отображаются аддитивно, поэтому попытка применения РЭМ к (8) представляется оправданной.
Поскольку функция (8) - комплексная, процедура РЭМ должна учитывать внутреннюю связь действительной и мнимой частей (или модуля и фазы) спектра. При нахождении РЭМ произвольной комплексной функции можно попытаться применить алгоритмы из [4, 6]. Для минимально-фазового эквивалента сигнала необходимые условия выполняются автоматически, поскольку модуль и фаза логарифмического спектра связаны преобразованиями Г ильберта [7, 8]
ln I s (iW) L
iW, inis t/WJI=— P I
W
j(w) = -1 p ¥ ln 1 S(iW) dw, ln|S(iW)|=1 p ¥ -^^dW p J W-w p J W-w
(9)
где символр означает, что интеграл вычисляется в смысле главного значения по Коши.
Если и(^) - минимально-фазовая функция, то согласно терминологии, принятой в теории аналитического сигнала, слагаемые в (8) образуют каноническую пару.
Применяя РЭМ (7) отдельно к каждому слагаемому в (8), комплексный спектр можно представить в виде суммы комплексных ХФ, которую следует рассматривать как аддитивное отображение сверточной модели, описывающей исходный сигнал во времени, в комплексную область.
На рис. 4 в левой и центральной колонках показаны действительные и мнимые компоненты комплексных ХФ для того же самого МФА, рассмотренного выше (рис. 3).
Complex EMD of logspectrum Deconvolution
Frequency, kHz Time, msec
Рис. 4. Действительные и мнимые компоненты комплексных ХФ для рассмотренного выше МФА
По аналогии с (7) комплексный РЭМ формирует отображение
K
SL [n] = X ck [n ] + r[nl (10)
k=1
где Ck[n] - любая из K комплексных пар или ХФ, изображенных на рис. 4. Здесь n = 1,...N, где N - это число отсчетов, учитываемых в логспектре. Последнее слагаемое в правой части (10) представляет собой комплексный остаток, полученный при разложении.
Каждую из комплексных ХФ в (10) можно использовать как начальный материал для деконволюции, т.е. обратного гомоморфного преобразования с целью возврата в исходную временную область и восстановления предполагаемых сверточных компонентов
обрабатываемого сигнала. Необходимая последовательность операций имеет вид
Uk(t) = FT_1[exp(ck)], (11)
где FT-1 - символ обратного преобразования Фурье.
Результаты подстановки в (11) комплексных ХФ (10) показаны в правой колонке на рис. 4. Эти функции представляют собой компоненты эмпирической сверточной модели
анализируемого сигнала. Термин “эмпирический" в данном контексте означает, что детальная структура системы, породившей этот сигнал, a priori неизвестна.
Эмпирические сверточные компоненты (11) можно по отдельности рассматривать как ИХ каскадно включенных звеньев, образующих указанную систему. Скорее всего, эмпирическая модель окажется существенно более сложной, чем можно было бы предположить до проведения описанного анализа. Этого следует ожидать, например, при решении задач, связанных с изучением процессов распространения волн различной физической природы. Для выяснения физической значимости и совместимости эмпирических ИХ, полученных с помощью (11), необходимо отдельное исследование, выходящее за рамки настоящей статьи.
С другой стороны, каждую комплексную ХФ в (10) можно рассматривать как исходный материал для нахождения элементарного комплексного псевдокепстра
Ck(t) = FT-1[ck]. (12)
Хорошо известно [2, 3, 9], что комплексный кепстр представляет собой действительную функцию времени (или задержки), как и в случае так называемого действительного кепстра. Но, в отличие от последнего, в нем сохраняются сведения о фазовых соотношениях в исследуемом сигнале. Благодаря этому комплексный кепстр может быть использован для восстановления сигнала путем деконволюции, основанной на обратных гомоморфных преобразованиях. Можно показать, что между частным комплексным псевдокепстром (12) и действительным псевдокепстром, рассмотренным в [5], существует определенная связь.
5. Заключение
В работе рассмотрено применение разложения по эмпирическим модам к минимальнофазовым аналогам некоторых сигналов. Показано, что результат этого разложения во временной области не удовлетворяет сверточной модели формирования сигнала. Предложен метод анализа структуры сигналов, основанный на РЭМ в области комплексных логарифмических спектров, позволяющий представить сигнал в виде разложения в эмпирическом сверточном базисе, найденном по этому же сигналу. В качестве примера рассмотрена соответствующая последовательность обработки и ее результаты применительно к нестационарному импульсному сигналу - атмос-ферику, порожденному электромагнитным излучением молниевого разряда. Разработаны необходимые вычислительные процедуры с применением среды Matlab.
ANALYSIS OF A SYSTEM PULSE RESPONSE BASED ON EMPIRICAL MODE DECOMPOSITION
Krasnitsky J.A.
The aim of this paper is to work out an approach to filter non-minimal phase transients before their complex spectra processing used empirical mode decomposition. The algorithm is based on zero-pole analysis of the signal with elimination of extraneous zeroes. Subsequent decomposition permits to expand the the logarithmic complex spectrum of minimum phase part into a little amount of complex intrinsic mode functions (IMFs). Deconvolution based on that IMFs allows to find the components of an empirical convolutional model of the signal considered. On the other hand, analysis of these IMFs is a way to produce a complex quasicepstrum. The latter points to an inherent structure of the signal reflected not only the amplitudes but phases relations in the signal too. It is worked out the set of Matlab codes to process the transient pulses of diverse physical nature.
Key words: non-stationary signal, minimum phase approximation, logarithmic spectrum, Hilbert transform, pseudocepstrum.
ЛИТЕРАТУРА
1. Huang N.E. and Shen S.S.P., Eds. Hilbert-Huang transform and its applications. - Singapore, World Scientific, 2005.
2. Hayes M.H. Statistical digital signal processing and modeling. - N.Y., Wiley, 1996.
3. Сильвиа M.T., Робинсон E^. Обратная фильтрация геофизических временных рядов при разведке на нефть и газ. - М.: Недра, 1983.
4. Tanaka Т., Mandic D.P. Complex empirical mode decomposition // IEEE Sign. Proc. Lett., V. 14, No 2, pp. 101-104, 2007.
5. Krasnitsky Yu.A. Pseudocepstral analysis of transient pulses based on empirical mode decomposition. - Transport and Telecommunication, V. 10, No 3, pp. 4-9, 2009.
6. Rilling G., Flandrin P., Goncalves P., Lilly J.M. Bivariate empirical mode decomposition // IEEE Sign. Proc. Lett. V. 14, No 12, pp. 936-939, 2007.
7. Hahn S.L. Hilbert transform in signal processing. - Artech House, 1996.
8. Сороко Л.М. Гильберт - оптика. - М.: Наука, 1981.
9. Г оноровский И.С. Радиотехнические цепи и сигналы. - М.: Радио и связь, 1986.
Сведения об авторе
Краснитский Юрий Александрович, 1938 г.р., окончил Ленинградский институт точной механики и оптики (1961), доктор физико-математических наук, хабилитированный доктор инженерных наук, профессор кафедры телекоммуникаций Института транспорта и связи (Рига, Латвия), автор более 150 научных работ, область научных интересов - радиофизика, радиолокация, цифровая обработка сигналов.
