# 2010/bmi10/EMG: BayesEMG.m

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

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 |

20 | samplerate = 1000; %samples per second |

21 | noutputs = 100; %output quantization levels |

22 | ratemax = 1; %rectified EMG is normalized to max value of 1 |

23 | inscale = 4; %arbitrary input scaling |

24 | alpha = 0.0001 / samplerate; %sets diffusion rate |

25 | beta = 1e-24 / (noutputs * samplerate); %sets probability of sudden jumps |

26 | |

27 | %load the data |

28 | load('E:\Users\w4i4v\TMR_01_Reformatted.mat'); %datafile has 3 columns: torque, bicepsEMG, tricepsEMG |

29 | biceps =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 |

35 | emg = abs(biceps - mean(biceps)); %essential to remove the mean before rectifying |

36 | emg = inscale * ratemax * emg / max(emg); %input prescaling to use full output range |

37 | emg(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 |

42 | x = linspace(ratemax/noutputs, ratemax, noutputs)'; %don't start with zero because requires n=0 exactly to match |

43 | MAP = zeros(length(emg),1); %store the bayes estimates |

44 | g = [(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; |

53 | prior = ones(noutputs,1) / noutputs; %start with uniform prior |

54 | |

55 | for 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 |

86 | end |

87 | |

88 | %show results |

89 | figure |

90 | hold on |

91 | %plot(torque/max(torque)); |

92 | plot(MAP/max(MAP),'r'); |

93 | hold off |

94 |