====== Řízení protetické ruky ======

===== Zadání =====
Realizujte ovládání protetické ruky na vývojové desce FRDM-K64F. K ovládání využijte EMG signál ze svalů na předloktí. K získání signálů sestrojte vhodný zesilovač a výstup přiveďte na AD převodníky. Vyhodnocení pohybu zobrazte na PC vhodnou animací.
===== Úvod =====
==== Popis problému ====
V moderních protézách se kromě lehčích a robustnějších materiálů hlavně ve formě karbonových vláken setkáváme stále častěji i s elektronikou. Myoelektrické protézy fungují na principu kontinuálního měření tzv. elektromyogramu (EMG), signálu, který vzniká ve svalu když se mění jeho konformace. Signál se snímá povrchovými elektrodami a když překročí danou mez, protéza se aktivuje. V roce 2013 se vyvinula první protéza ovládná čistě vůlí a permanentně upevněna na kost.

==== Popis realizace ====
Vstupní signály jsou snímany na šesti místech těla, vždy kladná, záporná a referenční svorka (3 elektrody pro biceps a 3 pro předloktí). To dovoluje získat a vyhodnotit dva druhy pohybu. Jako předzesilovače slouží dvě měřící jednotky MP-35 od firmy Biopac. Zesílené nefiltrované signály jsou přivedeny na analogové vstupy procesoru. V procesoru se EMG signály vzorkují frekvencí 1 kHz. Filtrace je realizována pomocí IIR filtrů. Výsledné vyhodnocení je posíláno přes sériovou linku do PC, kde jsou míry zvednutí a otočení ruky zobrazeny v čase.
===== Hardware =====
Umístění elektrod lze názorně vidět na obr. 1. Snímané signály jsou v rozsahu od -1,5 po 1,5 mV. Jsou přivedeny na vstupní kanál měřící jednotky MP-35, ta je zesílí 2000-krát, takže na výstupu se pohybují hodnoty od -3V až 3V, které se přivedou na vstupy AD převodníku procesoru FRDM-K35F, kanály A0 a A5 a na osciloskop pro vizualizaci. Měřící jednotka má pouze jeden analogový výstupní kanál, takže bylo potřeba využít dvě jednotky, jedna pro snímání signálu z bicepsu, druhá z předloktí. K jednotkám je dodáván software Student Lab, kde ve verzi Pro lze libovolně nastavit zesílení i filtraci vstupních dat. Měřená data jsou zobrazována na počítači.

Specifikace použitého procesoru je možné najít například zde: [[https://developer.mbed.org/platforms/FRDM-K64F/]]

{{ :2015:bio-arm:elektrody1.jpg |umístění elektrod}}
    Obr. 1: Umístění elektrod
{{ :2015:bio-arm:elektrody.jpg |Ukázka elektrody}}
    Obr. 2: Ukázka elektrody
{{ :2015:bio-arm:biopack_ekg.jpg?400 |Měřící jednotka}}
    Obr. 3: Měřící jednotka
{{ :2015:bio-arm:procesor.png |Vývojový kit s FRDM-K64F }}
    Obr. 4: Vývojový kit s FRDM-K64F 
{{ :2015:bio-arm:Set_popis.jpg |Pracovní plocha }}
    Obr. 5: Pracovní plocha     
{{ :2015:bio-arm:emg.png |Ukázka naměřeného EMG signálu }}
    Obr. 6: Ukázka naměřeného EMG signálu    
===== Software ===== 
==== Filtrace ====   
Filtrace obou EMG signálů je realizována pomocí IIR charakteristice založené na Butterworth pásmové propusti s mezními frekvencemi 20 a 50 Hz. Čtení i filtrování je po jednom vzorku v každé iteraci smyčky ve které je nastaven wait na 1 ms čímž je určena vzorkovací frekvece. Neni to tedy přesně 1 kHz, ale vzhledem k jednoduchosti příkazů ve smyčce dostatečně blízko této hodnotě. 

Funkce pro filtraci byly vygenerovány na stránce [[http://www.micromodeler.com/dsp/]]. Je zde možné vybrat typ filtru, různými způsoby upravit jeho charakteritisku a pak použít vygenerovanou hlavičku a funkce. Hlavičkový soubor se pouze připojí a funkce jsou nakopírovány na konci main.cpp programu. Ve free verzi je u IIR filtrů možné generovat kód pouze pro filtry do 4. řádu. 

Odkaz na interaktivní vytváření filtrů Micro modeler byl získán ze stránky [[https://community.arm.com/groups/embedded/blog/2014/02/04/introduction-to-digital-filters-2]] odkud taky pochází inspirace pro jejich použití. 
{{ :2015:bio-arm:prenosovka.png?800 |Přenosová charakteristika filtru}}
    Obr. 7: Přenosová charakteristika filtru
==== Učící fáze ==== 
Při spuštění programu má uživatel 8 vteřin na provedení jednoho nebo více reprezentativních přitažení předloktí k rameni. Po načtení dat proběhne výpočet obálky signálu, v okně 100 vzorků se počítá průměrná hodnota horní poloviny hodnot v okně. K této hodnotě obálky se přičte 1 a umocní na druhou (1 se přičte aby při umocnění všechny hodnoty vzrostly). Okno se posune o 1 vzorek a iterace se opakuje. Program dále vyhodnotí práh jako 25% hodnot obálky, podle kterého pozná, zda uživatel přitahuje ruku či nikoliv.

Obdobně se postupuje u signálu z předloktí: po 8 vteřinách učící fáze bicepsu je vteřina pauza signalizovaná led a poté zahájena učící fáze 8 vteřin pro předloktí. 
<code c>
  
//########################################################
//## ucici phase emg4 - biceps 8 vterin ... cca 4 flexy ##
//########################################################
led = 1;                //zhasnuti led
led2 =1;
led3 = 1;
pc.baud(115200);        //nastaveni baudrate

for(int i=0;i<(N);i++){          //ulozeni hodnoty do vzorek, filtrace a abs. hodnota
    vzorek[0] = ain.read_u16();
    nProcessedSamples = filter1_filterBlock( filter, vzorek, outputBuffer, inputBufferSize ); 
    emg4abs[i]=abs(outputBuffer[0]);
    wait(0.001f);                //fvz = 1kHz
}
led = 0;                //roznuti cervene led        

for(int i=0;i<N-M;i++){ //smycka pro vypocet obalky emg4 ucici faze
    for(int j=0;j<M;j++){
        Maktual[j]=emg4abs[i+j];  //nacteni do okna M vzorku (default 100)
    }
    qsort(Maktual,M,sizeof(float),compare); //serazeni od nejmensiho vzorku v okne po nejvetsi
    for(int j=(M/2);j<M;j++){               
        sum4=sum4+Maktual[j];    //secteni hornich defaultne 50 vzorku okna
    }
    emg4obal[i]=sum4/(M/2);      //vydeleni poctem defaultne 50 vzorku okna
    sum4=0;                      //tj. prumerna hodnota horni poloviny hodnot okna tvori obalku
}

for (int j=0;j<N-M;j++){
    emg4obal[j]=emg4obal[j]+1;                  //pricteni 1 pro omezeni hodnot mensich nez 1
    emg4obal[j]=emg4obal[j]*emg4obal[j];        //umocneni obalky (vsechny hodnoty se tak zvysi ptze jsou vetsi nez 1)
}
qsort(emg4obal,N-M,sizeof(float),compare);      //serazeni od min po max

meze2=0;                                        //inicializace meze pro biceps
for(int j=20;j<220;j++){                        //vyber 20. az 220. maxima
    maxs1[j-20]=emg4obal[(N-M)-(1*j)-1];        //a ulozeni do maxs1
    meze2=meze2+(maxs1[j-20]);                  //suma 200 submaxim
}      
meze2=meze2/200/4;      //podeleni poctem maxim a snizeni na ctvrtinu = definice prahu
</code>
==== Vyhodnocovací fáze ====
V této fázi program začíná v nekonečné smyčce kontinuálně meřit a vyhodnocovat naměřená data s oknem 100 vzorků. Na datech je provedena obálka podobně jako v učící fázi. Následuje kritérium, kdy se získaná hodnota porovná s prahy. Jsou-li hodnoty větší, rozsvítí se příslušná led a do čítače s rozsahem 0 až 400 hodnot pro oba signály se přičte 1 pokud již nejsou na maximální hodnotě. Jsou-li hodnoty menší, led zhasne a z čítače se odečte 1 pokud není hodnota 0. Oba signály používají stejnou proměnnou (procenta), biceps je kodován po jednotkách, předloktí po tisících. Reálná maximální hodnota je tak 400400.

Jelikož při zvedání ruky bicepsem vzniká signál, který je relativně dobře čitelný i na svodech předloktí, je vytvořena pojistka, při které se vyhodnocuje zda je biceps aktivní. Je-li aktivní, pak se práh pro aktuální vzorek otočení zvýší trojnásobně.
<code c>
led=1;                            //zhasnuti led indikujici start mereni
//####################################################
//## kriterium emg4 - biceps, kriterium emg3 - otoc ##
//####################################################

for(int i=0;i<M;i++){
    vzorek[0]=ain.read_u16();
    nProcessedSamples = filter1_filterBlock( filter, vzorek, outputBuffer, inputBufferSize ); 
    emg4[i]=abs(outputBuffer[0]); 
    
    vzorek[0]=otoc.read_u16();
    nProcessedSamples = filter1_filterBlock( filter, vzorek, outputBuffer, inputBufferSize ); 
    emg3[i]=abs(outputBuffer[0]); 
    wait(0.001f); 
}
for (int i=0;i<M;i++){
    data4[i]=emg4[i];       //ulozeni do data4 a data3 ktere se quicksortuje
    data3[i]=emg3[i];
}
while(1)                    //nekonecna smycka mereni
{
    qsort(data4,M,sizeof(float),compare);      //serazeni data4 od min po max
    qsort(data3,M,sizeof(float),compare);      //serazeni data3 od min po max
    for(int j=(M/2);j<M;j++){
        sum4=sum4+data4[j];                    //secteni horni poloviny dat
        sum3=sum3+data3[j];
    }
    obal4=sum4/(M/2);                          //prumer horni poloviny dat
    obal4=obal4+1;                                      
    obal4=obal4*obal4;                         //mocninne zvyseni dat
    sum4=0;                                    //reset sumy
    
    obal3=sum3/(M/2);                          //to same s predloktim
    obal3=obal3+1;
    obal3=obal3*obal3;
    sum3=0;
    
    if(meze2<obal4){            //jestlize bude hodnota obalky vetsi nez meze2               
        emg4krit=1;             //kriterium pro biceps = 1 ... ruka se zveda
        if(procenta<400){       //jestlize neni na max hodnote 400 pak se zvedne
           procenta++;  
        }
        led2=0;                 //zelena led signalizujici zvedani
    }else{
        emg4krit=0;
        if(procenta>0){         //jestlize neni na min hodnote 0 pak klesne
           procenta--;
        }  
        led2=1;                 //zelena led signalizuje pokles
    }
    
    meze1zaloha=meze1;          //zalohuje se puvodni meze
    if(emg4krit==1){            //jestlize se ruka zveda                 
        meze1=meze1*zvysenimeze;//a zvysi se prah (defaultne 3*)                       
    }
    
    if(meze1<obal3){            //jestlize bude hodnota obalky vetsi nez meze1     
        if(procenta<400400){    //jestlize neni na max hodnote
           procenta=procenta+1000; //tak se zvysi (mod 1000)
        }   
        led3=0;                 //modra led signalizujici otoceni   
    }else{
        led3=1;  
        if(procenta>400){      //jestlize neni na min hodnote mod 1000 0 pak klesne
           procenta=procenta-1000;
        }  
    }
    meze1=meze1zaloha;         //navrat na puvodni mez                  
    
    if (cyklus % 10 == 1){           //po seriove lince se posila kazda 10. hodnota
        pc.printf("%d\n",procenta);  //vypis procent po seriove lince
    } 
       
    for (int i=0;i<M-1;i++){   //posun v okne o 1 dozadu      
        emg4[i]=emg4[i+1];
        emg3[i]=emg3[i+1];
    }
    
    vzorek[0]=ain.read_u16();  //ulozeni vyfiltrovane hodnoty na posledni index okna
    nProcessedSamples = filter1_filterBlock( filter, vzorek, outputBuffer, inputBufferSize ); 
    emg4[M-1]=abs(outputBuffer[0]);
    vzorek[0]=otoc.read_u16();
    nProcessedSamples = filter1_filterBlock( filter, vzorek, outputBuffer, inputBufferSize ); 
    emg3[M-1]=abs(outputBuffer[0]); 
    wait(0.001f); 
    
    cyklus++;                  //dalsi cyklus     
    
    for (int i=0;i<M;i++){     //ulozeni do data4 a data3 ktere se quicksortuje
        data4[i]=emg4[i];
        data3[i]=emg3[i];
    }
}
</code>
{{ :2015:bio-arm:emgOO.png |Ukázka vyhodnocení dle prahu v MATLABu, hodnota 20 odpovídá aktivnímu otáčení při aktivním bicepsu}}
    Obr. 8: Ukázka vyhodnocení dle prahu v MATLABu, hodnota 20 odpovídá aktivnímu otáčení při aktivním bicepsu
==== Zobrazení stavu ====
Na sériovou linku se posílá každá 10. hodnota stavu. Ta je načtena programem "serialchart" pro vykreslování hodnot v čase. Program je volně dostupný z [[https://code.google.com/p/serialchart/]]. Data lze rovněž podobným způsobem vyhodnocovat pomocí MATLABu: 
<code matlab>
s = serial('COM3');
set(s, 'InputBufferSize', 1);
set(s, 'FlowControl', 'hardware');
set(s, 'BaudRate', 115200);
set(s, 'Parity', 'none');
set(s, 'DataBits', 8);
set(s, 'StopBit', 1);
set(s, 'Terminator', 'CR');
set(s, 'Timeout',5);
 
fopen(s);          
t=1;
x=0;
while(1) 
   a =fread(s);
   a=max(a); 
   x =[x a];  
   if length(x)>201
   plot(mod(x(end-200:end),1000));                     %vykresleni bicepsu (modulo 1000)
   hold on                                             %napr. z 66100 udela 100
   plot((x(end-200:end)-mod(x(end-200:end),1000)/1000),'r');  %vykresleni predlokti
   hold off                                            %od hodnoty 66100 se odecte 100 a podeli 1000
   axis auto;
   grid on;
   disp([num2str(t),'th iteration max= ',num2str(a)]);
   t=t+1;
   a=0; 
   drawnow;
   end
end
fclose(s);  
</code>

Zdrojový kód byl napsán pomocí vývojového prostředí mbed.org a je dostupný na [[https://developer.mbed.org/users/customer10123/code/Rizeni_ruky_K64F/]]
===== Demonstrační video a diskuze=====

{{youtube>9D8jCWQ_3tU?medium}}

Ve videu je možné shlédnout všechny fáze programu, na monitorech signály vstupující do procesoru z bicepsu a předloktí a křivku vyhodnocení stavu funkce bicepsu v čase. Místo rotace ruky lze použít stisk, u kterého je výsledný signál více zřetelný v závislosti na pozici elektrod.

Učící fáze je naprostou samozřejmostí, protože každý uživatel má hodnoty energie EMG signálu odlišné. Signály jsou filtrovány stejným filtrem. V kódu se za sebou pro filtraci střídají jednotlivé hodnoty signálu z bicepsu a předloktí. Původní představa toho, že se hodnoty jednotlivých vstupních kanálů budou mezi sebou hádat se ukázala jako mylná a oba výsledné vyfiltrované signály se jevily jako vyfiltrované skutečně správně. Nemusela se tak psát kopie filtru, který by filtroval druhý kanál, jak se původně smýšlelo.

Hodnotí se pouze zda-li je sval aktivní/neaktivní a po jakou dobu. Na základě těchto stavů lze pouze hodnotit zda se končetina hýbe k maximální nebo minimální vychýlce tzn. ustálený stav je pouze v maximu nebo minimu. K vychýlce 50% by uživatel musel kontinuálně daný sval stahovat a uvolňovat, výsledný charakter ustáleného pohybu se však jeví spíš jako vibrující kolem 50%, než ustálený stav. Tento postup vyhodnocování se však zdál nejvíce schůdný vzhledem k nevyspytatelnosti EMG signálu. Rychlost pohybu k maximu a minimu je rovněž konstantní, avšak v kódu lehce měnitelná vzhledem k použití operátoru %.
===== Závěr =====
Cílem projektu bylo realizovat ovládání protetické ruky s pomocí kitu FRDM-K64F a vhodně zobrazovat její pohyb na PC. Výsledný projekt je schopen detekovat EMG signál ze dvou míst a na základě jejich krátkodobých hodnot v čase vyhodnotit a realizovat dva různé pohyby. O tom v jakém stavu vychýlky v čase se zrovna končetina nachází zobrazuje křivka měřených hodnot na PC. Během návrhu se pracovalo se signály z bicepsu a předloktí. V praxi by se však použilo k ovládání jiné místo než předloktí, protože člověk, který by chtěl realizovat práci bicepsu by pravděpodobně předloktí neměl. Vhodným vylepšením by bylo vytvořit vlastní zesilovače, bylo by tak možné měření provádět kdekoliv s přístupem napájení. Celý algoritmus by se dal dále upravovat k lepšímu a plynulejšímu ovládání.