人行走的微多普勒效应的MATLAB仿真.zip

function [segment,seglength,T] = HumanWalkingModel(showplots, formove,... animove, genmovie, Height, rv, nt, numcyc) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Global Human Walking Model % % Based on "A global human walking model with real-time kinematic % personification", by R. Boulic, N.M. Thalmann, and D. Thalmann % The Visual Computer, vol.6, pp.344-358, 1990 % This model is based on biomechanical experimental data. % % By V.C. Chen and Yang Hai % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % total number of pulses np = nt*numcyc; fprintf('Total Number of Samples = %g\n',np); % spatial characteristics rlc = 1.346*sqrt(rv); % relative length of a cycle % temporal characteristics dc = rlc/rv; % duration of a cycle ds = 0.752*dc-0.143; % duration of support dsmod = ds/dc; % relative duration of support T = dc*numcyc; % total time duration % body segments' length (meter) headlen = 0.130*Height; shoulderlen = (0.259/2)*Height; torsolen = 0.288*Height; hiplen = (0.191/2)*Height; upperleglen = 0.245*Height; lowerleglen = 0.246*Height; footlen = 0.143*Height; upperarmlen = 0.188*Height; lowerarmlen = 0.152*Height; Ht = upperleglen + lowerleglen; sprintf('The Walking Velocity is %g m/s',rv*Ht); % time scaling dt = 1/nt; t = [0:dt:1-dt]; t3 = [-1:dt:2-dt]; % designate gait characteristic if rv <= 0 error('Velocity must be positive') elseif rv < 0.5 gait = 'a'; elseif rv < 1.3 gait = 'b'; elseif rv <= 3 gait = 'c'; else error('Relative velocity must be less than 3') end % Locations of body segments: % 3 translation trajectory coords. give the body segments location % relative to rv % calculate vertical translation: offset from the current height (Hs) of % the origin of the spine Os av = 0.015*rv; verttrans = -av+av*sin(2*pi*(2*t-0.35)); maxvl = max(verttrans); minvl = min(verttrans); diffvl = maxvl-minvl; if showplots == 'y' figure plot(t,verttrans) title('Vertical translation offset from the origin of the spine') xlabel('gait cycle') ylabel([sprintf('Position (m) [V_R = %g]',rv)]) axis([0 1 minvl-0.2*diffvl maxvl+0.2*diffvl]) drawnow end % calculate lateral translation: Os oscillates laterally to ensure the % weight transfer from one leg to the other. if gait == 'a' al = -0.128*rv^2+0.128*rv; else al = -0.032; end lattrans = al*sin(2*pi*(t-0.1)); maxvl = max(lattrans); minvl = min(lattrans); diffvl = maxvl-minvl; if showplots == 'y' figure plot(t,lattrans) title('Lateral translation offset from the origin of the spine') xlabel('gait cycle') ylabel([sprintf('Offset Position (m) [V_R = %g]',rv)]) axis([0 1 minvl-0.2*diffvl maxvl+0.2*diffvl]) drawnow end % calculate translation forward/backward: acceleration and deceleration % phases. When rv grows this effect decreases. if gait == 'a' aa = -0.084*rv^2+0.084*rv; else aa = -0.021; end phia = 0.625-dsmod; transforback = aa*sin(2*pi*(2*t+2*phia)); maxvl = max(transforback); minvl = min(transforback); diffvl = maxvl-minvl; if showplots == 'y' figure plot(t,transforback) title('Translation forward/backward from the origin of the spine') xlabel('gait cycle') ylabel([sprintf('Offset Position (m) [V_R = %g]',rv)]) axis([0 1 minvl-0.2*diffvl maxvl+0.2*diffvl]) drawnow end % two rotations of the pelvis % calculate rotation forward/backward: to make forward motion of the leg, % the center of gravity of the body must move. To do this, flexing movement % of the back relatively to the pelvis must be done. if gait == 'a' a1 = -8*rv^2+8*rv; else a1 = 2; end rotforback = -a1+a1*sin(2*pi*(2*t-0.1)); maxvl = max(rotforback); minvl = min(rotforback); diffvl = maxvl-minvl; if showplots == 'y' figure plot(t,rotforback) title('Rotation Forward/Backward') xlabel('gait cycle') ylabel([sprintf('Rotation forward/backward (degree) [V_R = %g]',rv)]) axis([0 1 minvl-0.2*diffvl maxvl+0.2*diffvl]) drawnow end % calculate rotation left/right: the pelvis falls on th side of the % swinging leg. a2 = 1.66*rv; temp1 = -a2+a2*cos(2*pi*(10*t(1:round(nt*0.15))/3)); temp2 = -a2-a2*cos(2*pi*(10*(t(round(nt*0.15)+1:round(nt*0.50))-0.15)/7)); temp3 = a2-a2*cos(2*pi*(10*(t(round(nt*0.50)+1:round(nt*0.65))-0.5)/3)); temp4 = a2+a2*cos(2*pi*(10*(t(round(nt*0.65)+1:nt)-0.65)/7)); rotleftright = [temp1,temp2,temp3,temp4]; maxvl = max(rotleftright); minvl = min(rotleftright); diffvl = maxvl-minvl; if showplots == 'y' figure plot(t,rotleftright) title('Rotation Left/Right') xlabel('gait cycle') ylabel([sprintf('Rotation left/right (degree) [V_R = %g]',rv)]) axis([0 1 minvl-0.2*diffvl maxvl+0.2*diffvl]) drawnow end % calculate torsion rotation: pelvis rotates relatively to the spine to % perform the step a3 = 4*rv; torrot = -a3*cos(2*pi*t); maxvl = max(torrot); minvl = min(torrot); diffvl = maxvl-minvl; if showplots == 'y' figure plot(t,torrot) title('Torsion Rotation') xlabel('gait cycle') ylabel([sprintf('Torsion angle (degree) [V_R = %g]',rv)]) axis([0 1 minvl-0.2*diffvl maxvl+0.2*diffvl]) drawnow end % leg flexing/extension: at the hip, at the knee, and at the ankle. % calculate flexing at the hip if gait == 'a' x1 = -0.1; x2 = 0.5; x3 = 0.9; y1 = 50*rv; y2 = -30*rv; y3 = 50*rv; end if gait == 'b' x1 = -0.1; x2 = 0.5; x3 = 0.9; y1 = 25; y2 = -15; y3 = 25; end if gait == 'c' x1 = 0.2*(rv-1.3)/1.7-0.1; x2 = 0.5; x3 = 0.9; y1 = 5*(rv-1.3)/1.7+25; y2 = -15; y3 = 6*(rv-1.3)/1.7+25; end if x1+1 == x3 x = [x1-1,x2-1,x1,x2,x3,x2+1,x3+1]; y = [y1,y2,y1,y2,y3,y2,y3]; else x = [x1-1,x2-1,x3-1,x1,x2,x3,x1+1,x2+1,x3+1]; y = [y1,y2,y3,y1,y2,y3,y1,y2,y3]; end temp = pchip(x,y,t3); flexhip = temp(nt+1:2*nt); maxvl = max(flexhip); minvl = min(flexhip); diffvl = maxvl-minvl; if showplots == 'y' figure plot(t,flexhip) hold on plot(x1,y1,'ro',x2,y2,'ro',x3,y3,'ro') title('Flexing at the Hip') xlabel('gait cycle') ylabel([sprintf('Flexing angle (degree) [V_R = %g]',rv)]) if x1 < 0 axis([x1 1 minvl-0.2*diffvl maxvl+0.2*diffvl]) else axis([0 1 minvl-0.2*diffvl maxvl+0.2*diffvl]) end drawnow end % calculate flexing at the knee: there are 4 control points. if gait == 'a' x1 = 0.17; x2 = 0.4; x3 = 0.75; x4 = 1; y1 = 3; y2 = 3; y3 = 140*rv; y4 = 3; end if gait == 'b' x1 = 0.17; x2 = 0.4; x3 = 0.75; x4 = 1; y1 = 3; y2 = 3; y3 = 70; y4 = 3; end if gait == 'c' x1 = -0.05*(rv-1.3)/1.7+0.17; x2 = 0.4; x3 = -0.05*(rv-1.3)/1.7+0.75; x4 = -0.03*(rv-1.3)/1.7+1; y1 = 22*(rv-1.3)/1.7+3; y2 = 3; y3 = -5*(rv-1.3)/1.7+70; y4 = 3*(rv-1.3)/1.7+3; end x = [x1-1,x2-1,x3-1,x4-1,x1,x2,x3,x4,x1+1,x2+1,x3+1,x4+1]; y = [y1,y2,y3,y4,y1,y2,y3,y4,y1,y2,y3,y4]; temp = pchip(x,y,t3); flexknee = temp(nt+1:2*nt); maxvl = max(flexknee); minvl = min(flexknee); diffvl = maxvl-minvl; if showplots == 'y' figure plot(t,flexknee) hold on plot(x1,y1,'ro',x2,y2,'ro',x3,y3,'ro',x4,y4,'ro') title('Flexing at the Knee') xlabel('gait cycle') ylabel([sprintf('Flexing angle (degree) [V_R = %g]',rv)]) if x1 < 0 axis([x1 1 minvl-0.2*diffvl maxvl+0.2*diffvl]) else axis([0 1 minvl-0.2*diffvl maxvl+0.2*diffvl]) end drawnow end % calculate flexing at the ankle: there are 5 control points if gait == 'a' x1 = 0; x2 = 0.08; x3 = 0.5; x4 = dsmod; x5 = 0.85; y1 = -3; y2 = -30*r