yxok.m

读入一段数据后

clear;
close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%            相关参数
fln=102;%文件名
channel1=1;%道数
channel2=4;


Nf=4096*0+1;     %起始点
Nl=4096*6;    %终点
pagesn=4096*4;  %假定为一次判别的点数，称其为一页


sample=1025;      %采样率
timesn=500;      %持时的点数


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 读文件
for chani=channel1:channel2                 %chani为道的变量
    flnm1='d:\data\gm2\c';
    flnm2='#';
    flc=chani;                          %将道变为文件名
    flnm3='.STS';
    filename=strcat(flnm1,num2str(fln),flnm2,num2str(flc),flnm3);
    fid=fopen(filename,'r');
    data=fread(fid,inf,'ushort');           
    dataok(chani,:)=data(Nf:Nl)/2^16*10-5;          %dataok读取了Nf到Nl个点
    
    figure(1);
    subplot((channel2-channel1)+1,1,chani);
    plot(dataok(chani,:));        
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%单个事件判别
%假定取得的40960个数据是由采集到的10页数据拼合而成；
%计算背景值--每页4096个点；



for chani=channel1:channel2
    mysum(chani)=0;
    for sn=1:pagesn
        mysum(chani)=mysum(chani)+abs(dataok(chani,sn));
    end
    noisevalue(chani)=mysum(chani)/pagesn;  %获得背景噪声值；  
end

%事件阀值设定为背景值的10倍，如果发现阀值符合时，间隔250个点，
for chani=channel1:channel2    
    eventnum(chani)=0;  %eventsn事件个数
    eventsn=1;          %进行事件判别的时间点数，本质对应于时间
    while(eventsn<pagesn)
        if(abs(dataok(chani,eventsn))>noisevalue(chani)*8)     %如果数值绝对值大于背景噪声值的10倍       
                eventtime(chani,eventnum(chani)+1)=eventsn;  %eventsn的值就是事件的位置，即时间点
                eventnum(chani)=eventnum(chani)+1;   %事件序号的值加1
                eventsn=eventsn+timesn;                          %样点加250，因为认为持时为250
        else
            eventsn=eventsn+1;                             %样点加1
        end         
    end
%    eventnum(chani)=eventnum(chani)-1;  %该道事件个数      
end
eventtime
%事件持时设定为250个点，则能量计算时总点数为250点。并设定能量阀值为3  

for chani=channel1:channel2
    if eventnum(chani)>0              %这里需要假定利用单点判定的事件个数大于1
        eventnlnum=0;
        eventok(chani)=0;
        for eventnumnl=1:eventnum(chani)  %单点事件个数    
            sumnl(chani,eventnumnl)=0;
            
            for snnl=eventtime(chani,eventnumnl)-timesn/10:eventtime(chani,eventnumnl)+timesn/10*9  %每个事件从阀值处前50点到后200点
                if snnl<Nl&snnl>0             %这个地方需要补充
                    sumnl(chani,eventnumnl)=sumnl(chani,eventnumnl)+abs(dataok(chani,snnl));      %250个点相加
                end
            end
            if(sumnl(chani,eventnumnl)>noisevalue(chani)*timesn*1.5)
                    %eventtimenl(chani,eventnumnl)=eventtime(chani,eventnumnl);
                    eventnlnum=eventnlnum+1;
                    eventnl(chani,eventnlnum)=eventtime(chani,eventnumnl); 
                    eventok(chani)=eventok(chani)+1;
                    
                    
            else 
                eventnl(chani,eventnumnl)=1;
            end
        end
    else
    end
end
eventnl
eventok(chani);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%对每个墓进行多道报警判别
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 %       mu1=[1 2 3 4];    
 %       mu2=[4 5 6 7];      
 %       mu3=[8 9 10 11];    
 %       mu4=[11 12 13 14];  
 %       mu5=[15 16 17 18];  
 %       mu6=[18 19 20 21];  
 %       mu7=[22 23 24 25];  
 %       mu8=[25 26 27 28];  
 %       mu9=[29 30 ];              
%

mualart=[1 4 8 11 15 18 22 25 29];  %墓的检波器编号
musg=3; %每次判别至少需要的事件次数
MualartNum=2;   %判定一个墓报警需要的最少道数
MuOk(1:9)=0;    %每个墓的判定结果
for MuNum=1:1
    mugroup(MuNum)=0;    
   for MuChan=0:3
%       mualart=strcat('mu',str2num(MuNum));
 %      if length(eventnl(mualart(MuNum)+MuChan,:))>musg
       if eventok(mualart(MuNum)+MuChan)>=musg           
           mugroup(MuNum)=mugroup(MuNum)+1;
       end
   end
   if mugroup(MuNum)>=MualartNum
       MuOk(MuNum)=1;
   end
end