Спасибо. Итак, давайте вначале разберемся со временн`ой нормировкой . Вопрос заключается в том, сколько точек (или сколько RR-интервалов) должно быть во входном массиве? Оказывается, что не имеет значение их количество, а важно, какое они образуют время обследования .

В сумме это время должно составлять около 333 секунд . И вот почему. Нижняя граница спектра VLF составляет 0.003 Hz . Расчитаем время снятия ритмограммы, которое будет соответствовать этой частоте. Для этого поделим 1 sec/0.003 Hz = 333 sec. Следует отметить также, что при спектральном преобразовании 0.003 Hz будет соответствовать еще и разности частот между соседними точками спектра, т.е. определять частотный шаг.

Здесь же отметим физичекую сущность первого элемента спектра. Это постоянная составляющая (нулевая частота). Она несет информацию о средней ЧСС за обследуемый период. В расчетах спектральной мощности ее учитывать не надо.

Теперь расчитаем, какая должна быть минимальная размерность спектрального массива. В общем случае она зависит от количества сердечных сокращений за 333 сек, а также алгоритма применемого спектрального преобразования.

К примеру, предполагая, что частота сердечных сокращений не будет превышать 90 уд.\мин, получим, что за 333 секунды произойдет около 333/60*90=~500 сердечных сокращений. Таким образом, размерность исходного массива с RR-интервалами должна быть не менее 500 точек .

Будем использовать алгоритм быстрого преобразования Фурье (БПФ). Его особенностью является не только малое время преобразования, но и то, что размерность спектрального массива должна быть кратна степени двойки. Таким образом, ближайшая размерность в данном случае составит 512 точек. Однако при БПФ всегда формируется симметричный спектр , состоящий из двух одинаковых, но зеркально расположенных относительно нулевой частоты спектральных частей. Следовательно, размерность частотного массива должна быть в два раза больше выбранной. Т.е. не в данном случае - не менее 1024 точек. Невыполнение этого условия приведет к катастрофическому искажению результатов пребразования.

Теперь, учитывая, что частотный шаг мы ранее расчитали и он составил 0.003 Hz, можно легко определить, где в спектральном массиве будут находиться частотные области VLF , LF и HF . Так, для VLF-части спектра (0.003-0.4 Hz) будуть соответствовать точки от 2 (Внимание! В нашем случае элементы нумеруются от 1. Поэтому 1 точке будет соответствовать нулевая частота) до 13 -ой спектрального массива, для LF-части спектра (0.4-0.15 Hz) -от 14 до 50 точки, а HF-части (0.15-0.4 Hz) - от 15 до 133 точки. Имея расчитанный спектр и границы частоных зон (VLF, LF, HF), легко получить спектральные мощности, подсчитав сумму квадратов амплитуд спектральных составляющих в каждой частотной области (см. текст программного модуля spec_lib ниже).

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

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

Из этого вытекает, что в качестве поверочного средства можно использовать обычные гармонические колебания с заданной частотой и амплитудой (Или их смесь). То есть, оцифровав, например, гармонику с и заданной амплитудой и частотой 0.01 Hz, мы получим в спектре спектральную составляющую с частотой 0.01 Hz и мощностью, соответствующей амплитуде входного сигнала. Все остальные гармоники спектра будут равны нулю. Среднеквадратическая мощность (P) синусоидального сигнала (гамоники) связана с ее амплитудой (A) как P=A*A/2 .

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

Таким образом, используя три гармоники, размещенные в частотых областях VLF, LF и HF мы имеем возможность моделировать любые виды частотных соотношений для проверки адекватности выбранного матаппарата частотного преобразования.

Здесь находится программа-иммитатор (HRV_Imit.EXE) , которая реализует описанный выше принцип. Она позволяет сформировать оцифрованную смесь из трех любых частот лежащих в диапазоне (VLF..HF) с произвольными мощностями вывести ее в текстовый файл export.rg . Структура файла отвечает российскому стандарту формата записи ритмограмм, где в каждой строке записывается по одному RR-интервалу в ASCII-формате. Такой файл легко импортировать в любую верифицируемую программу.

Ниже приведен код программного модуля, реализующего обработку ритмограммы с использованием алгоритма БПФ. Этот модуль использован в в тестовой программе (checkFFT.EXE ) , обрабатывающей текстовый файл ритмограммы с описанной выше структурой. Программу checkFFT.EXE можно также использовать для обработки ритмограмм, экспортированных и из других программ.

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

unit spec_lib;

{модуль, що реалізує перетворення ритмограми на спектр і обрахунок спектральних параметрів }
{написаний мовою Turbo Pascal, (c) Borland Inc.}
{використаний в тестовій програмі CheckFFT.EXE }
{(c) Pallar Ltd. Co, 2002, +(380)-432-523-041 http://www.pallar.com.ua/ }

INTERFACE

uses dos;

const EndArray = 1024; {розмiрнiсть спектрального масиву}
SpecWid : longint = 333000; {ширина спектру в мiлiсекундах}

var ECGArray : array [1..EndArray] of integer; {масив із вхідним сигналом}

VLF : longint; {спектр. потужнiсть в дiапiзонi 0.003-0.04 Hz}
LF : longint; {спектр. потужнiсть в дiапiзонi 0.04-0.15 Hz}
HF : longint; {спектр. потужнiсть в дiапiзонi 0.15-0.4 Hz}
nuLF : longint; {спектр. потужнiсть в дiапiзонi 0.04-0.15 Hz Normal Units}
nuHF : longint; {спектр. потужнiсть в дiапiзонi 0.15-0.4 Hz Normal Units}
TotalSpecPower : longint; {загальна спектральна потужнiсть}
LF_HF : real; {..........}
FL : real; {частота з НЧ максимумом спектр. потужностi}
N333 : longint; {к-сть точок до iнтервалу 333 сек вiд початку запису}

type
Complex = record re, im : real end;
MyArr = Array [ 1..EndArray ] of Complex;

procedure FFT( A : pointer; Dimension,Expon : word; mode : shortint);
procedure GetRRSpectrum;
procedure ResizeRRArray;

IMPLEMENTATION

procedure FFT ( A : pointer; Dimension,Expon : word; mode : shortint);
{--------------- процедура швидкого перетворення ФУР'Є ---------------------}
{ A - показчик на масив типу:

type complex=record re, im : real end;
myarr = array [ 1..Dimension ] of complex;

розмiрностi Dimension и 2**Expon;
mode = 1 - пряме перетворювання
мode = -1 - зворотнє перетворювання }


function exp( expon : byte) : integer;
begin exp := 2 shl (expon-1); end;

label 1;

var s,i,j,k,l,L1,L2,HalfOfDimension : integer;
SwapWord1,SwapWord2,v1,v2,w1,w2,my_phi : real;
A_ : ^myarr;
begin

A_ := A;

HalfOfDimension := Dimension div 2; j:=1;

for i:=1 to Dimension-1 do begin

if i >= j then k:=HalfOfDimension
else begin
SwapWord1 := A_^[j].re;
A_^[j].re := A_^[i].re;
A_^[i].re := SwapWord1;

SwapWord1 := A_^[j].im;
A_^[j].im := A_^[i].im;
A_^[i].im := SwapWord1;
k := HalfOfDimension;
end; {if }

1: if (k<j) then begin
j := j-k;
k := k shr 1;
goto 1;
end else j := j+k;

end; {do}


for l :=1 to Expon do begin
L1:=exp(l);
L2:=L1 shr 1;
v1:=1;
v2:=0;
my_phi:=pi/L2;
w1:=cos(my_phi);
w2:=mode*sin(my_phi);

for j:=1 to L2 do
begin
i:=j;
repeat
s:=i+L2;
SwapWord1:=A_^[s].re*v1-v2*A_^[s].im;
SwapWord2:=A_^[s].re*v2+v1*A_^[s].im;
A_^[s].re:=A_^[i].re-SwapWord1;
A_^[s].im:=A_^[i].im-SwapWord2;
A_^[i].re:=A_^[i].re+SwapWord1;
A_^[i].im:=A_^[i].im+SwapWord2;
i:=i+L1;
until i >= Dimension;
SwapWord1:=v1*w1-w2*v2;
v2:=v1*w2+v2*w1;
v1:=SwapWord1;
end; {for}
end; {for}
end; { FFT }


procedure GetRRSpectrum;
var p : pointer;
SSS : ^myarr;
NumPoint, i, k, repeater : integer;
x, TimeStep, AmplMasht : real;
begin

Mark(p); New(SSS);

for NumPoint := 1 to EndArray do begin {перенесення сигналу з масиву ECGArray в масив SSS}
SSS^[ NumPoint ].re := ECGArray[NumPoint]; {дійсна частина}
SSS^[ NumPoint ].im := 0; {неявна частина}
end; {for}

FFT(@SSS^, 1024, 10, 1); {перетворення Фур`є, одержання спектру сигналу}

{2 - коефiцiент втрати потужностi через симетричнiсть спектру }

{обрахунок спектральної потужності, що пропорційна сумі квадратів дійсної і неявної частин спектру}
AmplMasht := sqrt(2);
for NumPoint := 1 to EndArray do
SSS^[ NumPoint ].re :=
(sqr(SSS^[ NumPoint ].re/EndArray)+sqr(SSS^[ NumPoint ].im/EndArray))/AmplMasht;

{максимум спектру}
AmplMasht := SSS^[ 1 ].re;
for NumPoint := 2 to EndArray do if SSS^[ NumPoint ].re > AmplMasht then

AmplMasht := SSS^[ NumPoint ].re;
{нормування до 30000}
AmplMasht := AmplMasht + 0.0001;
AmplMasht := 30000/AmplMasht;
for NumPoint := 1 to EndArray do
SSS^[ NumPoint ].re := SSS^[ NumPoint ].re * AmplMasht;

for NumPoint := 1 to EndArray do
ECGArray[NumPoint] := round(SSS^[ NumPoint ].re);

Release(p);

repeater := 333 div SpecWid; if repeater = 0 then repeater := 1;

ECGArray[1] := 0; {обнулимо нульову частоту}

{приводимо масштаб спектру до амплiтуди вхiдного сигналу.
через те, що спектр симетричний, то i масштабний коеф. визначається
як 1024/2,
експозицiя складає 333 сек. =>
верхня частота складає : - 1.5375 Hz
крок по частотi - 0.003 Hz }

TimeStep := 3.0/AmplMasht; {3.0 ms per point}
x := 0;
{0.003-0.04 Hz}
for NumPoint := 2 to 13 do x := x + abs(ECGArray[NumPoint]);
VLF := trunc( x * TimeStep ); {потужність VHF-частини спектру}

{0.04-0.15 Hz}
x := 0;
for NumPoint := 14 to 50 do x := x + abs(ECGArray[NumPoint]);
LF := trunc( x * TimeStep ); {потужність LF-частини спектру}

FL := 15; x := abs(ECGArray[NumPoint]);
for NumPoint := 14 to 50 do
if x < abs(ECGArray[NumPoint]) then begin
x := abs(ECGArray[NumPoint]);
FL := NumPoint;
end;
FL := (FL-1)*0.003; {частота з НЧ максимумом спектр. потужностi }

{0.15-0.4 Hz}
x := 0;
for NumPoint := 51 to 133 do x := x + abs(ECGArray[NumPoint]); {!!! 51..133}
HF := trunc( x * TimeStep ); {потужність HF-частини спектру}

TotalSpecPower := VLF + LF + HF; {Загальна потужність}

if TotalSpecPower <> 0 then begin
nuLF := trunc(LF/(TotalSpecPower-VLF)*100); {нормалізовані потужності}
nuHF := trunc(HF/(TotalSpecPower-VLF)*100);
end else begin
nuLF := 0;
nuHF := 0;
end;

if HF > 0 then LF_HF := LF/HF {співвідношення LF/HF}
else LF_HF := 0;

end; {getrrspectrum}

procedure ResizeRRArray;
{процедура часового нормування (розтягування\стискування)первинного масиву до розмiрностi 1024}
{приведемо RR-масив до 333 точок, має вiдповiдати 5-хвилиннiй експозицiї}
var HorScaleStep : real;
i, CurrPoint : integer;
TmpArray : array [1..EndArray] of integer;
begin
HorScaleStep := EndArray / N333;

for i := 1 to EndArray do TmpArray[ i] := 0;

TmpArray[ 1] := ECGArray[1];
for i := 2 to N333 do begin
CurrPoint := trunc(i*HorScaleStep);
if CurrPoint <= EndArray then TmpArray[CurrPoint] := ECGArray[i];
end;

{заповнимо проущенi елементи масиву ...}
i := 0;
repeat
inc(i);
until ( i >= EndArray ) or (TmpArray[i] <> 0);
TmpArray[1] := TmpArray[i];

for i := 2 to EndArray do
if TmpArray[i] = 0 then TmpArray[i] := TmpArray[i-1];

{приберемо постiйну складову}
CurrPoint := TmpArray[1];
for i := 2 to EndArray do if CurrPoint > TmpArray[i] then CurrPoint := TmpArray[i];
for i := 1 to EndArray do begin
ECGArray[i] := TmpArray[i] - CurrPoint;
end;

end; {resizeRRarray}

END.