Интерполяция исходного сигнала

Перед исследованием спектральных и нелинейных характеристик сигнала ВСР необходимо провести его интерполяцию, т. е. спроектировать на равномерную временную сетку с помощью команды

signal (:,2)=interpl (NN (:,1), NN(:,2), signal (:,1), method)', где NN (:, 1) — столбец времени регистрации текущего (R-R) интервала с неравномерной сеткой; NN (:,2) — значения длительности текущего (R-R) интервала в соответствующих точках неравномерной сетки; signal (:, 1) — столбец времени с равномерной сеткой. Этот столбец создается следующим образом:

T1=NN(,); T2=NN(end,l); signal (:,!)= (ТІ: delta: T2);

здесь delta — шаг равномерной сетки. В настоящей работе используется шаг, равный 100 мс.

В среде MATLAB существует несколько типов интерполяции [52]:

  • linear — линейная интерполяция;
  • nearest — интерполяция методом ближайшего соседа;
  • next — интерполяция методом следующего соседа;
  • previous — интерполяция методом предыдущего соседа;
  • pchip — кубическая интерполяция;
  • spline — интерполяция с применением кубических сплайнов.

Согласно [10] при анализе сигналов ВСР рекомендуется использовать сплайновую кубическую интерполяцию.

Полный алгоритм получения интерполированного сигнала выглядит следующим образом:

T1=NN(,V); T2=NN(end,1);signal(:,!)= (TT.100: T2);

signal (:,2)=interpl (NN(:,), NN(:,2), signal (:,1), ‘spline’);

Пример выполнения интерполяции сигнала ВСР показан на рис. 6: здесь в левой колонке приведен фрагмент исходного сигнала ВСР с неравномерным временным шагом, в правой — 30

результат интерполяции первых трех отсчетов исходного сигнала ВСР на равномерную сетку с шагом 100 мс.

0

956

956

944

1900

916

2816

980

3796

988

4784

940

5724

844

6568

788

7356

844

8200

852

9052

972

10024

988

Пример интерполяции сигнала ВСР

Рис. 6. Пример интерполяции сигнала ВСР

0

956

100

955

200

953

300

952

400

951

500

950

600

948

700

947

800

946

900

945

1000

943

1100

940

1200

937

1300

934

1400

931

1500

928

1600

925

1700

922

1800

919

1900

916

2000

923

Практическое задание 1

Обработка сигнала ВСР статистическими методами

1. Краткие теоретические сведения

Статистические методы применяются для непосредственной прямой количественной оценки ВР сигналов ВСР. Статистические характеристики BP (R-R) интервалов включают следующие показатели [8; 25]:

  • М — среднее значение (R-R) интервалов, которое обратно пропорционально частоте сердечных сокращений (ЧСС):

М = — YNN., Nh •’

где N— количество элементов временного ряда NN; NNj— і-й отсчет этого ряда;

• SDNN — стандартное отклонение или среднеквадратическое отклонение (СКО) (R-R) интервалов на интервале наблюдения:

SDNN = ^-j-?(AW,. “ ;

• RMSSD — квадратный корень из средней суммы квадратов разности величин последовательных отсчетов временного ряда NN:

RMSDD =

  • 1 ЛМ П05
  • -^BNN^-NN„)2

  • PNN50(%) — отношение оценки NN50к общему количеству пар последовательных отсчетов временного ряда MV;
  • CV — коэффициент вариации, отношение СКО к среднему значению, выраженное в процентах:

cy = SDNN 100%

м

Ниже представлены значения статистических характеристик ВСР в норме [9]:

ЧСС, уд/мин.......................55-80

SDNN, мс.............................40-80

RMSSD, мс..........................20-50

cv,%

  • 3-12
  • 2. Методические рекомендации

Для получения некоторых статистических характеристик достаточно знать одну соответствующую команду, например:

  • М=теап (Х(:,2))
  • SDNN=std(NN(:,2))

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

RMSDD= (sum ((NN (1: end-1,2)-NN (2: end,2)).*2)/length (NN))N),5;

Величина NN50 рассчитывается с применением цикла for:

NN50=0;

for i=l: length (NN)-1

if (NN (i,2)-NN (i+1,2)) >50

NN50=NN50+l;

end

end

Наконец, PNN50 вычисляется путем простого деления:

PNN50 = NN50/(length (NN)-lflOO;

Порядок выполнения задания:

  • • запустите среду MATLAB и создайте новый скрипт;
  • • произведите импорт данных согласно подглаве 2.5;
  • • произведите очищение сигнала от артефактов в соответствии с прил. 4;
  • • заполните форму данными, полученными при проведении исследований в этом практическом задании;

Результаты статистического анализа ВСР

ФС

М, мс

ЧСС, уд/мин

SDNN, мс

СК,%

RMSSD, мс

pNN50, %

ф

н

п

  • • сохраните изображения импортированного сигнала, как указано в п. 3 прил. 2 (см. с. 89);
  • • сравните полученные результаты со значениями, приведенными на с. 33;
  • • результаты исследований внесите в отчет.

Практическое задание 2

Обработка сигнала ВСР вариационными методами

1. Краткие теоретические сведения

Вариационная пульсометрия предназначена для изучения закона распределения (R-R) интервалов как случайных величин. Для этого строятся гистограммы распределения (R-R) интервалов в координатах: количество (R-R) интервалов К — длительность (R- R) интервалов т [ 11; 32]. Из всего многообразия гистограмм распределения ВСР наиболее типичными являются [12]:

  • • нормальная гистограмма (рис. 7, я), близкая по виду к кривым Гаусса, типична для здоровых людей в состоянии покоя;
  • • асимметричная (рис. 7,0 — указывает на нарушение стационарности процесса, наблюдается при переходных состояниях;
Типы гистограмм распределения

Рис. 7. Типы гистограмм распределения:

а — нормальная; б — асимметричная; в — эксцессивная; г — полимодальная

  • • эксцессивная (рис. 7, в) — характеризуется очень узким основанием и заостренной вершиной, регистрируется при выраженном стрессе, патологических состояниях;
  • • полимодальная (рис. 7, г) — обусловлена наличием несинусового ритма (мерцательная аритмия, экстрасистолия), а также множественными артефактами при регистрации ВСР.

Основные характеристики вариационной пульсограммы — мода Мо, амплитуда моды АМ0, вариационный размах VR.

Мода Mq— это наиболее часто встречающееся в исследуемом ВР значение (R-R) интервала. В физиологическом смысле — это наиболее вероятное значение кардиоритма. При нормальном распределении ряда Л/Омало отличается от математического ожидания М.

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

Вариационный размах (VR) отражает степень вариативности значений (R-R) интервалов в исследуемом ВР. Он вычисляется по разности максимального (R-R)MaKC и минимального (R-R)MI1H значений (R-R) интервалов и поэтому при аритмиях или артефактах может быть искажен. При вычислении VR следует отбрасывать крайние значения (R-R) интервалов, если они составляют менее 3 % от общего объема анализируемой выборки. В физиологическом смысле VR отражает активность парасимпатического отдела ВНС [4].

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

• индекс напряжения регуляторных систем

ИН =

ЛМ.

2M0VR

отражает степень централизации управления ритмом сердца и характеризует в основном активность симпатического отдела ВНС;

• индекс вегетативного равновесия

ив-4^

зависит от соотношения между активностью симпатического и парасимпатического отделов ВНС;

• вегетативный показатель ритма

ВПР =

1

M.VR

позволяет судить о парасимпатических сдвигах вегетативного баланса: чем меньше величина В ПР, тем больше вегетативный баланс смещен в сторону парасимпатической активности;

• показатель адекватности процессов регуляции

ПАПР =

отражает соответствие изменений функции автоматизма синусового узла в ответ на изменение симпатических регуляторных влияний на сердце [10; 19].

Ниже представлены характеристики вариационной пульсо-метрии в норме [9]:

Мо, мс.........................700-1100

АМ0, %..............................30-50

VR, мс...........................150-300

ИН, у.е..........................80-150

ИВР, у.е.......................100-350

ВПР, у.е............................3-10

ПАПР, у.е.......................25-50

2. Методические рекомендации

Первым шагом в построении гистограмм является вычисление временного интервала А между максимальным (R-R)MaKC и минимальным (R-R)MI1H значениями (R-R) интервалов. Далее определяем количество 50-миллисекундных интервалов л, укладывающихся в интервале А. Построение гистограмм производят в новом графическом окне. Для этого в командной строке введите:

figure; n=round ((max (NN (:,2))-min (NN (:,2)))/50);

hist(X(:,2),n);

Для нахождения моды MOi амплитуды моды АМ0 и вариационного размаха VR примените следующие алгоритмы:

[am, m]=hist (NN (:,2), п);

ат= 100*am/sum (ат);

[AMO, Nmax]=max (ат);

M0=m (Nmax);

Y=m (ат>3);

VR=Y (end)-Y (1);

Примечание. При расчете ИН, ИВР, ВПР, ПАПР значения моды А/о и вариационного размаха VR следует подставлять в секундах.

Порядок выполнения задания:

  • • запустите среду MATLAB и создайте новый скрипт;
  • • произведите импорт данных согласно подглаве 2.5;
  • • произведите очищение сигнала от артефактов в соответствии с прил. 4;
  • • заполните форму данными, полученными при проведении исследований в этом практическом задании;

Результаты пульсометрического анализа сигнала ВСР

ФС

Тип распределения

Мо, мс

VR, мс

АМ0, %

ИН

ИВР

ВПР

ПАПР

ф

н

п

  • • сохраните изображения гистограммы распределения, как указано в п. 3 прил. 2;
  • • сравните полученные результаты со значениями, приведенными на с. 37;
  • • результаты исследований внесите в отчет.

Практическое занятие 3

Оценка спектральных характеристик ВСР [фурье анализ)

1. Краткие теоретические сведения

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

При спектральном анализе сигналов ВСР выделяют следующие спектральные компоненты: высокочастотная (High Frequency — HF), низкочастотная (Low Frequency — LF), очень низкочастотная (Very Low Frequency — VLF) и ультранизкоча-стотные (Ultra Low Frequency — ULF).

Согласно российским и западным стандартам исследования спектральных характеристик ВСР проводят в следующих диапазонах частот [10; 60]:

  • • диапазон HF— от 0,4—0,15 Гц (2,5—6,5 с);
  • • диапазон LF — от 0,15—0,04 Гц (6,5—25 с);
  • • диапазон VLF — от 0,04—0,003 до Гц (25—333 с);
  • • диапазон ULF— частоты ниже 0,003 Гц.

При анализе кратковременных 5-минутных записей диапазон ULFwt анализируется.

Активность парасимпатического звена вегетативной нервной системы и активность автономного контура регуляции характеризуются мощностью высокочастотной составляющей спектра HF.

Мощность низкочастотной составляющей спектра LF, по мнению Р. М. Баевского, характеризует преимущественно состояние симпатического центра регуляции сосудистого тонуса.

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

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

При спектральном анализе ВСР применяют прямое дискретное преобразование Фурье в частотный спектр дискретного сигнала х(?-Д/) [23]:

X (I • Am) = • A/) e~IIM“bt,

к=0

где N — количество отсчетов сигнала х(к • At); к = / = 0, 1, (У— 1) — номер отсчета; At — интервал времени между отсчетами; T=(N— 1 )At — временной интервал анализируемого сиг-

л 2л „

нала; Дсо = — — шаг спектра в частотной области.

Для преобразования частотного спектра Х(1-Аы) в дискретный сигнал x(k-At) используют обратное преобразование Фурье:

х ? А/) = (/ • Д<о)е-" д" 4 "Ат.

/=о

Под спектральной мощностью понимается следующее выражение [10]:

|ЛГ(/.Дт)|2

'

Далее вычисляют:

• абсолютную спектральную мощность в каждом частотном диапазоне HF, LF и VLF.

^HFf

HF=^PI-

!=LHFi

(1)

VLF = ? />;

где LHPf, LHFi, LLFf, LLFi, LyLFf , LyLFi номера граничных отсчетов спектральной мощности Р, в каждом частотном диапазоне HF, LF и VLF,

  • • суммарную мощность спектра (Total Power — ТР), которая определяется как сумма мощностей в диапазонах HF, LFn VLF,
  • • нормированное значение мощности в каждом частотном диапазоне в процентах от суммарной мощности спектра TP-HFn, LFnnVLFn;

„ LF

• индекс вегетативного баланса — ;

Иг

• индекс централизации

nu=HF+LF

(2)

VLF

• индекс активации подкорковых нервных центров

ИАП =

LF VLF'

Ниже приведены значения интенсивностей нормированных значений спектральных составляющих ВСР в норме [9], %:

HFn..........................10-30

LFn..........................15-45

VLFn........................20-60

2. Методические рекомендации

Алгоритм применения преобразования Фурье для оценки спектральных мощностей сигнала ВСР представлен в прил. 5.

Примечание. Для проведения спектрального анализа необходимо избавиться от постоянной составляющей (линейного тренда). Для этого в MATLAB используется функция detrend. Алгоритм быстрого фурье-преобразования вызывается командой fft.

Порядок выполнения задания:

  • • запустите среду MATLAB и создайте новый скрипт;
  • • произведите импорт данных согласно подглаве 2.5;
  • • произведите очищение сигнала от артефактов в соответствии с прил. 4;
  • • заполните форму данными, полученными при проведении исследований в этом практическом задании;

Результаты фурье-анализа сигнала ВСР

ФС

HF1

мс

мс

VLF, мс2

77»

мс2

HFn, %

LFn, %

VLFn, %

LF/HF

ИЦ

ИАП

ф

н

п

  • • сохраните полученные в результате исследований изображения, как указано в п. 3 прил. 2;
  • • сравните полученные результаты со значениями, приведенными на с. 42;
  • • результаты исследований внесите в отчет.

Практическое задание 4

Оценка спектральных характеристик ВСР [вейвлет-анализ]

1. Краткие теоретические сведения

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

Уравнение для непрерывного интегрального вейвлет-пре-образования анализируемого сигнала s(t) имеет следующий вид [20]:

W(a,b} = -= fs'(z)• —— ^,

daJ < a ) где a — масштабирующий параметр; b — параметр сдвига; у — анализирующий базисный вейвлет; s(t) — исходный сигнал.

В среде MATLAB связь масштабирующего параметра а с частотой [43] выражается формулой

delta ? f ’

где fc центральная частота вейвлета, вызываемая функцией centfrq: delta — время дискретизации исследуемого сигнала;/— анализируемая частота.

Результаты оценок применения вейвлет-преобразования для конкретного вида изучаемого сигнала могут зависеть от выбора базисной функции. Для оценок спектральных характеристик ВСР рекомендуется использовать следующие базисные вейвлеты [2]:

  • • Гаусса (gausN, где N = {5, 7, 8});
  • • Добеши (dbN, где N = {5, 7, 8, 9, 10});
  • • симлеты (symN, где N = {5, 7, 8});
  • • койфлеты (coifN, где N = {3, 4, 5});
  • • Морле (morl).

LF^

HFm

Информационные особенности вейвлет-преобразования анализируются по данным ВР индекса вегетативного баланса (/) , где LFm и HFm вейвлет-преобразования интен-

сивностей периодических составляющих ВСР в диапазонах ча

стот LFn Я/7соответственно [2|. На рис. 8 приведены графики

Функционал F

ла F

LF

(/) не является гладкой функцией,

а параметры его «выбросов» (локальных дисфункций) изменяются при функциональных нагрузках.

Для идентификации локальных дисфункций AF функциона

лу?

гг ЛМП применяется решающее правило [341:

HFm 'J где Д — порог принятия решения о наличии дисфункций функ

ционала F

LF^

HFm

в норме, Д = 10.

В качестве параметров локальных дисфункций AF функцио

нала F

LF*

HFm

(')

применяют количество дисфункций N, мак

симальное значение Лмакс их амплитуды и интенсивность 4ntense на интервале наблюдения.

2. Методические рекомендации

Алгоритм применения вейвлет-преобразования представлен в прил. 6. Для вычисления значений вейвлетной плотности мощности в диапазонах частот HF, LF и VLF формируют выборку вейвлет-коэффициентов W (а, Ь) согласно формуле (1) с использованием функции cwt. Затем производят ее оценки по формуле

6=/ _1 Y

../(/)= , b) • Дat ? at 2 ,

<‘=J 7

где jwf — индексы, соответствующие нижним и верхним границам исследуемого диапазона частот.

Примечание. В качестве базисного вейвлета в прил. 6 используется вейвлет Морле.

Порядок выполнения задания:

  • • запустите среду MATLAB и создайте новый скрипт;
  • • произведите импорт данных согласно подглаве 2.5;
  • • произведите очищение сигнала артефактов в соответствии с прил. 4;

• заполните форму данными, полученными при проведении исследований в этом практическом задании;

Результаты вейвлет-анализа сигнала ВСР

ФС

HF1 мс

МС

VLF, мс2

ГР мс2

HFn, %

LFn, %

VLFn, %

N

А •^макс

Д -^intense

ф

н

п

  • • сохраните изображения графика индекса вегетативного баланса, как указано в п. 3 прил. 2;
  • • результаты исследований внесите в отчет.

Практическое задание 5

Построение скаттерограммы сигнала ВСР

1. Краткие теоретические сведения

Скаттерограмма является компактным способом изображения последовательности (R-R) интервалов с помощью метода корреляционной ритмографии. Сущность метода корреляционной ритмографии заключается в графическом отображении последовательных пар (R-R) интервалов (предыдущего и последующего) в двухмерной координатной плоскости [35]. При этом по оси абсцисс откладывается (R-R)n — п-й (R-R) интервал, а по оси ординат откладывается (7?-7?)я+1(п+1 )-й (R-R) интервал.

График и область точек, полученных таким образом, называются корреляционной ритмограммой, или скаттерограммой. При построении скаттерограммы образуется совокупность точек (облако точек или пятна Пуанкаре), центр которых располагается на биссектрисе. Обычно скаттерограмма имеет форму эллипса. Расстояние от центра до начала осей координат соответствует наиболее ожидаемой длительности сердечного цикла (моде Мо). Величина отклонения точки от биссектрисы показывает, насколько п-й (R-R) интервал короче или длиннее (л+1)-го (R-R) интервала.

Количественные показатели скаттерограммы:

  • • длина L скаттерограммы определяется как размер длинной оси эллипса, выраженный в миллисекундах;
  • • ширина w скаттерограммы определяется как размер перпендикуляра к длинной оси, проведенный через ее середину, выраженный в миллисекундах;
  • • площадь Sскаттерограммы вычисляется по формуле оценки площади эллипса

с_ tiLw

4

На рис. 9 представлен пример скаттерограммы сигнала ВСР, на котором изображены аппроксимирующий эллипс и показатели ? и tv.

Скаттерограмма сигнала ВСР

Рис. 9. Скаттерограмма сигнала ВСР

2. Методические рекомендации

В среде MATLAB имеется специальная команда для построения скаттерограмм — scatter. Следующие команды дают возможность построить скаттерограмму и биссектрису угла:

figure; scatter (NN (1: end-1,2), NN (2: end,2)); hold on;

plot (min (NN (:, 2)): max (NN (:,2)), min (NN (:, 2)): max (NN

‘k’, 'Line Width ’,2); grid on;

В прил. 7 описывается алгоритм оконтуривания точек скаттерограммы аппроксимирующим эллипсом, основанный на ме тоде наименьших квадратов. Этот алгоритм позволяет рассчитать длину ?, ширину w и площадь эллипса S.

Порядок выполнения задания:

  • • запустите среду MATLAB и создайте новый скрипт;
  • • произведите импорт данных согласно подглаве 2.5;
  • • произведите очищение сигнала от артефактов в соответствии с прил. 4;
  • • заполните форму данными, полученными при проведении исследований в этом практическом задании;

Результаты оценки скаттерограммы сигнала ВСР

ФС

L, мс

W, мс

S, мс2

ф

н

п

  • • сохраните изображения скаттерограммы, как указано в п. 3 прил. 2;
  • • результаты исследований внесите в отчет.

Практическое задание 6

Комплексная оценка функционального состояния

1. Краткие теоретические сведения

Одним из методов комплексной оценки вариабельности сердечного ритма является вычисление показателя активности регуляторных систем (ПАРС). Он вычисляется по специальному алгоритму [24]

ПАРС = ?|ЛГ,|.

/=1

В табл. 2 представлены соответствия между оценками, полученными статистическими, вариационными и спектральными методами обработки сигналов ВСР, и весовыми коэффициентами индекса ПАРС, используемыми в формуле (2).

Таблица 2

Связь оценок ВСР и весовых коэффициентов ПАРС [24]

Название весового коэффициента

Значение весового коэффициента

Оценки временного ряда ВСР

Суммарный эффект регуляции

+2 + 1 0 -1 -2

М, мс <660 <800

>1000 >1200

Функция автоматизма N2

+2 + 1 0 -1 -2

SDNN,mc <20 >100

>100

VR/M <0,10, >0,30, (0,10-0,30) >0,45, >0,60

cv,% <2,0 >8,0

>8,0

Окончание табл. 2

Название весового коэффициента

Значение весового коэффициента

Оценки временного ряда ВСР

Вегетативный гомеостаз А3

+2 + 1 0 -1 -2

VR, мс <60 <150

>300 >500

АМ0, % >80 >50

<30 <15

ИН, у.е.

>500

>200

<50

<25

Устойчивость регуляции N4

+2 0 -2

cv,% <3,0

>6,0

Активность подкорковых нервных центров n5

+2 + 1

0 -1 -2

VLF, % >70 >60

<40 <20

LF,% >25

HF, % <5 <20

>30 >40

Оценка ПАРС позволяет дифференцировать различные степени напряжения регуляторных систем и оценивать адаптационные возможности организма. Оценка функционального состояния человека формируется в соответствии с данными, приведенными ниже [10; 24]:

норма.................................................................0—2

функциональное напряжение:

умеренное.....................................................3—4

выраженное..................................................5—6

резкое...........................................................7—8

перенапряжение регуляторных систем............8—9

истощение регуляторных систем.........................10

2. Методические рекомендации

Значения весовых коэффициентов определяются оценками сигнала ВСР, полученными в результате выполнения прак тических заданий 1—3 (табл. 2), столбец «Оценки временного ряда ВСР». Для коэффициентов N2, N3) N5 достаточно выполнения двух условий из трех.

Порядок выполнения задания:

  • • проведите вычисление ПАРС согласно рекомендациям этого практического задания, используя данные из форм на с. 34, 39, 43;
  • • на основании полученных выше значений ПАРС оцените функциональное состояние человека;
  • • результаты исследований внесите в отчет.

Практическое занятие 7

Построение странного аттрактора

1. Краткие теоретические сведения

Одним из самых простых и наглядных способов оценки хаотического поведения системы является построение аттрактора фазовой траектории движения системы. Данный метод применяется для анализа статистических и фрактальных свойств аттракторов фазовых траекторий [51].

В фазовом пространстве колебательной системы предельное множество точек, притягивающее фазовые траектории, называется аттрактором (от глагола to attract — притягивать). Аттракторы могут быть простыми, регулярными, с установившимся режимом и так называемые «странные» аттракторы с непериодическими колебаниями. Странный аттрактор в отличие от регулярного не является многообразием (т. е. не является кривой или поверхностью), его геометрическое устройство очень сложно, а его структура фрактальна. Такие аттракторы обладают геометрической (масштабной) инвариантностью, или, другими словами, скейлинговой структурой.

Как показано в [35], критериями «странного» аттрактора являются:

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

Способ определения размерности аттракторов описан в работах [29; 54].

Для таких сложных систем, как живые организмы, по измеренному ВР одной наблюдаемой динамической переменной можно сконструировать странный аттрактор. Согласно теореме Такенса, основные свойства аттрактора будут такими же, как у исследуемого объекта, и на основе подобия можно определить его характеристики и попытаться построить математическую модель исходной системы [44]. Моделирование не является целью данного учебного пособия, его задача — показать возможности методов нелинейной динамики оценивать информационно значимые параметры сложных систем.

В общем случае фазовое пространство является TV-мерным. Однако для практических случаев можно считать пространство отображения трехмерным, т.е. N= 3. На рис. 10 представлена трехмерная визуализация переменной (сигнала ВСР), ее первой и второй производной, что дает хорошее представление о динамике процесса [13].

Аттрактор в трехмерном фазовом пространстве

Рис. 10. Аттрактор в трехмерном фазовом пространстве

Анализируя вид аттрактора ВСР, мы анализируем интегральные целостные процессы, базирующиеся на взаимодействиях между отдельными компонентами ВСР или же между этими компонентами и сопряженными изменениями в других органах и системах, протекающих как феномен самоорганизации или синергии [48].

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

Первый шаг в общем определении размерности аттрактора состоит во введении оценки размерности Хаусдорфа — размерности подмножества в метрическом пространстве [28]. Размерность Хаусдорфа D странного аттрактора можно вычислить следующим образом:

In (А (5))

Z) = lim v 7 ,

In (1/5)

где N (б) — количество ячеек с размером 6, покрывающих кривую.

2. Методические рекомендации

Для построения трехмерного странного аттрактора необходимо знать сигнал и две его первые производные. В MATLAB операция взятия производной сводится к применению функции diff к анализируемому массиву данных. При этом длина вектора производной меньше длины исходного массива данных на 1.

Для построения трехмерного графика используется команда plot3. Алгоритм построения трехмерного аттрактора сигнала ВСР выглядит следующим образом:

Dl=diff(NN(:,2)); D2=diff(Dl);

plot3 (NN(1: end-2,2), DI (1: end-1), D2,’b’)

Алгоритм нахождения размерности Хаусдорфа D приведен в прил. 8.

Примечание. Для более корректной работы алгоритма применена перенормировка исходного ВР.

Порядок выполнения задания:

  • • запустите среду MATLAB и создайте новый скрипт;
  • • произведите импорт данных согласно подглаве 2.5;
  • • произведите очищение сигнала от артефактов в соответствии с прил. 4;
  • • выполните построение странного аттрактора согласно рекомендациям п. 2 настоящего задания;
  • • сохраните изображения странного аттрактора, как указано в п. 3 прил. 2;
  • • проведите расчет размерности Хаусдорфа D согласно рекомендациям п. 2 настоящего задания;
  • • результаты исследований внесите в отчет.

Практическое задание 8

Оценка старшего показателя Ляпунова

1. Краткие теоретические сведения

Важной характеристикой хаотических пульсаций являются показатели Ляпунова, которые определяют скорость экспоненциального роста малых возмущений (разбегания двух изначально близких траекторий аттрактора) [6]. Нарис. 11 представлен пример разбегания двух изначально близких траекторий аттрактора для фрагмента реального сигнала ВСР.

Разбегание двух изначально близких траекторий аттрактора реального сигнала ВСР

Рис. 11. Разбегание двух изначально близких траекторий аттрактора реального сигнала ВСР

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

Расстояние между ними можно записать как

|Х,(О)-^(О)| = 8о«1.

Стечением времени расстояние будет меняться. Обозначим это расстояние через время t как 8 (z) = |(/) - Х2 (/)|. Тогда старший показатель Ляпунова X можно определить по соотношению

8(/) = 80 ехр(Х/).

Стоит отметить, что величина, обратная показателю Ляпунова, в общем случае не может рассматриваться как средний горизонт предсказуемости или время предсказуемости процесса [15].

В случае хаотического поведения старший показатель Ляпунова всегда положителен (Z > 0). Более того, используя этот критерий, нетрудно определить величину неупорядоченности системы, т. е. степень ее хаотичности: чем больше старший показатель, тем глубже хаос. Это означает, что система мгновенно забывает то, что происходило с ней ранее, всякая детерминированность в движении такой системы отсутствует, совпадение траекторий или даже просто их сближение исключено.

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

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

Полный спектр показателей Ляпунова получается при определении темпов разбегания траекторий в разных направлениях в A-мерном фазовом пространстве. Сумма всех положительных показателей Ляпунова часто называется энтропией Колмогорова [28].

Алгоритмы оценивания старшего показателя Ляпунова по одной реализации исследуемого процесса основаны на ис пользовании того факта, что с течением времени расстояние между двумя траекториями увеличивается со скоростью, определяемой X [41].

Зафиксируем один отсчет х(/) из наблюдаемой реализации х. Все отсчеты данной реализации x(i), для которых выполняется условие |х(/)-х(/)|<є, будем считать началами е-близких траекторий. Исходную и соседние траектории сформируем путем последовательной записи отсчетов, начиная с t и і соответственно. Расстояние между данной и соседней траекториями через интервал времени т после начала сравнения определяется выражением

дІ8і(х(/),х(/);т) = |х(ґ+т)-х(/+т)|.

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

Zdist(xW>40;T) •

Тогда крутизна кривой S (т) на линейном участке определяется показателем Ляпунова X.

2. Методические рекомендации

В прил. 9 представлена реализация алгоритма нахождения старшего показателя Ляпунова.

Примечание. После выполнения алгоритма будет получен график зависимости S (т) (вектор vectau — набор интервалов т, вектор stat — значение статистики для каждого интервала). Предлагается оценить номер точки L, на которой заканчивается линейный участок. Старший показатель Ляпунова находится из угла наклона этой прямой, т. е. после выполнения команд

C=polyfit (vectau (1: L), stat (1: L)1);

Lambda=C (1);

Переменная closejtr содержит TV ячеек. Каждая /-я ячейка содержит номера всех точек, находящихся в є окрестности от 1-й точки.

Порядок выполнения задания:

  • • запустите среду MATLAB и создайте новый скрипт;
  • • произведите импорт данных согласно подглаве 2.5;
  • • произведите очищение сигнала от артефактов в соответствии с прил. 4;
  • • выполните расчет показателя Ляпунова согласно рекомендациям п. 2 настоящего практического задания;
  • • сохраните изображения зависимости статистики, как указано в п. 3 прил. 2;
  • • результаты исследований внесите в отчет.

Практическое задание 9

Расчеты спектра размерностей Реньи

1. Краткие теоретические сведения

Отображение динамических процессов в фазовом пространстве позволяет ввести вероятностное описание анализируемого процесса. Поскольку вся реализация отображается в некоторый объем фазового пространства, то при разбиении всей области фазового пространства на ячейки подходящего размера каждая из них может быть охарактеризована своей заселенностью — целым ЧИСЛОМ точек Пі (б) в ячейке с номером і объемом 6J. Величина d характеризует геометрическую размерность фазового пространства, а величина б является линейным размером ячейки. В этом случае относительная заселенность ячейки, т.е. вероятность того, что в z-й ячейке находится хотя бы одна точка отображающего множества, представляет собой предел [ 14]:

пі

W

где п — число точек в z-й ячейке; і = 1, 2, 3,...; N (б) — суммарное число занятых ячеек, в которых есть хотя бы одна точка. При этом выполняется условие нормировки вероятности

ZV(8)

1>(8) = 1.

Удобной мерой оценки неопределенности состояния динамической системы является спектр энтропии Реньи, связанный с f степенями Р-.

М8)=Г

^(8)

In .

/=о J

В дополнение к размерности Хаусдорфа D для характеристики неоднородности статической структуры аттрактора вводится спектр обобщенных размерностей Реньи

df =-lim

(3)

В частном случае при f 0 d0 сводится к размерности Хаусдорфа и информации о статистических свойствах отображения не несет:

ЛГ(8) !

dn = - lim In У--- = - lim

s^° I yZo In 5 I 8->°

В случае анализа сигналов ВСР геометрическая размерность, являющаяся показателем структурной неоднородности

системы, отражает структурную сложность процессов регуляции сердца [3]. При f 1 после раскрытия неопределенности по правилу Лопиталя из формулы (3):

d,=- lim

  • 1 8->0
  • *.(8)

In 5

*(8)

где Кх (8) = ~^Pj In Pj — совпадает с формулой информации Шс-7=0

нона[1 ].

Величина Кх (8) определяет приращение информации при всех и тестных вероятностях попадания точек аттрактора в z-ю ячейку и при условии прохождения траектории через эту же ячейку. Поэтому величина d{ показывает, как возрастает получаемая информация при стремлении размеров ячейки к нулю. В нашем случае информационная размерность отражает информационную сложность характеристик исследуемого сигнала.

Рассмотрим еще один частный случай обобщенной размерности dq, имеющий важное прикладное значение, при:/= 2. Выражение для d2 имеет вид

= lim

8-»0

In 5

где К2 (5) = In

— корреляционная энтропия [31].

При анализе сигналов ВСР корреляционная размерность отражает динамическую неоднородность функционирования сердца.

Как известно, временные ряды могут представлять собой фрактальные объекты. Под фракталом понимают множество, части которого подобны целому в некотором смысле [45]. Основной характеристикой таких объектов является фрактальная размерность D, которая принимает дробные значения. Для одномерных ВР /)є[1...2]. Мандельброт установил связь между фрактальной размерностью D ряда и его показателем Херста Н:

D = 2-H .

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

Оценка корреляционной размерности через корреляционный интеграл

Корреляционный интеграл [5] определяется по формуле где є — размер разрешающей ячейки; N— количество точек аттрактора; 0 () — функция Хевисайда; |х, -ху| — абсолютное расстояние между і-й иj-й его точками в /77-мерном пространстве.

(4)

По сути говоря, Ce(e,AQ — зависимость количества точек аттрактора в /77-мерном пространстве, расстояние между которыми меньше 8, от размера разрешающей ячейки, отнесенная к полному количеству пар точек, т. е. -N2 (в знаменателе формулы стоит N (N — 1), поскольку поставлено условие і* j).

Полученные зависимости Се(є,7У) откладываются в двойном логарифмическом масштабе на плоскости (для наглядности лучше брать по основанию 10). Затем выделяют линейные участки отложенных кривых и по методу наименьших квадратов производят поиск аппроксимирующих их прямых. Далее, для всех полученных кривых Се (є, У) вычисляют первую производную от аппроксимирующих их прямых dc (фрактальные размерности) и откладывают ее как функцию от т. Теоретически производная dc определяется из предела [5]:

dr = lim lim

V ?—>0 7V->00

JlgCe(e,7V) digs

Данный алгоритм вычисления dc связан с тем, что при сравнительно малых значениях є должен соблюдаться степенной закон

Ce(e,#)~?rfc,

где dc корреляционная размерность. Поскольку корреляционная размерность идет под индексом f= 2 в спектре Реньи df , то она является нижней оценкой размерности Хаусдорфа-Безиковича (которая идет под индексом f= 0), так как спектр Реньи является ниспадающим с ростом индекса/. Также нужно учитывать, что в численных экспериментах N всегда конечно и оба предела є -> 0 и N -> оо являются бессмысленными.

На полученном графике dc (m) ищут точку тс, когда зависимость dc (/и) достигнет насыщения. Значение этой точки тс будет соответствовать независимой оценке размерности про странства вложения, а соответствующее ей значение dc соответствует корреляционной размерности исследуемого псевдо-фазового аттрактора, восстановленного из исследуемого ряда.

На рис. 12 представлен пример нахождения корреляционной размерности для реального сигнала ВСР.

Нахождение корреляционной размерности

Рис. 12. Нахождение корреляционной размерности

Количественную оценку точности нахождения dc предложили Экман и Рюэль [44]. Данное условие выполняется при соблюдении неравенства

dc<2gN ,

где dc вычисленная корреляционная размерность; N — длина исследуемого ряда.

Методические рекомендации

Прямая оценка размерностей Реньи

В прил. 10 представлена реализация алгоритма прямого подсчета трех основных размерностей Реньи — геометрическая (d0), информационная (dl) и корреляционная (d2) с использованием алгоритма boxcounting.

Примечание. Для более корректной работы алгоритма применена перенормировка исходного ВР.

Оценка корреляционной размерности

через корреляционный интеграл

Алгоритм расчета корреляционного интеграла и размерности dc ВР ВСР представлены в прил. 11.

Примечание. В MATLAB имеется встроенная функция Хевисайда: heaviside. Однако использование ее для прямой подстановки в формулу (4) не целесообразно для длительных ВР из-за долгого времени вычисления. Поэтому в прил. 11 используется комбинация других быстрых математических операций, дающих аналогичный результат и позволяющих использовать операцию определения знака sign и операцию сложения:

Corrlnt (е, т) = Corrlnt (е, т) + N*0.5*...

(sign (Epsilon (e)-norm (data (i,:)-data (j,:),!))+!);

В прил. 11 значения параметра е выбираются исходя из значения стандартного отклонения SDNN ряда и варьируются в пределах от 0.1* SDNN до 1* SDNN.

Порядок выполнения задания:

  • • запустите среду MATLAB и создайте новый скрипт;
  • • произведите импорт данных согласно подглаве 2.5;
  • • произведите очищение сигнала от артефактов в соответствии с прил. 4;
  • • реализуйте алгоритмы для расчета размерностей Реньи согласно рекомендациям настоящего практического задания;
  • • результаты исследований внесите в отчет.

Практическое задание 10

Вычисление аппроксимированной энтропии

1. Краткие теоретические сведения

Другая важная характеристика нелинейной динамики — энтропия — величина, определяющая степень предсказуемости системы или, иначе говоря, определяющая относительную скорость накопления информации в системе с течением времени. Если энтропия достигает нуля, то информационная полнота системы достигает максимума и система становится полностью предсказуемой [29].

Траектория движения динамической системы, измеренная равноотстоящими отсчетами через интервалы времени т в TV-мерном фазовом пространстве, разделенном на ячейки размером д, представляется вектором х(/) = [jq (/),...,xd (/)].

Если Pj — вероятность нахождения точки траектории движения динамической системы в ячейке і, то энтропию Колмогорова в этом случае определяют по формуле Шеннона:

^ = -Хд1пй.

/

Для странного аттрактора энтропия положительна, но имеет конечное значение, и это значение является количественной характеристикой степени хаотичности системы. Практическое вычисление К2 энтропии Колмогорова выполняется на основе обобщенного корреляционного интеграла [1].

Энтропия представляет собой простой показатель для оценки общей сложности и предсказуемости временных рядов. Для практической реализации расчета энтропии при анализе ограниченных и зашумленных временных рядов (R-R) интервалов в ряде работ используется алгоритм расчета аппроксимированной (АрЕп), предложенный S. М. Pincus [62].

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

ApEn(m,r,N) = Y In [С" (г)/С"'*1 (г)],

где С/” (г) и С/”+1 (г) определяются суммами

і N-m+l

g eH* W-* wi)>

Оценка Ар Ей определяет количественную вероятность того, что точки ВР, которые изначально близки (|x(z)-x(j) < f), остаются так же близки и для последующих сравнений. Высокие значения АрЕп указывают на высокую неравномерность и сложность распределения ВР.

Методические рекомендации

При расчете АрЕп обычно используется значение т = 2, что позволяет интерпретировать АрЕп как различие вероятностей обнаружить сходные векторы при размерностях вложения т = 2 и т = 3 соответственно с толерантностью, задаваемой параметром г = 0.2*SDNN, где SDNN — стандартное отклонение.

Описанный выше алгоритм реализован в прил. 12.

Порядок выполнения задания:

  • • запустите среду MATLAB и создайте новый скрипт;
  • • произведите импорт данных согласно подглаве 2.5;
  • • произведите очищение сигнала от артефактов в соответствии с прил. 4;
  • • выполните расчет аппроксимированной энтропии согласно рекомендациям п. 2 настоящего практического задания;
  • • результаты исследований внесите в отчет.

Практическое задание 11 Расчеты показателя Херста

1. Краткие теоретические сведения

Понятие обобщенного броуновского движения введено Мандельбротом через обобщения случайной функции на случай фрактальной гауссовой функции Вн (/) с параметром Не [0,1], называемым показателем Херста [45].

Обобщенный броуновский процесс имеет нулевое среднее приращение:

<ДВя) = (Вя(ґ)-Дяо)) = О.

Из уравнения (1) следует, что дисперсия приращений S(t— /0) может быть записана в виде

s(/-roMMO-Mo)]2H2k-'oi>

где о — положительная константа. Тогда Ви (/) имеет гауссово распределение [20]: где /(, t2 отчеты ВР в моменты времени 1 и 2.

Р(ДВ„<г) =

  • 1
  • 72(г2-/1)И

Для анализа ВР традиционно применяется метод Херста, называемый также классическим методом нормированного размаха или А/З’-методом [5].

Для имеющегося временного ряда %(/) вычисляется среднее значение (^(0) на интервале времени т, имеющем ту же размерность, что и время t:

т /=1

Затем рассчитывается зависимость накопленного отклонения X (t, т) на интервале времени т:

W=1

По накопленному отклонению вычисляется функция абсолютного размаха R:

R (т) = max X (t, т) - min X (t, т).

Размах зависит от длины интервала т и может расти с ее увеличением. Далее вычисляется зависимость безразмерной функции R/S от длины временного интервала т делением R на стандартное отклонение S ряда ? (/):

V ~ и=1

Херст по результатам исследования многих природных процессов установил эмпирическую связь между нормированным размахом R/S и длиной интервала т через показатель Н [30]:

R/S = (t/2)h ,

где Н может принимать значения от 0 до 1. Это наблюдение Херста интересно потому, что если отсутствует долговременная статистическая зависимость (случайный ряд), данное отношение должно асимптотически стремиться к т = 1/2 при стремлении длины выборки к бесконечности. Значения же Н > 0,5 характеризуют сохранение тенденции ряда к росту или убыванию как в прошлом, так и в будущем. Если Н < 0,5 — это означает склонность ряда к смене тенденции: рост сменяется убыванием и наоборот [30].

2. Методические рекомендации

В среде MATLAB при анализе временных рядов используются несколько методов оценки показателя Херста с помощью функции wfbmesti (X), которые основаны на вейвлет-преобразованиях: в первых двух методах используется дискретное вейвлет-преобразование (ДВП) по базису Хаара для первого и второго уровней детальности разложения [56]. Третий основан на оценке линейной регрессии уровней детальности разложения ДВП в двойных логарифмических координатах [57].

Для этих же целей лучше подходит метод накопленной дисперсии. Дисперсия имеет вид

var|(A'(/,)-jr(/1))2| = a2|z,-Z1|2",

где Н — показатель Херста, определяется как угловой коэффициент из отношения logcj,.m(A?)=c + /7Tog|5|, здесь а,ТО5 (АуТ) — среднеквадратичное отклонение приращений АТ, соответствующих временному интервалу 5; с — константа. При использовании этого метода требуется выполнение условия нормального распределения первых разностей (приращений) временного ряда [26].

В прил. 13 приведен алгоритм реализации данного метода оценки показателя Херста.

Порядок выполнения задания:

  • • запустите среду MATLAB и создайте новый скрипт;
  • • произведите импорт данных согласно подглаве 2.5;
  • • произведите очищение сигнала от артефактов в соответствии с прил. 4;
  • • выполните оценку показателя Херста с помощью функции wfbmesti;
  • • реализуйте алгоритм нахождения показателя Херста методом накопленной дисперсии согласно рекомендациям п. 2 настоящего практического задания;
  • • сохраните изображения полученных в ходе расчетов графиков, как указано в п. 3 прил. 2;
  • • результаты исследований внесите в отчет.

Практическое задание 12

Применение метода максимумов модулей коэффициентов вейвлет-преобразования

1. Краткие теоретические сведения

Введем обобщенную статистическую сумму Z(s,g) [5], характеризуемую показателем степени q, который в общем случае может принимать любые значения в интервале -оо < q < +со 5 следующим образом:

*(4

^(м)=2>’(г),

/=1

где Рі (є) — вероятность нахождения произвольной точки исследуемого ВР в /-й ячейке размером е.

Распределение плотности вероятностей спектра обобщенных размерностей Dq определяется с помощью соотношения

Функция т (q) имеет вид

lnZ(e,^) т(<7) = lim-------.

In є

Показатель т(

Традиционно рассматриваются ряд моментов q в диапазоне значений от —5 до 5 с шагом Дд = 0,1 [59; 61]. Варьирование показателя q позволяет рассматривать различные масштабы флуктуации исходного сигнала: при q < 0 основной вклад в стати стическую сумму вносят флуктуации малого порядка, при q > О флуктуации больших масштабов вносят больший вклад в статистическую сумму.

Алгоритм метода WTMM представлен на рис. 13 в виде структурной блок-схемы. В ходе работы алгоритма используется ВР, интерполированный на равномерную сетку.

Исходный сигнал

Спектрограмма

Локальные максимумы ММ(а,Ь)

Статистическая сумма

Z(a,q) = ^(MM(a,t)q)

Скейлинговый показатель

Z(a,q) ~ax(q)

Мультифрактальный спектр

а = ;?)(а) = qa- т(

dq

Мультифрактальный показатель Херста #2= <4=2

Структурная блок-схема метода WTMM

Рис. 13. Структурная блок-схема метода WTMM

На рис. 14 представлены характерные величины, вычисляемые с помощью методов мультифрактального анализа.

Характерные мультифрактальные величины теоретического спектра

Рис. 14. Характерные мультифрактальные величины теоретического спектра:

Но показатель, соответствующий наиболее вероятным флуктуациям во временном окне сигнала; Н2 степень корреляции; amin — наименьшие флуктуации в спектре; amax — наибольшие флуктуации в спектре; W — ширина муль-тифрактального спектра, вариабельность флуктуаций в спектре

Мультифрактальные показатели являются количественной мерой оценки самоподобия ВР и могут характеризовать функциональные изменения в регуляторных механизмах организма [61].

2. Методические рекомендации

В пакете MATLAB коэффициенты непрерывного вейвлет-преобразования рассчитываются с использованием функции cwt (у, a, w), где у — вектор значений анализируемого ВР, а — вектор масштабирующих параметров, w — название анализирующего базисного вейвлета (подробнее см. практическое задание 4).

Поскольку анализируется сигнал ВСР, то целесообразно использовать частоты, имеющие физиологическое значение, а именно диапазоны LF [0,04—0,15], VLF [0,003—0,04]. Диапазон IIFhc рассматривается из-за зашумленности [27; 61].

Алгоритм, представленный на блок-схеме полностью реализован в прил. 14.

Примечание. В MATLAB локальные экстремумы модулей коэффициентов ММ (а,Ь) находятся с использованием функции findpeaks'.

[pks, locs]=findpeaks (abs (sc (і,:)));

В MATLAB нахождение коэффициентов наклона производится с помощью команды polyfit'.

const=polyfit (log2 (A), log2 (Z(і,:)), 1);

tau (i)=const (1);

Имеет интерес посмотреть результат для различных базисных функций. Для метода WTMM можно использовать следующие вейвлеты [43]:

  • • вейвлеты Гаусса (gausN, где N є {1.. 8});
  • • Морле (morl);
  • • мексиканская шляпа (mexh);
  • • вейвлет Мейера (meyr);
  • • вейвлеты Добеши (dbN, где Ng{1..10});
  • • симлеты (symN, где Ne{2..8});
  • • койфлеты (coifN, где Ng{1..5}).

Порядок выполнения задания:

  • • запустите среду MATLAB и создайте новый скрипт;
  • • произведите импорт данных согласно подглаве 2.5;
  • • произведите очищение сигнала от артефактов в соответствии с прил. 4;
  • • реализуйте алгоритм получения мультифрактального показателя Херста по методу WTM М для двух частотных диапазонов LF и VLF согласно рекомендациям п. 2 настоящего задания;
  • • сохраните изображения мультифрактальных спектров, как указано в п. 3 прил. 2;
  • • результаты исследований внесите в отчет.

Практическое задание 13

Реализация метода мультифрактального детрендированного флуктуационного анализа

1. Краткие теоретические сведения

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

Различие методов WTMM и MFDFA в получении скейлин-гового показателя. Исходные сигналы с помощью линейной интерполяции предварительно проектируются на равномерную сетку. На рис. 15 приводится блок-схема для расчета мультифрактального показателя Херста аналогично методу WTMM.

Исходные сигналы с помощью линейной интерполяции предварительно проектируются на равномерную сетку. Полученные ВР избавляются от тренда методом скользящего окна, а значения интерполированного ряда/ (/), і = 1, 2,..., N, делятся по непересекающимся сегментам длины s, число которых равно целому значению Ns= [АД].

2. Методические рекомендации

Полная реализация вышеописанного алгоритма представлена в прил. 15.

Примечание. Длина сегмента s определяется схоже с масштабирующим параметром в методе WTMM:

  • 1
  • 5--------,

delta • f

где delta — шаг дискретизации,/ — исследуемая частота.

Исходный сигнал y(t)

Полиномиальный тренд

л(0 = ЁФ’-‘

____________________к=0_______________________

Момент дисперсии

F2(y,s) = у ZH(V - + ~ ^v(0|2

Флуктуационная функция (

z,«=k-zkM2}

Угловой коэффициент log2 Zq(s) = h(q) ? log2s + const

Скейлинговый показатель т(^) = q ? h(q) -1

Мультифрактальный спектр а _ dt(q) =qa- x(q) dq

Мультифрактальный показатель Херста ^2 = a l9=2

Структурная блок-схема метода MFDFA

Рис. 15. Структурная блок-схема метода MFDFA

Изменение частот происходит в диапазонах, аналогичных методу WTMM — диапазоны LF[0,04—0,15], VLF[0,003—0,04].

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

C=polyfit (Index, Y (Index), m);

fit=polyval (C, Index);

Порядок выполнения задания:

  • • запустите среду MATLAB и создайте новый скрипт;
  • • произведите импорт данных согласно подглаве 2.5;
  • • произведите очищение сигнала от артефактов в соответствии с прил. 4;
  • • реализуйте алгоритм получения мультифрактального показателя Херста по методу MFDFA для двух частотных диапазонов LFh VLF согласно рекомендациям п. 2 настоящего задания;
  • • сохраните изображения мультифрактальных спектров, как указано в п. 3 прил. 2;
  • • результаты исследований внесите в отчет.

  • [1] NN50 — количество пар последовательных отсчетов временного ряда NN на интервале наблюдения, различающихся более чем на 50 мс;
 
Посмотреть оригинал
< Пред   СОДЕРЖАНИЕ   ОРИГИНАЛ   След >