В сумі цей час має складати біля 333 секунд. І ось чому. Нижня границя спектру VLF складає 0.003 Hz. Розрахуємо час зняття ритмограми, який буде відповідати цій частоті. Для цього поділимо 1sec/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.04 Hz) їм будуть відповідати точки від  2 (Увага! В нашому випадку елементи нумеруються від 1. Тому 1-ій точці буде відповідати нульова частота) до 13-ой спектрального масиву, для LF-частини спектру (0.04-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.