2010/bmi10/EMG: BayesEMG.m

File BayesEMG.m, 3.8 KB (added by adutta1, 8 years ago)

A NONLINEAR FILTER FOR ONLINE BAYESIAN ESTIMATION OF MYOELECTRIC SIGNALS

Line 
1% % A NONLINEAR FILTER FOR ONLINE BAYESIAN ESTIMATION OF MYOELECTRIC SIGNALS
2% %
3% %
4% % Terence D. Sanger, MD PhD
5% % Div. Child Neurology and Movement Disorders
6% % Stanford University Medical Center
7% % 300 Pasteur, room A345
8% % Stanford, CA 94305-5235 USA
9% % (650)736-2154  fax: 725-7459
10% % sanger@stanford.edu
11
12
13%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
14%
15% Preliminary section: load data, set constants, initialize variables
16%
17%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
18
19%set parameters
20samplerate = 1000;  %samples per second
21noutputs = 100;      %output quantization levels
22ratemax = 1;        %rectified EMG is normalized to max value of 1
23inscale = 4;        %arbitrary input scaling
24alpha = 0.0001 / samplerate;                %sets diffusion rate
25beta = 1e-24 / (noutputs * samplerate);     %sets probability of sudden jumps
26
27%load the data
28load('E:\Users\w4i4v\TMR_01_Reformatted.mat');    %datafile has 3 columns:  torque, bicepsEMG, tricepsEMG
29biceps =Data(:,4)';
30%triceps = v(3,:)';
31%torque = v(1,:)';
32
33%calculate rectified EMG after removing the mean, and normalize
34%here, we use only biceps and thus cannot approximate extensor torque
35emg = abs(biceps - mean(biceps));               %essential to remove the mean before rectifying
36emg = inscale * ratemax * emg / max(emg);       %input prescaling to use full output range
37emg(emg>ratemax) = ratemax;                     %make sure we don't go over
38
39%initialize variables
40%   x is the latent variable (the driving rate)
41%   MAP is the output estimate
42x = linspace(ratemax/noutputs, ratemax, noutputs)';  %don't start with zero because requires n=0 exactly to match
43MAP = zeros(length(emg),1);                     %store the bayes estimates
44g = [(alpha/2) (1 - alpha) (alpha/2)];          %approximate spatial second derivative operator
45
46%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
47%
48% Following is the main section of the algorithm; steps are numbered as in the text
49%
50%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
51
52%1.     Initialize p(x,0) = 1;
53prior = ones(noutputs,1) / noutputs;            %start with uniform prior
54
55for t=1:length(emg) %iterate for each sample of EMG
56   
57    %2. Forward propagate p(x,t-) Å
58    %           ?p(x-?,t-1)+(1-2?)p(x,t-1)+?p(x+?,t-1)+?+(1-?)p(x,t-1);
59    prior = filtfilt(g, 1, prior);              %drift term by convolving with second derivative operator
60    prior = beta + (1-beta) * prior;            %sets probability of a sudden jump
61   
62    %3. Measure the rectified emg;
63    emgval = emg(t);                            %if this were online, would read a new sample here
64   
65    %4. Calculate the posterior likelihood function
66    %       P(x,t) Å P(emg|x)p(x,t-);
67    measurement_model = exp(-emgval./x) ./ x;   %exponential model for P(emg|x)
68    posterior = measurement_model .* prior;     %calculate posterior density using Bayes rule
69   
70    %5. Output the signal estimate MAP(x(t)) = argmax P(x,t);
71    pp = min(find(posterior == max(posterior)));    %find the maximum of the posterior density
72    if (pp > 1 && pp < length(posterior)),          %interpolate to find the zero
73        dL = posterior(pp-1) - posterior(pp); 
74        dR = posterior(pp) - posterior(pp+1);
75        PeakIndex = (pp - .5 - (dL / (dR - dL)));   %index runs from 1 to noutputs
76    else
77        PeakIndex = pp;    %if maximum occurs at an endpoint do not interpolate
78    end
79    MAP(t) = (ratemax / (noutputs-1)) * PeakIndex;  %convert index of peak value to scaled EMG value
80   
81    %6. Divide p(x,t) by a constant C so that  º p(x,t) dx = 1;
82    posterior = posterior / sum(posterior);        %normalize the density
83   
84    %7. Repeat from step 2;
85    prior = posterior;                          %prior for next iteration is posterior from this iteration
86end
87
88%show results
89figure
90hold on
91%plot(torque/max(torque));
92plot(MAP/max(MAP),'r');
93hold off
94