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.
koef
je maximální délka filtrulength
je délka výsledného vektoruAdaptivní 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 φ. |
Adresář | Obsah |
---|---|
pcary zcary |
Data s rovnými čárami - horizontální, vertikální, diagonální |
pspiral zspiral |
Data se spirálami |