CC BY
1710. Черчиньяни К. Математические методы в кинетической теории газов. М.: Мир, 1973. - 247 с.
11. Карлеман Т. Математические задачи кинетической теории газов. - М.: ИЛ, 1960. - 124 с.
12. Фирсов А.Н. Исследование решений уравнения Больцмана, близких к равновесным : дисс. ... канд. физ.-мат. наук : 01.02.05 : защищена 18.12.1975 : утв. 21.04.1976 / Фирсов Андрей Николаевич. - Л.: ЛГУ, 1975. - 129 с.
13. Maslova N.B. Nonlinear evolution equations. Kinetic approach. // Series on Advances in Mathematics for Applied Sciences. . - Vol. 10. - Singapore: World Scientific Publishing Co., 1993. - 193 p.
14. Маслова Н. Б. Математические методы исследования уравнения Больцмана // Алгебра и анализ. . - Т. 3, вып. 1. - Л., 1991. - С. 3-56.
15. Вычислительные методы в динамике разреженных газов. Сборник статей под ред. В.П. Шидловского. - М.: Мир, 1969. - 277 c.
16. Muntz ER. Rarefied gas dynamics // Ann. Rev. Fluid Mech. - 1989. - Vol. 21. -Pp. 387-417.
17. Деч Г. Руководство к практическому применению преобразования Лапласа и Z-преобразования. - М.: Наука, 1971. - 288 с.
УДК 303.732.4
doi:10.18720/SPBPU/2/id21-152
Малыхина Галина Фёдоровна 1,
профессор, д-р техн. наук, профессор; Семенютин Владимир Борисович 2, профессор, д-р техн. наук; Мухидинова Наргис Шухратовна 3,
магистр, аспирант
ФУРЬЕ И ВЕЙВЛЕТ АНАЛИЗ В МЕДИЦИНСКОЙ СИСТЕМЕ КОНТРОЛЯ АУТОРЕГУЛЯЦИИ МОЗГОВОГО
КРОВООБРАЩЕНИЯ
1 3
' Россия, Санкт-Петербург, Санкт-Петербургский политехнический университет Петра Великого, 1 g_f_malychina@mail.ru; Россия, Санкт-Петербург, Национальный медицинский исследовательский центр им. В. А. Алмазова
Аннотация. В исследовании выполнен анализ результатов интеллектуальной обработки сигналов системного артериального давления и скорости кровотока и в средней мозговой артерии головного мозга в диапазоне волн Майера. Анализ использует вычисление кросс-когерентности и угла сдвига фаз между сигналами. Для вычисления показателей кросс-когерентности и фазового сдвига разработаны алгоритмы, основанные на дискретном преобразовании Фурье и на непрерывном вейвлет-преобразовании дискретного сигнала с использованием вейвлет-функций Морле. Выполнено сопоставление этих двух подходов.
Ключевые слова: анализ Фурье, вейвлет-анализ, мозговое кровообращение, волны Майера.
Galina F. Malykhina\
Doctor of Technchnic Sciences, Professor;
Vladimir B. Semenyutin2, Doctor of Technchnic Sciences, Professor;
Nargis Sh. Mukhidinova , Master of Science, Postgraduate Student
FOURIER AND WAVELET ANALYSIS IN THE MEDICAL SYSTEM FOR MONITORING THE AUTOREGULATION OF CEREBRAL CIRCULATION
1 3 Peter the Great St. Petersburg Polytechnic University, St. Petersburg, Russia, 1 g_f_malychina@mail.ru;
Almazov National Medical Research Centre, Ministry of Health of Russia,
St. Petersburg, Russia
Abstract. The study analyzed the results of intelligent processing of signals of systemic arterial pressure and blood flow velocity in the middle cerebral artery of the brain in the Mayer wavelength range. The analysis uses the calculation of cross-coherence and phase angle between signals. To calculate the indicators of cross-coherence and phase shift, algorithms have been developed based on the discrete Fourier transform and on the continuous wavelet transform of a discrete signal using the Morlet wavelet functions. Comparison of these two approaches is made.
Keywords: Fourier analysis, wavelet analysis, cerebral circulation, Mayer waves.
Введение
Современным направлением развития медицинских систем является создание единой технологической платформы для сбора, накопления данных и применения к ним новых цифровых методов анализа, машинного обучения и искусственного интеллекта. Технология анализа ауторе-гуляции мозгового кровотока является современным методом исследования состояния пациентов при различных поражениях головного мозга: при черепно-мозговой травме, при ишемическом инсульте, внутричерепном кровоизлиянии, при церебральной артериовенозной мальформации. Создание цифровой системы для измерения и анализа сигналов, выполняющей диагностирование и прогнозирование состояния пациента в реальном времени, является актуальной задачей цифровой медицины.
Система церебральной ауторегуляции проявляется согласованностью фазового сдвига между колебаниями CАД и скоростью кровотока (ЛСК) в средней мозговой артерии. Ауторегуляция нарушается при некоторых неврологических заболеваниях [1, 2]. Одним из наиболее информативных показателей является сдвиг фаз между колебаниями CАД и ЛСК в М-диапазоне волн 0.080-0.120 Гц. Сложность состоит в том, что
полезный сигнал, характеризующий ауторегуляцию, наблюдается в смеси с другими физиологическими сигналами, поэтому получение диагностической информации связано с обработкой сигналов.
Для анализа выполнены эксперименты и получены результаты измерения АД и ЛСК в условиях принудительного дыхания с частотой, соответствующей средней частоте М-волн (0,10 Гц), и в условиях спонтанного дыхания, оказывающего внешнее воздействие на систему мозгового кровообращения. Таким образом, время-частотный кросс-анализ АД и ЛСК в диапазоне волн Майера в режиме реального времени при форсированном и спонтанном дыхании является адекватным неинвазивным методом оценки состояния системы мозгового кровообращения. Этот метод может использоваться для диагностики и прогнозирования нарушения церебральной ауторегуляции.
Целью исследования является анализ интеллектуальной системы диагностирования и прогнозирования нарушений мозгового кровотока, построенной на основе Фурье и вейвлет преобразования сигналов АД и ЛСК в базальных артериях головного мозга.
1. Кратковременное преобразование Фурье
Для выполнения исследования были получены результаты измерения системного артериального давления (САД) и линейной скорости кровотока (ЛСК) в средней мозговой артерии для двух режимов: для произвольного и для принудительного дыхания. Для дискретных измерительных сигналов, полученных с периодом дискретизации 0,01 с, был произведен анализ Фурье и непрерывный вейвлет анализ. По результатам анализа вычислены значения когерентности и фазового сдвига по Фурье и вейвлет-когерентности и фазового сдвига в пространствах вейвлет-разложения.
Традиционный анализ выполняют в частотной области на основе преобразования Фурье, при этом предполагается, что сигналы локально стационарны [3, 4]. Для нестационарных сигналов обычно применяют оконное преобразование Фурье, которое в случае гауссова окна называется преобразованием Габора [5]. Размер окна выбран таким образом, чтобы получить достаточную точность разложения сигнала в области М-волн, в частности, выбрано окно длительностью 20,48 секунды, содержащее И^п = 2048 отсчетов.
Анализируемые сигналы САД и(п) и скорости кровотока ЛСК р(п) присутствуют в смеси с другими физиологическими сигналами и с шумами, поэтому при вычислении их характеристик требуется выполнять сглаживание. Для этой цели было предложено из нескольких окон составить фрейм, и определить средние для всех окон значения спектральных характеристик сигналов внутри данного фрейма. Величина смещения зависит от вида окна. Прямоугольное окно рекомендуется смещать на %
длины, а окно Хемминга и окно Ханна — на % длины. Смещение фрейма выбрано равным смещению окна, которое составляет 512 отсчета.
При цифровой обработке сигналов САД и(п) и ЛСК v{ri) использованы фреймы длиной Nframe, состоящие их окон длиной Nwin =2048, сме-
1
щающиеся на величину Nshift =-Nwin = 512 с перекрытием
з
Nover=- Nwin. Количество окон в одном фрейме зависит от требуемой величины сглаживания L= frame I .
/ "shift
Для уменьшения искажений частотных спектров последовательность значений сигналов САД и ЛСК умножают на окно Ханна
й(п) = 0.5 (l -cos(2n—-ri—-JJu(ri), ri=1...Nwin;
\ \ "win '
v(n) = 0.5 il — cos [2n-и v(n),
\ V NWin — 1//
центрированных последовательностей: и(п) = u(ri) — mean(u), v(ri) = vt — mean(v), вычислены взаимные спектры путем Nwin - точечного преобразования Фурье комплексной последовательности суммы и(п) + jv{ri) сигналов — Suv(k) и собственные спектральные плотности сигналов и(п) и v(n):
W*) = U(k)U*(k)/Nwin ; SViV(k) = V(k)V*(k)/Nwin. Вычисления повторяем для каждого смещенного положения окна, в результате получаем K значений взаимных и собственных спектральных плотностей.
Функцию когерентности рассчитываем по сглаженным спектральным плотностям для каждого фрейма:
п= 1 N ■ ■
IL ± ... iVWin ;
Yu,v{k) —
(Ю
Средний фазовый сдвиг &иу между сигналами получаем исходя из соотношения к) = ехр(-)&иу).
После вычисления функции когерентности и фазового сдвига сигналов для всех окон текущего фрейма выполнен сдвиг фрейма.
Функцию когерентности и фазовый сдвиг для каждого следующего фрейма получаем исходя из соответствующих значений предыдущего фрейма путем удаления характеристик первого окна и добавления характеристик последнего окна нового фрейма. Тогда для каждого т-го фрейма получаем функцию когерентности уи1У(к,т) и функцию фазового сдвига 0и1р(к, т.), где к — номер отсчета по частоте, т. — номер отсчета
по времени. Дискретность отсчета времени равна Nshi^t = win/ 2. Поэтому результатом алгоритма предварительной обработки является двумерная время-частотная функция когерентности сигналов Yu,v(k>т) и двумерная время-частотная функция фазового сдвига &UiV(k, т).
Недостатком преобразования Фурье с окном является то, что размер окна не подстраивается под спектральные свойства сигнала [6, 7]. Во многих случаях теория вейвлетов может заменить традиционный подход, основанный на оконном преобразовании Фурье [8]. В теории вейвлетов используется адаптивное окно, которое меняет свой размер в зависимости от частоты исследуемого сигнала. Вейвлет-преобразование сигнала [5, 6], позволяют обнаруживать локализованные прерывистые периодичности, связанные с определенными нарушениями ауторегуляции головного мозга.
Сопоставление непрерывного и дискретного вейвлет-преобразования показало, что для обнаружения гармонических составляющих сигналов САД и ЛСК лучше подходит непрерывное вейвлет-преобразование, использующее гармонический вейвлет-базис.
2. Непрерывное вейвлет-преобразование
Непрерывное вейвлет-преобразование (НВП) дискретных сигналов САД u(t) и ЛСК v(t) определяется как дискретная свертка с масштабированным и нормализованным вейвлетом [9,10]:
ZN~X f(n-n')6t\ л v^1 f(n-n')6t\ -— j, V(s, n)= > -— j,
nr=0 ^ ' nr=0 ^ '
t! I * f(n—n')St\ где n= у— параметр дискретного времени, у (—--) — нормализованная комплексно-сопряженная вейвлет-функция, смещенная на величину n'St и масштабированная на величину s.
По теореме о свертке вейвлет-коэффициенты могут быть более эффективно вычислены в виде обратного ДПФ от произведения Фурье-образов сигнала и вейвлет-коэффициентов.
В качестве базиса рассмотрены неортогональные вейвлеты, которые эффективны для анализа временных рядов и апериодических сдвигов: комплексный вейвлет Морле, комплексный вейвлет Поула и вещественный вейвлет, представляющий разность Гауссианов (DOG) [11, 12]. Эти вейвлет-базисы хорошо локализованы как по частоте, так и по времени. Предварительные расчеты показали, что все вейвлет-базисы дают неплохие результаты, но несколько лучшие результаты получаются с использованием вейвлета Морле, графики для которого показаны на рисунках 1.a и 1.b.
a)
b)
Рис. 1. Вейвлет Морле: a) во временной области; b) в частотной области
Алгоритм для вычисления когерентности и фазового сдвига сигналов посредством НВП применяется к двум дискретным центрированным
сигналам ип = %п- и Vn = (n- Zn.
Для работы в реальном времени вычисление выполняется для последовательных фреймов длиной: N = 2/¿oor(iog2(N0)+i, где N' — первоначальная длина фрейма анализируемых сигналов.
НВП представляет собой свертку сигнала и вейвлет- функции, которую лучше вычислять в частотной области. Параметры вейвлет-разложения, такие как минимальный масштаб s0 , максимальный уровень разложения J, конкретные уровни разложения j = 0,1 ... / и вектор масштабов Sj определены по формулам:
Sq = 2St,
NuSt
J = Sj~1 log 2^——), Sj = s0 2 W.
При выборе безразмерного параметра частоты (ш0) следует учесть, что при увеличении ш0 увеличивается частотная локализация, но уменьшается временная локализация сигнала. Для анализа сигналов был выбран параметр вейвлета Морле ш0 = 6, приблизительно соответствующий частоте дыхания.
ДПФ аналитического вейвлета Морле определено по формуле:
ijj(sMk) = п 2
Это выражение содержит функцию Хэвисайда:
1, ык > 0 ык< 0
На каждой шкале коэффициентов по формуле:
выполняется нормализация вейвлет-
\
2ns ~8t
Поэлементное умножение Фурье-образа анализируемого сигнала и Фурье-образов вейвлетов на каждом уровне разложения и последующее обратное ДПФ — преобразование, позволяет получить коэффициенты вейвлет-разложения сигналов ип и vn в следующем виде: си{п, s) = Е^й^^е^8*, cv(n, s) = Z%=¡Vki¡¡\sú>k)ei*>*nSt■
Выполнение вычисления вейвлет-коэффициентов области Фурье дало возможность уменьшить объем вычислений. Поскольку сигналы и(п) и v(n) присутствуют в смеси с другими физиологическими сигналами и с шумами, при вычислении вейвлет-преобразования сигналов необходимо выполнять сглаживание [13-15].
Сглаживание по времени более эффективно выполнить в частотной области с использованием окна Гаусса. Сглаживание по масштабам выполнено путем осреднения с прямоугольным окном. В результате анализа двух сигналов и(п) и v(n) получены вейвлет-коэффициенты си(п, s) и cv(n, s), сглаженные по времени и масштабам.
Вейвлет-кросс-спектр характеризует общую энергию двух сигналов, которая отлична от нуля, если два сигнала коррелируют друг с другом, и исчезает, если два сигнала независимы:
Cu,v(n>= £u(n, s)cv(n, s).
В общем случае вейвлет-кросс-спектр является комплексной функцией, которая представлена в виде амплитуды и фазы:
/ Im(cu>v(n, 5)) сиАп>= [£иАп> s)\exp\arctg-73-ту Re[cuv(n, s))
Когерентность или согласованность двух сигналов может быть определена как модуль нормализованного кросс-спектра. Когерентность определяет линейную связь между двумя сигналами. Величина когерентности варьируется от нуля до единицы. Квадрат нормированного значения когерентности определяется по формуле:
[си(п, s)cv(n,
Hu,v(n>s) —
12"
[c£(n, s)J \cv(n, s)J' При наличии когерентности сигналов, заданной условием Hu,vin's) > 0.6 локальный фазовый сдвиг сигналов может быть получен по следующей формуле:
0UlV(n, s) = arctg——ï—--—.
нe{cUtV{n, s))
Значение фазового сдвига позволяет судить о степени синхронизации сигналов при наличии временной задержки между ними.
3. Результаты эксперимента
Сигнал разделен на фреймы, размер фрейма составляет 2048*8=16384.
Для проверки работы программы были вычислены значения вейвлет-когерентности и угла сдвига в пространстве разложения, соответствующем 0.1 Гц, для фрейма длиной 2048*8=16384 отсчетов. При проведении эксперимента на этом фрейме было воспроизведено воздействие на волонтера, которое должно было повлиять на ауторегуляцию. На рисунке 2 показан график зависимости от времени САД и ЛСК в центральных мозговых артериях в диапазоне М-волн на этом участке.
На рисунке 3 отражена зависимость от частоты квадрата вейвлет-когерентности и фазового сдвига между АД и ЛСК в центральных мозговых артериях в диапазоне М-вол (от 0.08 до 0.12 Гц).
Максимальная когерентность сигналов в диапазоне М-волн достигается на частоте 0.10 Гц и фазовый угол межу сигналами на этой частоте достигает от 2 до 3 радиан.
Высокая когерентность сигналов на частоте 1.0-1.6 Гц обусловлена частотой сердечных сокращений и фазовый угол межу сигналами на этой диапазоне частот составляет 1 радиан.
500 450
cz
< 400
О
с§ 350
Q.
о 300
с;
S 250
g 200 150 100
0 20 40 60 80 100 120 140 160 180 Время (сек)
Рис. 2. Зависимость от времени САД и ЛСК в центральных мозговых артериях в диапазоне М-волн (от 0.08 до 0.12 Гц)
Вейвлет-когерентность 2 (Время=0.01*Фрейм/2)
0.1 0.2 0.3 0.4 0.5 0.С а) Р (ГЦ)
' Фаза (Время=0.01*Фрейм/2)
с[2 аз
б)
Р (Гц)
Рис. 3. Зависимость от частоты квадрата вейвлет-когерентности и фазового сдвига между АД и ЛСК в центральных мозговых артериях в диапазоне М-волн
(от 0.08 до 0.12 Гц)
Были предложены два алгоритма: алгоритм, основанный на преобразовании Фурье с окном (КПФ), и алгоритм, основанный на НВП. Различие алгоритмов состоит в том, что алгоритм КПФ использует окно постоянного размера, что и поэтому имеет одинаковое разрешение во всем диапазоне частот. Алгоритм НВП использует базис, размер которого изменяется, имеет разное значение в разных диапазонах частот. Это свойство позволяет получить более высокое разрешение по времени и по частоте в диапазоне М-волн.
Для подавления шумов алгоритм КПФ использует осреднение по времени результатов вычисления, полученных по нескольким соседним фреймам. В то время как алгоритм НВП для подавления шумов использует низкочастотные фильтры по двум переменным, по времени и по частоте, поэтому шумы подавляются более полно. Для ускорения вычислений алгоритм НВП применяет вычисление сверток в частотной области.
Для сравнения алгоритмов КПФ и НВП выполнено вычисление когерентности и фазового сдвига сигналов АД и ЛК с использованием одинаковых данных при одинаковой длине фрейма. Результаты вычисления когерентности показаны на рис. 4. При вычислении по алгоритму КПФ использован размер фрейма 2048*8 отсчетов, перекрытие фреймов 2048 отсчетов, число отсчетов для вычисления преобразования Фурье 2048*2. Результат показан на рис. 4а. При вычислении по алгоритму НВП использован размер фрейма 2048*8 отсчетов, перекрытие фреймов отсутствует, число отсчетов для вычисления преобразования Фурье 2048*2. Результат показан на рис. 4б.
3
0
Рис. 4. Зависимости когерентности от частоты, полученные по алгоритму КПФ (а) и
алгоритму НВП (б)
Алгоритм НПВ позволяет получить более гладкие оценки характеристик сигналов.
Благодарности
Работа выполнена при финансовой поддержке гранта РФФИ (договор 19-29-01190\19 по теме «Разработка и применение методов математического анализа физиологических параметров для оперативной диагностики патологических состояний пациентов с оценкой рисков их перехода в критические состояния в режиме реального времени»).
Список литературы
1. Semenyutin V.B., Asaturyan G.A., Nikiforova A.A., Aliev V.A., Panuntsev G.K., Iblyaminov V.B., Savello A. V., Patzak A. Predictive value of dynamic cerebral autoregulation assessment in surgical management of patients with high-grade carotid artery stenosis. // Front. Physiol. - 2017. - Vol. 8. - URL: https://doi.org/10.3389/fphys.2017.00872.
2. Aliev V., Semenyutin V., Panuntsev G. Predictive Value of Cerebrovascular Reserve in Patients with Carotid Artery Stenosis for Choosing Treatment Strategy. // Int. J. Pathol. Clin. Res. - 2019. - Vol. 5. - URL: https://doi.org/10.23937/2469-5807/1510086.
3. Sejdic E.; Djurovic I.; Jiang J. Time-frequency feature representation using energy concentration: An overview of recent advances. // Digital Signal Processing. - 2009. -Vol. 19 (1). - Pp. 153-183. doi:10.1016/j.dsp.2007.12.004.
4. Jacobsen E., Lyons R. The sliding DFT. // Signal Processing Magazine. - Vol. 20, issue 2. - Pp. 74-80 (March 2003).
5. Jont B. Allen (June 1977). Short Time Spectral Analysis, Synthesis, and Modification by Discrete Fourier Transform. // IEEE Transactions on Acoustics, Speech, and Signal Processing. ASSP-25 (3). - Pp. 235-238. doi:10.1109/TASSP.1977.1162950.
6. Malykhina G.F., Merkusheva A.V. Classes of transformations of a non-stationary signal in information-measuring systems. VI. Correspondence of the form of covariance and the type of time-frequency transformation. // Sci. instrumentation. - 2007. - 17. -Pp. 75-87.
7. Merkusheva A.V., Malykhina G.F. The generalized Fourier transform method for time-frequency transformations, multiplexing and filtering of non-stationary signals in information systems. // Sci. Instrum. - 2006. - 16. - Pp. 85-96.
8. Malykhina G.F., Merkusheva A. V. Metrological aspects of non-stationary signal transformation for dynamic spectrum analysis // 10th IMEKO TC7 Symposium on Advances of Measurement Science. - 2004.
9. Stefanovska A., Bracic M., Kvernmo H.D. Wavelet analysis of oscillations in the peripheral blood circulation measured by laser Doppler technique // IEEE Trans. Biomed. -1999. - Eng. 46. https://doi.org/10.1109/10.790500.
10. Malykhina G., Salnikov V., Semenyutin V., Tarkhov D. Digitalization of medical services for detecting violations of cerebrovascular regulation based on a neural network signal analysis algorithm. // ACM International Conference Proceeding Series. - 2020. https://doi.org/10.1145/3444465.3444526.
11. Bassani T., Cesar Nievola J. Brain-computer interface usingwavelet transformation and naïve bayes classifier. // Advances in Experimental Medicine and Biology. - 2010. https://doi.org/10.1007/978-0-387-79100-5_8.
12. Birch A.A., Dirnhuber M.J.F.R.C.A., Hartley-Davies R., Iannotti F., Neil-Dwyer G. Assessment of autoregulation by means of periodic changes in blood pressure. // Stroke. - 1995. - Vol. 26. https://doi.org/10.1161/01.STR.26.5.834.
13. Bracic M., Stefanovska A. Wavelet-based analysis of human blood-flow dynamics. // Bull. Math. Biol. - 1998. - Vol. 60. https://doi.org/10.1006/bulm.1998.0047.
14. Mallat S. A Wavelet Tour of Signal Processing, 2009. https://doi.org/10.1016/B978-0-12-374370-1.X0001-8.
15. Morlet J.: Sampling theory and wave propagation. // Issues Acoust. signal/image Process. Recognit. - 1983. https://doi.org/10.1007/978-3-642-82002-1_12.
УДК 338.246.2
аоЫ0.18720/8РВРШМ21-153
Шаталова Ольга Михайловна1,
ст. науч. сотрудник, д-р экон. наук, доцент; Касаткина Екатерина Васильевна2,
доцент, канд. физ.-мат. наук, доцент
КЛАССИФИКАЦИЯ ПРОМЫШЛЕННЫХ РЕГИОНОВ РФ В ИССЛЕДОВАНИИ УСЛОВИЙ УСТОЙЧИВОГО РАЗВИТИЯ: АДАПТИРОВАННЫЙ АЛГОРИТМ Л-MEANS
1 Россия, Ижевск, Удмуртский филиал Института экономики УрО РАН,
oshatalova@mail.ru; 2 Россия, Ижевск, Ижевский государственный технический университет имени М.Т. Калашникова, e.v.trushkova@gmail.com
Аннотация. Статья представляет разработанный алгоритм кластерного анализа, применимый в целях классификации промышленных регионов по признаку их экономической специализации. Предложенный алгоритм основан на использовании метода ^-средних, адаптированного с учетом теоретического положения концепции Фишера-Кларка, условий организации государственной системы статистических из-
