Интерполяционный полином Ньютона

Имея табличные данные, построим интерполяционный полином степени п в виде, предложенном Ньютоном

Р_n(x)=А₀+А₁(х–x₀)+А₂(х–x₀)(х–x₁)+×××+А_n(х-x_о)(х-х₁)...(х-x_n-1). (10.4)

Равносильный вариант полинома можно записать при симметричной перенумерации узлов исходной таблицы узлов и значений функции 0 n, 1 n-1, 2 n-2,...

P_n(x)=B_n+B_n-1(x-x_n)+B_n-2(x-x_n)(x-x_n-1)+…+B₀(x-x_n)(x-x_n-1)(x-x₁). (10.5)

Коэффициенты полиномов (10.4) и (10.5) определяются из условий Лагранжа

Pn(х_i) = f_i , 0 £ i £ n. (10.6)

Полагаем х = x₀, тогда в формуле (10.4) все слагаемые, кроме А₀, обращаются в нуль, следовательно

а₀ = f₀ . (10.7)

Затем полагаем x=x₁, тогда по условию (10.6) имеем

f₀ +A₁(x₁ – x₀) = f ₁,

откуда находим коэффициент

, (10.8)

который называется разделенной разностью первого порядка.

Величина f₀₁ близка к первой производной функции f(x) при малом расстоянии между узлами x₀ и x₁.

При x=x₂ полином (10.4) принимает значение

P_n (x₂)=f₀ +f₀₁ (x₂-x₀)+A₂(x₂-x₀)(x₂-x₁),

из условия Лагранжа (10.6) находим

, (10.9)

где

.

Величина называется разделенной разностью второго порядка, которая при близком расположении x_0, x_1, x₂ будет пропорциональна второй производной функции f(x). Аналогичным образом при x = x₃ находим коэффициент полинома Ньютона

, (10.10)

где

, .

Для коэффициента A_k методом математической индукции запишем следующее выражение:

. (10.11)

Полученные результаты сведем в табл. 1.

Таблица 1

x	F(x)

Для построения интерполяционного полинома Ньютона используются только диагональные элементы приведенной таблицы, остальные элементы являются промежуточными данными. Поэтому в программе, реализующей вычисление коэффициента полинома, разделенные разности для экономии памяти разместим в массиве, где первоначально хранились значения функции f(x) в узлах. Этот массив будет частично обновляться при вычислении разделенных разностей очередного порядка. Так, при вычислении разностей первого порядка элемент остается неизвестным (коэффициент А₀ (10.7)), элемент f₁ заменяется на f₀₁ (коэффициент A₁ (10.9)), f₂ - на f₀₂ и т.д. При вычислении разделенных разностей второго порядка первые два элемента массива f_i, где размещены коэффициенты А₀ и А₁ полинома, оставляем неизменными, остальные элементы заменяем разделенными разностями.

Таким образом, после вычисления все коэффициенты полинома Ньютона будут размещены последовательно в массиве узловых значений функции f(x).

Заметим, что добавление новых узлов в табл. 1 не изменит уже вычисленных коэффициентов, таблица будет дополнена новыми строками и столбцами разделенных разностей.

После определения коэффициентов полинома Ньютона вычисление его значений при конкретных аргументах х наиболее экономично проводить по схеме Горнера, получаемой путем последовательного вынесения за скобки множителей (х-х_i) в формуле (10.4)

(10.12)

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

Интерполяция сплайнами

Полиномиальная интерполяция не всегда дает удовлетворительные результаты при аппроксимации зависимостей. Так, например, при представлении полиномами резонансных кривых колебательных систем большая погрешность возникает на концах ("крыльях") этих кривых. Несмотря на выполнение условий Лагранжа в узлах, интерполяционная функция может иметь значительное отклонение от аппроксимируемой кривой между узлами. При этом повышение степени интерполяционного полинома приводит не к уменьшению, а к увеличению погрешности. Возникает так называемое явление волнистости.

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

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

j ^(IV) (х) = 0 . (10.13)

Функцию j(х) будем использовать для аппроксимации зависимости f(x), заданной в узлах х_о, х₁..., х_n значениями f₀, f₁,..., f_n.

Если в качестве функции j(x) выбрать полином, то в соответствии с уравнением (10.13) степень полинома должна быть не выше третьей. Этот полином называют кубическим сплайном, который на интервале хÎ[x_i-1 ,х_i] записывают в виде

(10.14)

где а_i, b_i, с_i, и d_i, - коэффициенты сплайна, определяемые из дополнительных условий; i = 1,2,..., n - номер сплайна.

В отличие от полиномиальной интерполяции, когда вся аппроксимируемая зависимость описывается одним полиномом, при сплайновой интерполяции на каждом интервале [x_i-1, х_i] строится отдельный полином третьей степени (10.14) со своими коэффициентами.

Коэффициенты сплайнов определяются из условий сшивания соседних сплайнов в узловых точках:

1) равенство значений сплайнов j(x) и аппроксимируемой функции f(x) в узлах - условия Лагранжа

j_i(x_i-1)=f_i-1, j_i(x_i)=f_i; (10.15)

2) непрерывность первой и второй производных от сплайнов в узлах

(10.16)

(10.17)

Кроме перечисленных условий, необходимо задать условия на концах, т.е. в точках х₀ и х_n. В общем случае эти условия зависят от конкретной задачи. Довольно часто используются условия свободных концов сплайнов. Если линейка не закреплена в точках вне интервала [х₀, х_n], то там она описывается уравнением прямой, т.е. полиномом первой степени. Следовательно, исходя из условий (10.17) непрерывности вторых производных сплайнов на концах интервала, запишем соотношения

(10.18)

(10.19)

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

Получим алгоритм определения коэффициентов кубических сплайнов из условий (10.15)–(10.19). Условия (10.15) в узлах x_i-1 и х_i после подстановки i-го сплайна принимают вид

a_i=f _i-1, (10.20)

(10.21)

где

.

Продифференцируем дважды сплайн (10.14) по переменной х:

(10.22)

(10.23)

Из условий непрерывности производных (10.16) и (10.17) при переходе в точке х_i от i-го к (i + 1)-му сплайну, с учетом выражений (10.22) и (10.23), получим следующие соотношения:

(10.24)

(10.25)

И, наконец, из граничных условий (10.18) и (10.19) на основании выражения для второй производной (10.23) получим, что

c₁ =0, (10.26)

c_n+3d_nh_n=0 . (10.27)

Соотношения (10.20), (10.21), (10.24)-(10.27) представляют собой полную систему линейных алгебраических уравнений относительно коэффициентов сплайнов а_i, b_i, с_i и d_i. Но, прежде чем решать эту систему, выгодно преобразовать ее так, чтобы неизвестными была только одна группа коэффициентов с_i.

Из уравнения (10.25) коэффициенты d_i выразим через коэффициенты с_i:

d_i= (c_i+1-c_i)/(3h_i) (10.28)

Объединяя уравнения (10.20), (10.21) с соотношением (10.28),представим коэффициенты b_i также через коэффициенты с_i:

b _i=(f_i-f_i-1)/h_i-(c_i+1+2c_i)h_i /3 . (10.29)

После подстановки выражений (10.28) и (10.29) в соотношение (10.24) получим уравнение, в которое входят только неизвестные коэффициенты с_i. Для симметричности записи в полученном уравнении уменьшим значение индекса i на единицу

(10.30)

где 2 £ i £ п.

При i = n, учитывая условие свободного конца сплайна, в уравнении (10.30) следует положить

с_n+1 = 0 . (10.31)

Таким образом, п – 1 уравнение вида (10.30) вместе с условиями (10.26) и (10.31) образует систему линейных алгебраических уравнений для определения коэффициентов с_i. Коэффициенты d_i и b_i вычисляются после нахождения с_i по формулам (10.28) и (10.29), коэффициенты а_i равны значениям аппроксимируемой функции в узлах в соответствии с формулой (10.20).

В каждое из уравнений типа (10.30) входят только три неизвестных с последовательными значениями индексов c_i-1, с_i, с_i+1. Следовательно, матрица системы линейных алгебраических уравнений относительно с_i является трехдиагональной, т.е. имеет отличные от нуля элементы только на главной и двух примыкающих к ней диагоналях. Для решения систем с трехдиагональной матрицей наиболее эффективно применять так называемый метод прогонки, являющийся частным случаем метода исключения Гаусса.

Рассмотрим алгоритм метода прогонки применительно к решаемой задаче. Для сокращения записи уравнение (10.30) представим в виде

h_i-1 c_i-1 + s_ic_i + h_ic_i+1 = r_i , (10.32)

где

s_i = 2 (h_i-1 + h_i),

. (10.33)

Метод прогонки разделяется на два этапа - прямой и обратный ходы. В процессе прямого хода метода прогонки вычисляют значения коэффициентов линейной связи каждого предыдущего неизвестного с_i с последующим с_i+1.

Из уравнения (10.32) при i = 2, с учетом граничного условия (10.26), установим связь коэффициента c₂ с коэффициентом c₃

c₂ = k₂-I₂ c₃, (10.34)

где k₂ и I₂ - прогоночные коэффициенты,

k₂ = r₂/s₂, I₂ = h₂/s₂.

Затем, подставляя выражение (10.34) в уравнение (10.32) при i = 3, получим линейную связь коэффициента c₃ c коэффициентом с₄:

c₃ = k₃ – I₃ c_{4 .}

Поступая аналогичным образом для любых соседних коэффициентов с номерами i и i + 1, можно установить линейную связь в виде

c_i = k_i –I_i c_i+1 . (10.35)

В процессе выполнения прямого хода метода прогонки необходимо вычислить значения всех прогоночных коэффициентов k_i и I_i, для которых получим рекуррентные соотношения. Подставим формулу связи (i – 1)-го и i-го коэффициентов

c_i-1 = k_i-1 – I_i-1 c_i

в уравнение (10.32), в результате получим

Сравнение последнего соотношения с формулой (10.35) позволяет записать рекуррентные формулы для определения прогоночных коэффициентов k_i и I_i:

(10.36)

Учитывая граничное условие (10.26) и соотношение

c₁ = k₁ – I₁ c₂ ,

а также полагая c₂ ¹ 0, задаем начальные коэффициенты k₁ = 0 и I₁ = 0. Затем по формуле (10.36) вычислим все n пар прогоночных коэффициентов k_i и I_i.

На основании соотношения

c_n = k_n–I_nc_n+1

и граничного условия (10.31) получим, что

c_n = k_n. (10.37)

Далее последовательно применим формулу (10.35) при i = n – 1, n – 2,..., 2и вычислим значения искомых величин c_n-1, c_n-2,...,c₂. Эта процедура называется обратным ходом метода прогонки.

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

После определения всех коэффициентов с_i другие коэффициенты сплайнов вычисляются по формулам (10.20), (10.28) и (10.29), после чего аппроксимирующую функцию j(x) можно рассчитать с помощью соотношения (10.14) в любой точке х на интервале [х₀, х_n ].

Лекция 11. Статистическая обработка данных в медицинской диагностике.

Общие понятия

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

Кардиоинтервалография (КИГ) использует графическое представление данных об изменениях СР при помощи регистрации последовательных кардиоинтервалов. Они измеряются по длительности RR-интервалов электрокардиограммы (ЭКГ) в секундах или миллисекундах (1 c = 1000 мс). Методика КИГ представлена на рис.11.1.

Рис. 11.1(а)

Рис. 11.1(б)

Рис. 11.1(в)

Вначале, как это видно на рис. 11.1(а), при регистрации ЭКГ измеряются последовательные RR-интервалы (первые пять из рис. 11.1(а)). Методика построения ритмограммы и гистограммы RR-интервалов отображена на рис. 11.1(б) и 11.1(в):

(б) – ритмограмма (Тrr - величина длительности RR-интервалов, N - порядковый номер RR-интервала). Ритмограмма строится путем откладывания значения RR-интервалов по оси ординат последовательно, согласно их порядкового номера.

(в) - гистограмма (N - порядковый номер RR-интервала в отрезке числовой оси. RR-интервалы обозначены как RR1 ... RR5). Построение гистограммы осуществляется следующим образом. Отрезок оси длительностей RR-интервалов (обычно от 0,4 до 1,4 с) разбивают на короткие участки (обычно длительностью 0,05 с). Затем по кардиоинтервалограмме для каждого участка подсчитывается количество кардиоинтервалов, длительность которых принадлежит этому участку. Полученное число откладывают на графике в виде столбика. Совокупный график называют гистограммой. Подсчет количества кардиоинтервалов, попадающих в отрезки числовой оси, можно производить по-разному. На рис.11.1 (в) приведен вариант графического представления гистограммы при обработке последовательного ряда RR-интервалов, представленного на рис. 11.1 (а), (б).

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

Числовые характеристики КИГ

Числовые характеристики КИГ включают следующие показатели:

1) Мода (МО) - наиболее часто встречающееся значение длительности кардиоинтервалов в гистограмме. МО указывает на наиболее вероятный уровень функционирования синусового узла. При симметричном распределении (гистограмме) величина МО совпадает с математическим ожиданием М.

2) Амплитуда моды (Амо) - показатель, равный числу RR-интервалов, попавших в отрезок числовой оси, соответствующий моде. Лучше всего нормировать этот показатель, указывая долю RR-интервалов, соответствующих моде, от всей выборки. Как правило, Амо выражается в этих случаях в процентах. Принято считать, что данный показатель отражает стабилизационный эффект централизации управления СР и зависит в основном от симпатических влияний.

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

, (11.1)

где - значения величин максимального диапазона гистограммы,

- минимального.

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

Вторичные показатели, получаемые при гистографическом анализе, включают индекс вегетативного равновесия (ИВР), вегетативный показатель ритма (ВПР), показатель адекватности процессов регуляции (ПАПР), индекс напряжения регуляторных систем (ИН). Показатели ВПР и ИН были предложены Г.И. Сидоренко в 1973 г., в дальнейшем Р.М. Баевский несколько модифицировал формулу ИН.

Индекс вегетативного равновесия указывает на соотношение между активностью симпатического и парасимпатического отделов нервной системы и рассчитывается по формуле:

(11.2)

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

ПАПР = Амо / Мо (11.3)

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

(11.4)

Индекс напряжения регуляторных систем свидетельствует о степени централизации управления СР и вычисляется по формуле:

(11.5)

Увеличение симпатического тонуса приводит к приросту Амо и уменьшению Мо и и, следовательно, нарастанию ИН. Усиление парасимпатических влияний, наоборот, ведет к уменьшению Амо, увеличению величины Мо и, следовательно, ИН. У здоровых людей в состоянии физического и психического покоя значения ИН составляют 80-140.

Статистический анализ СР

При статистическом анализе рассматривают последовательные кардиоинтервалы как ряд величин, производят расчет их статистических характеристик, но без наглядного гистографического отображения.

К статистическим показателям относятся математическое ожидание М, среднеквадратичное отклонение &, коэффициент вариации V, коэффициент асимметрии Аs, эксцесс Еx.

Математическое ожидание - величина, обратная средней частоте пульса (HR) или частоте сердечных сокращений (ЧСС) за 1 мин:

(11.6)

При анализе последовательности RR интервального ряда М вычисляется по формуле:

(11.7)

где RR_i -значение i-го интервала ритмограммы, N - число интервалов.

Соответственно, HR = 60/М (уд/мин).

Величины HR и ЧСС совпадают, но разница между ними состоит в том, что первая величина определяется при исследовании пульсаций артерий, а вторая - при исследовании работы сердца.

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

Среднеквадратичное отклонение динамического ряда RR-интервалов - вычисляется по формуле:

(11.8)

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

Коэффициент вариации V представляет собой показатель, который получается при нормировании среднеквадратичного отклонения по частоте пульса.

Показатели асимметрии Аs и эксцесса Еx определяют степень отличия случайного процесса - последовательности RR-интервалов - от так называемого нормального.

Коэффициент асимметрии Аs вычисляется по следующей формуле:

(11.9)

Коэффициент эксцесса Ех вычисляется по формуле:

(11.10)

Указанные показатели вариабельности СР, как и многие другие, обладают определенным недостатками, связанными с тем, что они должны применяться в условиях стабильного или квазистабильного состояния. Формально они могут быть вычислены и для нестационарных условий, при наличии трендов и т.п. Однако при этом существенно увеличиваются трудности в интерпретации этих показателей. Cтатистический и гистографический анализ СР имеет существенные недостатки. Так, эти методики не позволяют получить информацию о волновой структуре ритмограммы и охарактеризовать их динамические свойства при переходных процессах. В значительной мере такую оценку позволяют провести методы исследования корреляционной ритмографии автокорреляционного и спектрального анализа ритмограмм.

Корреляционная ритмография (скаттерография, двумерная гистография) - это отображение динамики СР в виде точек на прямоугольной системе координат. Проекция каждой точки на ось ординат представляет собой длительность последнего RR-интервала (RR1), а проекция на ось абсцисс - предшествующего (RR).

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

Рис. 11.2

При равенстве соседних RR-интервалов соответствующие им точки лежат на прямой, идущей по биссектрисе угла, образованного абсциссой и ординатой. В случае увеличения последующего RR-интервала точка располагается над биссектрисой, а при уменьшении - под ней. Изменения СР приводят к разбросу точек. По корреляционной ритмограмме можно определить М (центр совокупности точек), х (расстояние между наиболее удаленными точками по оси абсцисс или ординат). Чем медленнее период колебаний длительностей RR-интервалов, тем более вытянут эллипс вдоль биссектрисы. Поэтому одним из способов описания корреляционного облака является расчет отношения величин продольной его оси (а) и поперечной (б): а/б. Чем больше выражена медленная периодика, тем больше величина отношения а/б. С учетом величин М, а/б можно вычислить индекс функционального состояния (ИФС) по формуле:

(11.11)

Величина показателя тем выше, чем меньше напряжение регуляторных механизмов.

Лекция 12. Методы распознавания образов в биомедицинских сигналах.

Теоретическое введение

Алгоритмы распознавания желудочкового комплекса ЭКГ (QRS-комплекса) решают следующие основные задачи: обнаружение комплекса (т.е. установление факта его наличия на анализируемом участке) и определение характерных точек комплекса (опорной точки, служащей для измерения RR-интервала, точек начала и конца комплекса, а также крайних точек и вершин его зубцов). Можно выделить несколько основных групп методов распознавания QRS-комплекса при оперативном анализе ритма сердца по электрокардиосигналу (ЭКС):

· простейшие пороговые методы;

· структурные методы;

· методы сравнения с образцами (корреляционные методы);

· методы на основе цифровой фильтрации.

Простейшие пороговые методы

Простейшие пороговые методы основываются на применении несложных логических правил по отношению к исходному ЭКС или к его первой производной, в качестве оценки которой обычно используется первая разность отсчетов сигнала. Факт обнаружения комплекса фиксируется при превышении сигналом (или модулем сигнала) некоторого порога. Такие методы отличаются относительной простотой, но обладают невысокой устойчивостью к помехам и к изменению ЭКС. Кроме того, для обеспечения надежной работы этих алгоритмов необходима подстройка порога обнаружения QRS-комплекса для каждого пациента. Из-за этих недостатков простейшие пороговые методы находят ограниченное применение.

Метод сравнения с образцом

При использовании метода сравнения с образцами предполагается вычисление в текущем режиме взаимной корреляционной функции между входным ЭКС и одним или несколькими образцами желудочковых комплексов. Эти образцы могут представлять собой либо усредненные модели различных видов ранее обнаруженных комплексов, либо заранее определенных типовых комплексов. Обнаружение желудочкового комплекса может осуществляться по превышению полученной функцией взаимной корреляции заданного порога, что должно свидетельствовать о высокой степени линейной зависимости анализируемого фрагмента ЭКС и соответствующего образца. Такой алгоритм может дать хорошее качество обнаружения QRS-комплекса даже в условиях значительных помех. Кроме того, одновременно с обнаружением комплексов при этом решается и задача классификации их форм.

Достоинством таких алгоритмов является их адаптируемость в ходе анализа к форме сигнала каждого конкретного пациента. Однако реализация корреляционных методов распознавания QRS-комплекса в системах оперативной обработки ЭКС связана с чрезвычайно высокими требованиями к производительности используемого процессора и может быть осуществлена с применением специализируемых быстродействующих вычислительных устройств. В связи с этим часто предлагаются упрощенные методы получения оценок взаимной корреляционной функции, хотя результаты анализа в таких случаях оказываются несколько ниже.