stepdetection.m

一个matlab实例程序

%function Acc = process(fname)
fname='anttifootgai.TXT';

STARTP = 2000; % start position of the plotted TS window
ENDP = 3000;   % end position of the plotted TS window
SCA3000 = 1333;  % SCA3000 sensitivity 

% read input data file
fid = fopen(fname);
data = fgetl(fid);
if length(data) < 26
 for i = 1:26-length(data)
  data = [data,' '];
 end;
end;

while ~feof(fid)
  l = fgetl(fid);
  data = [data;l];
end;

% cell array B contains raw acceleration values in hex format
B = cell(length(data),4);        
for i = 1: length(data)
Aux = textscan(data(i,:),'%s','delimiter',',');
  for j = 1:4
    B{i,j} = Aux{1}{j};
  end;
end;

Ln = length(B);

% array Acc contains accelerations converted to the g-scale
Acc = zeros(Ln,3);

% VTI conversion algorithm
  Threshold = 2^16/2;
for i = 1:Ln
 for j = 1:3
  if length(B{i,j}) == 1
    val = hex2dec(num2str(B{i,j}(1:1)));
  elseif length(B{i,j}) == 3
    val = hex2dec(num2str(B{i,j}(3:3)));
  else
    val = hex2dec(num2str(B{i,j}(3:6)));
  end;
  if val >= Threshold
    Acc(i,j) = -(2^16 - val)/8;
  else
   Acc(i,j) = val/8;
  end;
  Acc(i,j) = Acc(i,j)/SCA3000*1000;
 end;  
end;

%plotting 
L = length(Acc);
T = 1:L;
T = T ./ 200;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%产生需要的数组形式
Lv=length(Acc);
numd=zeros(Lv,1);%差分处理后的numd为3D合成加速度
numa=zeros(Lv,1);%处理后的numa为3D合成加速度
numb=zeros(Lv,1);%处理后的numb为3D合成加速度

for i=1:length(Acc)
    numd(i,1)=((Acc(i,1))^2+(Acc(i,2))^2+(Acc(i,3))^2)^(0.5);
    numd(i,1)=numd(i,1)./1000;
    numa(i,1)=numd(i,1);
    numb(i,1)=numd(i,1);
end;

count=0;%初始化计数变量，用于统计周期个数
max=1000;%最高统计周期个数
time=zeros(max,1);%周期统计数组，统计每个跨步所用时间，最多统计,max个。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%预处理


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%判断信号类型
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%一种方法是互相关
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%一种方法是方差
kind=0;%初始化kind变量用于表示类型判断的结果
%kind=1;%表明为脚部信号
%kind=2;%表明为头部信号
%kind=3;%表明为肩部信号
%kind=4;%表明为荷包处信号
va=var(numd(:,1));
if va>=0.1
    kind=1;
else
    kind=2;
end;
  


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%画出初始图形
%figure;
%set(gca,'FontSize',14);
%plot(numb(:,1),'g-');
%xlabel('[n]');%横坐标使用点数作为计数，200个点表示一秒钟/200还是100个点每秒有疑问
%ylabel('[3D]');
%title('3D Acc');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%如果是脚部信号
if kind==1
    figure;
    set(gca,'FontSize',14);
    plot(numb(:,1),'g-');
    xlabel('[n]');%横坐标使用点数作为计数，200个点表示一秒钟/200还是100个点每秒有疑问
    ylabel('[3D]');
    title('3D Acc');
    i=1;%控制变量
    j=1;%用于控制对于周期统计的赋值
    k=1;%用于寻找一步的开始点
    m=0.025;%限制统计和的范围
    n=0.1;%限制加速度值的范围
    p=0;%用于标志一个跨步行为是否被检测到了
    q=0.5;%用于标示一个跨步的起始点
    l1=0;%用于标示一步的起始点横坐标
    l2=0;%用于标示运动的结束点横坐标
    l3=0;%用于标示一个周期的结束点的横坐标
    l4=0;%用于传值
    while i<Lv-1
        if i<10
            i=i+1;
        else
            a=(abs(numb(i-9,1)-1)+abs(numb(i-8,1)-1)+abs(numb(i-7,1)-1)+abs(numb(i-6,1)-1)+abs(numb(i-5,1)-1)+abs(numb(i-4,1)-1)+abs(numb(i-3,1)-1)+abs(numb(i-2,1)-1)+abs(numb(i-1,1)-1)+abs(numb(i,1)-1))./10;
            if a>q&p==0&numb(i,1)>1
                p=1;
                k=i;
                while numb(k,1)>1
                    k=k-1;
                end;
                if k<l4
                    k=l4;
                end;
                text(k,numb(k,1),' \leftarrow 一步的起始点','FontSize',12);
                l1=k;
                i=i+1;
            end;
            if a<m&numb(i,1)<(1+n)&numb(i,1)>(1-n)&p==1
                text(i,numb(i,1),' \leftarrow 进入stance阶段','FontSize',12);
                l2=i;
                while a<m&numb(i,1)<(1+n)&numb(i,1)>(1-n)&p==1&i<Lv&i<(2.*l2-l1)
                    i=i+1;
                end;
                l3=i-1;
                if l2>l1&l3>l2&(2.*(l3-l2))>(l2-l1)%如果走完这一步后经过了一段时间停止
                    l3=fix(l2+(l2-l1)./3);
                    l4=l3;
                    text(l3,numb(l3,1),' \leftarrow 一个周期的结束','FontSize',12);
                    count=count+1;
                    time(j,1)=(l3-l1)./200;
                    j=j+1;
                    i=i+1;
                    p=0;
                    l1=0;
                    l2=0;
                    l3=0;
                else
                    l3=fix(l2+(l2-l1)./3);
                    l4=l3;
                    text(l3,numb(l3,1),' \leftarrow 一个周期的结束点','FontSize',12);
                    count=count+1;
                    time(j,1)=(l3-l1)./100;
                    j=j+1;
                    i=i+1;
                    p=0;
                    l1=0;
                    l2=0;
                    l3=0;
                end;
            else
                i=i+1;
            end;
        end;
    end;
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%如果是其余信号
elseif kind==2
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%使用峰值探测所需要的预处理
    n=0;
    L=600;%预处理中取600个点,代表6秒
    for i=(1+L):(Lv-L)
        k=i-L;
        for k=(i-L):(i-1)
            n=n+numd(k,1);
        end;
        n=n/L;
        numa(i,1)=numd(i,1)-n;
        n=0;
    end;
    figure;
    set(gca,'FontSize',14);
    plot(numa((1+L):(Lv-L),1),'r-');
    xlabel('[sec]');
    ylabel('[3D]');
    title('3D disposaled acceleration');
    m=0;
    i=1;
    for i=(1+L):(Lv-L)%取L是代表6秒钟之内的点
        if numa(i,1)>m
            m=numa(i,1);
        end;
    end;  
    peak=zeros(Lv,2);%定义这组数据是为了存储整个信号上所有的极大值点，下一步即从里面选出peak的极大值点
    k=1;
    for i=(1+L):(Lv-L)%取601和6000是因为上面预处理的时候是对这一段进行的
        if numa(i,1)>numa(i-1,1)&numa(i,1)>numa(i+1,1)
            peak(k,1)=numa(i,1);
            peak(k,2)=i;
            k=k+1;
        end;
    end;
    p=m.*0.4;%这里0.5是实验取得的，实验使用了三个数据，结果均为无误判点和漏判点
    count=0;
    ii=0;
    j=1;%用于控制周期计数
    for k=2:(Lv-L-1)
        if peak(k,1)>peak(k-1,1)&peak(k,1)>peak(k+1,1)&peak(k,1)>p
            i=peak(k,2);
            text(i-L-1,numa(i,1),' \leftarrow TOP POINT','FontSize',12);            
            while numa(i,1)>0
                i=i-1;
            end;
            text(i-601,numa(i,1),' \leftarrow BEGIN POINT','FontSize',12);
            a=peak(k,2)-i;
            b=i-ii;
%            b/(a+b)
            if (b/(a+b))>=0.9
                time(j,1)=4.*a./100;
                j=j+1;
                count=count+1;
            else
                time(j,1)=(a+b)./100;
                j=j+1;
                count=count+1;
            end;
%            time(j,1)                
            ii=peak(k,2);
            i=peak(k,2);
        end;
    end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%显示结果
%%%%%%%结果显示中第一种得到的一个周期为人实际走两步的时间，第二种得到的为一步的时间。
%if time(1,1)<(2/3).*time(3,1)
%    time(1,1)=0;
%end;
%if time(2,1)<(2/3).*time(3,1)
%    time(2,1)=0;
%end;
sum=0;
k=0;
for j=1:1000
    if time(j,1)>0
        sum=sum+time(j,1);
        k=k+1;
    end;
end;
sum=sum/k;
for j=1:1000%从此第三个周期开始显示，避免前两个喜好的不规则对整体的影响（缺乏理论依据）
    if time(j,1)<(2/3).*sum
        time(j,1)=0;
    end;
    if time(j,1)>(4/3).*sum
        time(j,1)=0;
    end;
    if time(j,1)>0
        time(j,1)
    end;
end;

count
                                                                                                                         




%figure;
%set(gca.'FoneSize',14);
%plot(numa(2000:3000,1),'g-');
%xlabel('[sec]');
%ylabel('[3d]');
%title('disposaled signal');


%%%point 2000 to point 3000
%figure;
%set(gca,'FontSize',14);
%plot(numd(:,1),'r-');
%xlabel('[sec]');
%ylabel('[3D]');
%title('3D acceleration');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fclose all;

return