Diagnostika pacientů s Parkinsonovou chorobou

Testy 3 - shrnutí výsledků

V průběhu minulých (1, 2) experimentů se skripty a jednotlivé funkce postupně vyvíjely. Nakonec jsem dospěl k této verzi funkcí:

function [v, vx, vy, a, ax, ay] = values(X, Y)
  x = sgolayfilt(X, 3, 5);
  y = sgolayfilt(Y, 3, 5);
  dt = 5;

  dx = filter([.5 0 -.5], 1, x);
  dy = filter([.5 0 -.5], 1, y);
  dx(1) = 0;
  dy(1) = 0;
  dx(2) = (x(2) - x(1));
  dy(2) = (y(2) - y(1));

  ds = sqrt((dx .* dx) + (dy .* dy));
  v = ds ./ dt;
  vx = dx ./ dt;
  vy = dy ./ dt;

  dvx = filter([1 -2 1], 1, x);
  dvy = filter([1 -2 1], 1, y);
  dvx(1) = 0;
  dvy(1) = 0;
  dvx(2) = dvx(3) / 2;
  dvy(2) = dvy(3) / 2;

  dvs = sqrt((dvx .* dvx) + (dvy .* dvy));
  a = dvs ./ dt;
  ax = dvx ./ dt ./dt;
  ay = dvy ./ dt ./dt;
return

Funkce values vypočítá ze souřadnic v polích X, Y rychlosti a zrychlení ve složce x i y. Nejdříve lehce přefiltruje SG filtrem nízkým řádem a krátkým oknem pro odstranění kvantizační chyby a šumu. Protože implementace filtru je s fázovým posuvem, je nutné provést korekci v prvních dvou prvcích vypočtených polí.

function [len, l] = getLength(koef, length)
  l = koef*ones(length,1);

  brk = round(length/5);
  start = koef / 6;
  stop = koef;
  l(1:brk) = start : (stop-start)/(brk-1) : stop;

  brk = round(length/40);
  start = koef;
  stop = koef / 2;
  l(end-brk:end) = start : (stop-start)/(brk) : stop;

  len = 2 * round(l / 2) + 1;
end

Funkce getLength vypočítá ze vstupních údajů adaptivní délku pro SG filtr.

Adaptivní délka má lichobežníkový půběh. V první části (1/5 celkové délky) lineárně stoupé z jedné šestiny na maximální délku. V poslední čtyřicetině pak lineárně kledá na jednu polovinu.

Maximlní délka (koef) pro adaptivní filtrování se vypočítá jako koef = 1.1 * (max(x) + max(y) - min(x) - min(y)) / 200.

Délka filtru pro filtrování v časové oblasti se vypočítá jako len = 2 * ceil(700 / mean([v' 10000])) + 1.

function [fx] = myfilt(x, len)
  length = max(size(x));
  for i = 1 : length
    if (i == 1) || (len(i) ~= len(i-1))
      b = sgolay(2, len(i));
    end

    row = (len(i) + 1) / 2; % center
    start = max(1, min(i - row, length - len(i)));
    if i < row
      row = i;
    elseif (i+row) > length
      row = len(i) - (length - i);
    end

    fx(i) = b(row, :) * x(start : start+len(i)-1)';
  end
end

Funkce myfilt provádí adaptivní filtrování. Pro každý prvek pole x je spočítána SG matice řádu 2 a délky len[i] vypočítané v minulé funkci. Touto maticí je potom přefiltrováno okolí x[i].

function [newX, newY] = myResample(oldX, oldY, distance)
  length = max(size(oldX,1), size(oldX,2));

  totalLength = 0;
  ds = zeros(1, length);
  for i = 2 : length
    ds(i) = sqrt((oldX(i-1)-oldX(i))*(oldX(i-1)-oldX(i)) + (oldY(i-1)-oldY(i))*(oldY(i-1)-oldY(i)));
    totalLength = totalLength + ds(i);
  end

  newX = zeros(1, floor(totalLength / distance));
  newY = zeros(1, floor(totalLength / distance));
  newX(1) = oldX(1);
  newY(1) = oldY(1);

  cds = 0;
  i = 1;
  j = 2;
  while i < length
    while (cds < distance) && (i < length)
      i = i + 1;
      cds = cds + ds(i);
    end
    while (cds >= distance) && (i <= length)
      a = 1 - (cds - distance)/ds(i);
      newX(j) = oldX(i-1) + a*(oldX(i) - oldX(i-1));
      newY(j) = oldY(i-1) + a*(oldY(i) - oldY(i-1));
      cds = cds - distance;
      j = j + 1;
    end
  end
  newX = newX(1:j-1);
  newY = newY(1:j-1);
return

Funkce resample přepočítá vektory oldX a oldY tak, aby dva následující body od sebe byly vzdáleny přesně distance.

Jméno souboru Test na Popis
_test.01 čáry

Převzorkování z času na vzdálenost. Lineární regrese.

Tiskne graf x-y a pruh 0.5cm široký, do kterého by se data měla vejít.

_test.02 spirály

Převzorkování z času na vzdálenost. Adaptivní filtrování.

Tiskne graf x-y a vyhlazené x-y, grafy x a vyhlazeného x v závislosti na délce.

_test.03 spirály

Filtrování v časové oblasti.

Tiskne graf x-y a vyhlazené x-y. FFT a spektrogram na rozdíl x, y.

_test.04 spirály

Filtrování v časové oblasti. Výpočet φ.

Tiskne graf x-y, závislodt φ na čase. FFT a spektrogram na rozdíl φ.

Data

Adresář Obsah
pcary
zcary
Data s rovnými čárami - horizontální, vertikální, diagonální
pspiral
zspiral
Data se spirálami