Содержание к диссертации
Введение
Глава 1. Анализ предметной области. Постановка задач исследования 13
1.1. Общие сведения об анализе временных рядов 13
1.2. Краткие сведения о современных методах анализа НВР
1.2.1. Мгновенный спектр сигнала 15
1.2.2. Вейвлет-анализ 18
1.2.3. Метод «ryceHHH,a»-SSA 23
1.3. Преобразование Хуанга-Гильберта 28
1.3.1. Исходный метод преобразования Хуанга-Гильберта 28
1.3.2. Модифицированный метод ННТ 47
1.4. Постановка задач исследования 56
Глава 2. Разработка методики использования метода ННТ для анализа ВР 58
2.1. Обоснование выбора входных параметров метода ННТ 59
2.1.1. Исследование влияния значения ансамблевого числаNE на результаты декомпозиции методом ННТ 60
2.1.2. Анализ зависимости погрешности декомпозиции ВР от SNE 63
2.2. Проверка рекомендаций по выбору входных параметров на основе анализа результатов применения метода ННТ к детерминированным ВР 68
2.2.1. Исследование ВР, представляющего сумму двух периодических ВР
2.2.2. Исследование особенностей декомпозиции ВР, содержащего мгновенные значения сигнала со скачкообразно изменяющейся частотой 74
2.2.3. Исследование точности нахождения функции AM периодических ВР с AM методом DQ 80
2.2.4. Исследование точности нахождения функции мгновенной частоты ВР, содержащих отсчеты ЧМ дискретных сигналов, методом DQ 90
2.2.5. Выводы по результатам проверки обоснованности рекомендаций по выбору входных параметров ННТ и результатам исследования детерминированных ВР 98
Проверка рекомендаций по выбору входных параметров ННТ
на основе статистического моделирования 99
2.3.1. Исследование точности декомпозиции ВР, представляющего собой смесь белого шума и периодической составляющей, на основе статистического моделирования 102
2.3.2. Исследование точности декомпозиции ВР, представляющего собой смесь белого шума и нелинейного тренда,
на основе статистического моделирования 104
2.3.3. Исследование точности нахождения значений функции AM ВР, представляющего собой смесь отсчетов AM сигнала и белого шума 105
2.3.4. Исследование точности нахождения значений МЧ ВР, представляющего собой смесь отсчетов ЧМ сигнала и белого шума 106
2.3.5. Выводы по результатам проверки обоснованности рекомендаций по выбору входных параметров ННТ на основе статистического моделирования 108
2.4. Методика использования метода ННТ для анализа BF 109
2.5. Выводы 111
Глава 3. Применение методики анализа ВР, основанной на преобразовании Хуанга-Гильберта, для обработки ВР, полученных экспериментально 112
3.1. Анализ ВР, содержащего среднемесячные числа Вольфа .112
3.2. Анализ ВР, содержащего альфа-ритм ЭЭГ 127
3.3. Выводы 137
Глава 4. Программный инструмент NSDA для MATLAB 138
4.1. Общее описание функций программного инструмента NSDA 138
4.2. Графический интерфейс пользователя программного инструмента NSDA для MATLAB 139
4.3. Оптимизация программного кода NSDA 143
Заключение 146
Список литературы
- Вейвлет-анализ
- Исследование влияния значения ансамблевого числаNE на результаты декомпозиции методом ННТ
- Анализ ВР, содержащего альфа-ритм ЭЭГ
- Графический интерфейс пользователя программного инструмента NSDA для MATLAB
Вейвлет-анализ
Основные сведения о вейвлет-анализе изложены в книге И. Добеши «Десять лекций по вейвлетам» [9]. Здесь же приводятся только те моменты, которые необходимы для сравнения вейвлет-анализа с преобразованием Хуанга-Гильберта.
Термин «вейвлет» ввели Гроссман и Морле в середине 80-х годов в связи с анализом свойств сейсмических и акустических сигналов [69]. Вейвлет-пре-образование (wavelet transform) является инструментом, разбивающим данные на составляющие с разными частотами, каждая из которых затем изучается с разрешением, подходящим масштабу. Подобно оконному преобразованию Фурье, вейвлет-преобразование является инструментом для частотно-временной локализации особенностей сигнала [9]. Для непрерывного вейвлет-анализа (Continuous Wavelet Transform - CWT) используется выражение вейвлет-преобразования где а - задает масштабирование и называется параметром растяжения, Ь - соответствует временному сдвигу, и называется параметром положения,?/; - материнский вейвлет или базис вейвлет-преобразования. Непрерывно варьируя эти параметры а и Ь мы получим разные частотно-временные отображения для исходного сигнала.
Для выражения (1.6) справедливо обратное преобразование л а (1.7) может трактоваться с двух точек зрения: как способ восстановления исходной функции по ее вейвлет-преобразованию, или же как способ записи функции x(t) в виде суперпозиции вейвлетов фа . Стоит отметить, что выражение для вейвлетов (1.6) во многом схоже с оконным преобразованием Фурье (1.5), отличаясь только видом функции1 , которую мы можем определить любым образом с учетом отдельных ограничений [9].
Аналогично определяется прямое (1.8) и обратное (1.9) дискретное вей-влет-преобразование (Discrete Wavelet Transform - DWT): На практике для анализа ВР используется разновидность DWT, основанная на таком выборе параметров ао, bo и функции ф, чтобы функции фт п образовывали ортонормированный базис в L (Ж). В частности, для ао = 2 и bo = 1 существует целое семейство функций ф, обладающих хорошими свойствами частотно-временной локализации. Примером такого базисного вейвле-та может служить хорошо известный базис Хаара [47]
Ортонормированный базис для DWT позволяет проводить декомпозицию временного ряда применяя набор фильтров [9, 47]. Сначала сигнал x(t) пропускается через НЧ фильтр с импульсным откликом g и через ВЧ фильтр с импульсным откликом /г, в результате чего получаются свертки: у [п] = (х д)[п\= Е х М д[п-к]
Эти фильтры, связанные между собой, и называются квадратурными зеркальными фильтрами (QMF). В результате сверток получаются детализирующие коэффициенты cDi (после ВЧ-фильтра) и коэффициенты аппроксимации САІ (после НЧ-фильтра). Так как половина частотного диапазона сигнала была отфильтрована, то, согласно теореме Котельникова [15], отсчёты сигналов можно проредить в 2 раза (т.е. провести «децимацию» = downsampling):
Схема, иллюстрирующая данное разложение, представлена на рис. 1.2. На следующем шаге к детализирующим коэффициентам cDj применяется свертка, в результате чего выстраивается двоичное дерево фильтров (рис. 1.3). Это дерево имеет несколько уровней, его листья и узлы соответствуют пространствам с различной частотно-временной локализацией. Варьируя число каскадов, можно получать разное число детализирующих коэффициентов cDi и, меняя уровень декомпозиции исходного временного ряда, коэффициентов аппроксимации сА{, .
В настоящее время наиболее часто применяемым на практике средством вейвлет-анализа является подвид DWT, называемый вейвлетной пакетной НЧ прореживание Трехуровневое дерево вейвлет-декомпозиции [9] декомпозицией (Wavelet Packet Decomposition - WPD) [46]. Отметим главные отличия WPD от DWT: 1. повторная свертка применяется не только к детализирующим коэффициентам cDi, но и к коэффициентам аппроксимации сА{, в результате чего строится полное двоичное дерево декомпозиции (см. рис. 1.4); 2. из полного двоичного древа выбирается только одна ветвь декомпозиции, таким образом, чтобы конечная ветвь имела наименьшее значение энтропии среди всех ветвей двоичного древа. Более подробно этот алгоритм поиска оптимального базиса рассмотрен в [46]. Отметим, что результаты, полученные в области вейвлет-анализа (в частности, построение полной библиотеки всех известных дочерних вейвле-тов после их прохождения через гребенку двоичных фильтров), позволили существенно снизить вычислительные затраты WPD до сложности порядка 0(7Vlog2A ), что сделало этот метод весьма эффективным и распространенным для анализа временных рядов.
Исследование влияния значения ансамблевого числаNE на результаты декомпозиции методом ННТ
Из приведенного в Главе 1 описания ННТ следует, что основными параметрами модернизированного метода являются ансамблевое число NE (число итераций внешнего цикла добавления шума) и амплитуда добавляемого к сигналу шума (в отличие от автора метода [83] будем использовать соотношение сигнал/добавленный шум SNE)- На практике наиболее часто руководствуются рекомендациями автора, приведенными в [83], которые, однако, не имеют строгого теоретического обоснования. Принимая во внимание результаты критического рассмотрения существующих модификаций ННТ, представляется целесообразным провести самостоятельное исследование данной проблемы.
Выбор выделенной компоненты, наиболее близкой к исходному ВР, будем осуществлять в автоматическом режиме в соответствии с процедурой, описанной в [82]. 2.1.1. Исследование влияния значения ансамблевого числа NE на результаты декомпозиции методом ННТ
Рассмотрим применение ННТ к ВР (2.1) с параметрами / = 100 Гц, Т = 1, N = 2000, безразмерная частота сигнала fs = f/fd = 0.05, который ранее был проанализирован с помощью изначальной версии метода ННТ в [53]. Зависимость величины ошибки (2.2) от ансамблевого числа NE при различных значениях SNE приведена на рис. 2.1. Аналогичная зависимость для ВР с добавленным белым шумом (соотношение сигнал/шум в 10 дБ) приведена на рис. 2.2. Из рис. 2.1 и 2.2 видно, что:
Величина ошибки Error остается практически постоянной для ансамблевых чисел NE больших 100. Данный результат противоречит рекомендации из [83, с. 26] о том, что при уменьшении SNE ДЛЯ повышения точности декомпозиции требуется увеличивать ансамблевое число NE 2. Наиболее важным параметром, влияющим на точность декомпозиции методом ННТ, оказывается значение SNE (отметим, что Н. Хуанг рекомендует использовать значение SNE равным 14 дБ, однако, каких либо обоснований такого выбора не приводит).
Указанные зависимости от ансамблевого числа NE не меняются для различных значений частоты исходного сигнала /, длины отрезка Т и числа отсчетов N, поэтому результаты этих исследований в данной работе не приводятся.
Отметим, что хотя значение NE слабо влияет на точность декомпозиции ВР методом ННТ, данный параметр существенно влияет на скорость работы программной реализации метода ННТ. В связи с тем, что теоретическая сложность алгоритма EEMD неизвестна [53], были проведен вычислительный эксперимент, результаты которого представлены на рис. 2.3. 450
Таким образом, полученные результаты позволяют сделать следующие выводы: 1. Минимально возможное значение параметра NE равняется 100. 2. Рекомендации, данные в [83], в соответствии с которыми следует выбирать NE равным 200, обеспечивают выигрыш в точности не более 2%, поэтому далее мы ограничиваемся случаем NE равным 100. 3. Рекомендации из [83, с. 26] увеличивать NE при уменьшении соотношения сигнал/шум SNE не подтверждаются результатами вычислительного эксперимента. 2.1.2. Анализ зависимости погрешности декомпозиции ВР от SNE
Рассмотрим результаты применения метода ННТ к ВР (2.1), рассмотренному в предыдущем разделе. Зависимость величины ошибки (2.2) от SNE В диапазоне от -15 до 40 дБ для трех различных значений NE (100, 200, 400) приведена на рис. 2.4.
Причины возникновения максимумов в зависимости величины ошибки Error от SNE ПОНЯТНЫ ИЗ рис. 2.6, на котором отмечены номера значащих мод, отобранных автоматическим алгоритмом. Из рис. 2.6 видно, что максимумы обсуждаемой зависимости возникают при смешивании мод, номера которых отличаются друг от друга на единицу. Полученные результаты согласуются с известными фактами метода ННТ, что номер искомой компоненты определяется ее частотой и длиной ВР [53, 82]. Отметим, что эффект смешивания мод был обнаружен в [78] при анализе ВЧ сигналов, выявляющих 1 г
Смешивание 6 и 5 мод Соотношение сигнал/добавленный шум, дБ Номер значащей компоненты для ВР (2.1) с собственной частотой fs = 0.00625 кисту брюшной полости, однако, не получил своего объяснения. Определены значения абсцисс точек минимума зависимостей, представленных на рис. 2.5 при различных значениях безразмерной частоты сигнала fs (см. табл. 2.1).
Анализ ВР, содержащего альфа-ритм ЭЭГ
В данной главе рассмотрен опыт применения исследуемого метода ННТ для анализа некоторых известных временных рядов, в том числе: ВР, содержащий среднемесячные значения чисел Вольфа [76] и ВР, содержащий а-ритм электроэнцефалограммы (ЭЭГ) [51]. Оба ряда являются нестационарными согласно результатам тестов на стационарность ADF [48] и KPSS [64].
Число Вольфа R или относительное цюрихское число солнечных пятен является одним из главных индексов солнечной активности. Суточный индекс активности пятен R [73], определяется как R = k(10g + s), где S - число отдельных пятен, д - число групп пятен и к - фактор обсерватории (обычно меньший 1), который позволяет учесть суммарный вклад условий наблюдений, тип телескопа, и привести наблюдаемые величины к стандартным цюрихским числам.
В настоящее время накопление и распространение чисел Вольфа осуществляет Королевская обсерватория Бельгии (Royal Observatory of Belgium), на сайте [76] которой помещены таблицы и графики среднемесячных чисел Вольфа, начиная с 1749 г.
В связи с тем, что процессы, протекающие на Солнце, оказывают непосредственное влияние на геофизические и биологические процессы на Земле [5, 13, 35, 36, 52, 73], задача исследования временного ряда, описывающего динамику числа солнечных пятен, является актуальной. Данный вывод, в том числе подтверждается, продолжающимися в настоящее время активными поисками методов обработки данного ВР, содержащего среднемесячные числа Вольфа [19, 23, 38, 40, 41, 65]. Напомним, что основная трудность анализа обсуждаемого ВР обусловлена его нестационарностью, поэтому большинство известных сегодня методов анализа ВР, как стационарных так и нестационарных, применялись к данному ВР, в том числе методы спектрального оценивания (классические и непараметрические [13, 17, 38]), вейвлет-анализ [65], метод SSA [37, 52]. В тоже время универсальных методов, пригодных для анализа частотно-временных характеристик нестационарных временных рядов, на сегодняшний день не существует.
Необходимо также отметить особое место рассматриваемой задачи в истории развития методов спектрального оценивания [18], поскольку именно на ней проходили апробацию все известные на сегодняшний день методы частотного анализа временных сигналов (классические, параметрические, непараметрические). Большинство этих работ рассматривается в [18].
Ранее метод ННТ в его исходной модификации (EMD+HT) применялся к анализу ВР, содержащего среднемесячные значения чисел Вольфа за период с января 1749 г. по апрель 2010 г. (N = 3136) [40, 41]. В [41] автор использовал декомпозицию анализируемого ВР с помощью простого метода EMD с отсеиванием шумовых мод в соответствии с [82]. Далее он вычислил и визуализировал Гильбертов спектр, следуя [53]. Отметим, что в целом, полученные результаты которых для исследования обсуждаемого ВР использовался вейвлет-анализ. В связи с появлением новой модификации ННТ, включающей метод EEMD вместо метода EMD и метод DQ вместо Преобразования Гильберта, а также в связи с разработкой новой методики анализа ВР, описанной в разделе 2.4, представляется целесообразным провести самостоятельный анализ ВР, содержащего среднемесячные значения чисел Вольфа за период с января 1749 г. по октябрь 2012 г. (N = 3166) (см. рис. 3.1). 300 р 250 Т750 1800 1850 1900 1950 2000 Год
Для разделения выделенных компонент на информативные и шумовые были вычислены значения логарифмов средних приведенных энергий компонент от логарифмов соответствующих средних периодов (рис. 3.2). Здесь же изображены границы критерия значимости, соответствующие уровням доверительной вероятности 90%, 95% и 99%, вычисленные в соответствии с [82].
Из рис. 3.2 видно, что значения логарифма средней приведенной энергии выделенных компонент №1, JY52, JYe3, №4 и №10 находятся ниже доверительной границы, соответствующей вероятности 90%, что позволяет считать данные компоненты относящимися к шуму. Моды №5, №6, №7, №8 и №11 находятся выше доверительной границы критерия значимости, соответствующей вероятности 99%, поэтому будем считать их информативными. Отметим, что средний логарифм средней приведенной
Значения логарифма средней приведенной энергии выделенных компонент от логарифма их среднего периода, и границы критерия значимости энергии компоненты №9 оказывается близким к доверительной границы критерия значимости, а потому требуется ее дополнительный анализ. Также необходимо отметить, что наша методика позволяет отсеять 5 выделенных компонент, в то время как в [41] удалось отсеять только компоненту №1. Это объясняется тем, что в методе EEMD, в отличие от метода EMD, помимо шумовых компонент, присутствующих в самом ВР, возникают шумовые компоненты, обусловленные добавлением шума в 10.7 дБ при ансамблевой декомпозиции.
Компоненты анализируемого ВР, отнесенные к информативным, представлены на рис. 3.3 - 3.8. Здесь же для сравнения приведены соответствующие компоненты, выделенные с помощью вейвлет-анализа на основе WPD [46, 65]. Отметим, что полученные результаты согласуются с результатами, получаемыми с помощью метода SSA [37]. Это объясняется тем, что и EEMD метод и метод вейвлет-анализа (модифика 116 ция WPD) подобны частотным фильтрам, обсуждавшимся в главе 1 [49, 53]. Отличия компоненты №11, выделенной методом EEMD и соответствующей составляющей, выделенной при вейвлет-анализе, объясняется ошибкой декомпозиции метода ННТ, которая обсуждалась ранее в разделе 2.2.2.
Вычисление значений отсчетов функции, описывающей AM информативных компонент. В связи с тем, что безразмерная частота основной компоненты анализируемого ВР fs = 0.00758, для выделения функции, описывающей AM, будем использовать метод DQ. Вычисленные зависимости для каждой информативных компонент представлены на рис. 3.9 - 3.14 (зависимости a(t)). Здесь же для сравнения приведены соответствующие зависимости функций, описывающих AM, выделенные с помощью Преобразования Гильберта.
Вычисление зависимостей МЧ информативных компонент от времени. В связи с тем, что значения безразмерных частот всех информативных компонент оказываются меньше 0.1, для расчета значений МЧ использовался метод DQ. Вычисленные зависимости МЧ информативных компонент от времени представлены на рис. 3.9 - 3.14 (зависимости f(t)). Здесь же для сравнения приведены соответствующие зависимости значений МЧ информативных компонент, рассчитанные с помощью Преобразования Гильберта.
Также для сравнения степени близости результатов, полученных методом DQ и с помощью Преобразования Гильберта, были вычислены мгновенные значения периода каждой из информативных компонент, по которым далее были рассчитаны средние значения периодов ВР,
Графический интерфейс пользователя программного инструмента NSDA для MATLAB
Из рис. 3.23 видно, что энергия остаточного ряда оказывается существенно меньше энергии выделенной информативной компоненты. Оценочные вычисления показывают, что энергия остаточного ряда составляет не более 0.1% от энергии информативной компоненты.
Для сравнения результаты вычисления огибающей с помощью метода DQ, Преобразования Гильберта и метода SSA (L = 60, группируем компоненты №1, №4, №7, №10) представлены на рис. 3.24. Из рис. 3.24 видно, что огибающая, вычисленная методом SSA, значимо отличается от аналогичных зависимостей, вычисленных с помощью Преобразования Гильберта, в то время как аналогичная зависимость огибающей методом DQ практически совпадает с кривой, вычисленной с помощью Преобразования Гильберта. Следователь 134 но, использование разработанной методики для анализа ЭЭГ представляется более целесообразным, чем метода SSA.
Наличие остаточного ряда позволяет количественно оценить качество представления исходного ряда в виде суммы информативной компоненты и шума. Для этого достаточно проверить статистическую гипотезу о принадлежности распределения остаточного ВР нормальному распределению 7V(/i, о"), используя для этого, например, критерий согласия Колмогорова, позволяющий проверить сложную гипотезу Но [25, 26].
Гистограмма остаточного ВР и плотности нормального распределения с параметрами, оцененными по остаточному ряду, представлены на рис. 3.25. Вычисление критерия Колмогорова показало, что нулевая гипотеза Щ должна быть отвергнута на уровне значимости 95%.
Данный результат с нашей точки зрения объясняется тем, что не удалось добиться полной декомпозиции ВР из-за наличия шумов в исходном анализируемом ВР. Таким образом, результаты применения метода ННТ в задаче анализа ВР, содержащего среднемесячные значения чисел Вольфа [76], а также ВР, содержащего а-ритм электроэнцефалограммы (ЭЭГ) [51], показывают работоспособность предложенной методики анализа:
Результаты анализа ВР, содержащего среднемесячные значения чисел Вольфа [76], согласуются с результатами анализа данного ВР, полученными другими методами анализа НВР (вейвлет-анализ [23], метод SSA [37]), а также согласуются с современными представлениями о частотно-временных характеристиках данного ВР (цикл Швабе-Вольфа, минимум Дальтона, цикл Хале, цикл Гляйсберга, цикл Зюсса [73]).
Результаты анализа ВР, содержащего а-ритм ЭЭГ, показывают близость огибающей, выделенной методом DQ и Преобразованием Гильберта. Отличия огибающей, выделенной методом SSA, объясняются отсутствием обоснованного выбора группировки компонент, выделяемых данных методом [51].
Для возможности анализа ВР был разработан программный инструмент в математическом пакете MATLAB [66] фирмы The MathWorks Inc, который позволяет использовать следующие методы: метод SSA и модификации метода ННТ. Программный инструмент имеет свидетельство о регистрации программы №2014616281 от 19 июня 2014 г. Пакет MATLAB является хорошо известным и подробно документированным [2, 10, 22, 24] программным средством для инженерных и научных расчетов. MATLAB содержит встроенные средства для применения спектрального анализа и вейвлет-анализа, а также эффективные библиотечные функции построения эмпирических кривых, что в полной мере соответствует сформированной методике анализа НВР.
Общее описание функций программного инструмента NSDA Программа позволяет анализировать ВР с помощью разработанной методики. В основе данного инструмента лежит m-файл NSDA.m с описанием окна приложения NSDA.fig, реализующий весь графический интерфейс пользователя.
Данный файл также может порождать несколько специальных диалоговых окон: окно со значениями ВР (VarTable.fig), окно импорта данных из файла (ImportFile.fig), окно «О Программе» (About.fig). Данные файлы содержат набор вызываемых функций из соответствующих m-файлов, реализующих те или иные конкретные алгоритмы обработки ВР. Прототипы функций приведены в приложение А.
Программа NSDA имеет простой интерфейс, состоящий из главного окна (см. рис. Б.1 в приложении Б) и нескольких вспомогательных окон. Для загрузки исходных данных требуется требуется нажать кнопку «Открыть ...», после чего появится диалог с просьбой выбрать файл, содержащий ВР для анализа. Программа распознает следующие форматы файлов данных: .mat (MATLAB файлы), .dat .txt (бинарные или текстовые файлы данных), .xls .xlsx (файлы таблиц Excel). После выбора появится окно ImportFile (рис. Б.2), позволяющее предварительно визуально оценить загружаемые данные, а также выбрать разные переменные из файла из выпадающего списка ниже. Кроме того, эта форма обладает функционалом вариации импортирования временных рядов в зависимости от знака дробного делителя ( . или , ) (рис. Б.З), и возможности выбора одномерного вектора из многомерного входного массива (рис. Б.4).