tic
clear;
clc;

fp = fopen('G:\gzb現(xiàn)地單元數(shù)據(jù)\PD071227094503109.dat','r');
count = 600000;
threshold1 = 0.05;
threshold2 = 2.48e-3;
threshold3 = 0.0005;
N = 10;

status = fseek(fp,offset,'bof');
raw_data = fread(fp,count,'int16');
fclose(fp);
float_data = (raw_data - 2048)*10.0/4096.0;
fs = 30000000.0;
f_filter = 100000.0;
[b,a] = cheby1(7,0.5,f_filter/fs,'high');
filter_data = filter(b,a,float_data);
pulse_data = zeros(1,length(raw_data));
j = 0;
begin_post = 0;
end_post = 0;
last_end_post = 0;
last_begin_post = 0;
pulse_infor = zeros(100,2);
pulse_num = 1;
i = 1000;
while i<=length(raw_data)
    if abs(filter_data(i))> threshold1
      begin_post = 0;
      end_post = 0;
      for j = i:length(raw_data)
         temp_power = 0.0; 
         c = 0.0;
         c = filter_data(j-2) + filter_data(j+2) - 2*filter_data(j);
         for k = j-N:j+N
             if j+N>length(raw_data);
                 break;
             else
            temp_power = temp_power + filter_data(k)*filter_data(k);
             end
         end
         if temp_power > threshold2&&(abs(c)>threshold3)
             continue;
         else
             end_post = j;
             begin_post = i;
             if begin_post>last_end_post+256
             pulse_data(begin_post:end_post) = filter_data(begin_post:end_post);
             pulse_infor(pulse_num,:) = [begin_post,end_post];
             pulse_num = pulse_num + 1;
             last_begin_post = begin_post;
             last_end_post = end_post;
             else
             pulse_data(last_begin_post:end_post) = filter_data(last_begin_post:end_post);
             pulse_infor(pulse_num-1,:) = [last_begin_post,end_post];
             last_end_post = end_post;
             end
             i = j + 1;
             break;    
         end
      end
    else
        i = i + 1;
    end
end
all_value = zeros(pulse_num-1,256);
for i=1:pulse_num-1
   array_pulse = filter_data(pulse_infor(i,1):pulse_infor(i,2));
   if length(array_pulse)>256
   array_pulse = array_pulse(1:256); 
   else
   temp = zeros(256 - length(array_pulse),1);
   array_pulse(1:256) = [array_pulse;temp];   
   end   
   max_value = max(abs(array_pulse));
   array_pulse = array_pulse/max_value;
   all_value(i,:) = array_pulse'; 
end

start = pulse_infor(1:pulse_num-1,1);
width = pulse_infor(1:pulse_num-1,2)-pulse_infor(1:pulse_num-1,1);
threshold = 0.85;
indication(1:pulse_num-1) = 0;
m = 1;
n = 2;
cluster_no = 1;
cluster = [];
while m<pulse_num
    no = 1;
    if isequal(indication(m),1)
        m = m+1;
        continue
    else
        plusem = all_value(m,:);
%         plusen = plusem;
        cluster(cluster_no,no) = m;        
        no = no+1;
         indication(m) = 1;
         n = m + 1;
    end
    while n<pulse_num
        if isequal(indication(n),1)
            n = n+1;
            continue;
        else
            plusen = all_value(n,:);
        end

%          ffm = fft(plusem,256);
%          ffn = fft(plusen,256);
%          ffm = ffm(1:128);
%          ffn = ffn(1:128);
%          fm = ffm.*conj(ffm);
%          fn = ffn.*conj(ffn);
%          for i=1:128;
%             fm(i) = sqrt(fm(i));
%          end
%          for i=1:128;
%             fn(i) = sqrt(fn(i));
%          end
               ppp = 128;
         pluset = plusem(1:2:255);
         pluseq = plusen(1:2:255);
         ffm = fft(pluset,ppp);
         ffn = fft(pluseq,ppp);
         ffm = ffm(1:(ppp/2));
         ffn = ffn(1:(ppp/2));
         fm = ffm.*conj(ffm);
         fn = ffn.*conj(ffn);
         for i=1:(ppp/2);
            fm(i) = sqrt(fm(i));
         end
         for i=1:(ppp/2);
            fn(i) = sqrt(fn(i));
         end
        cc = fm*fn'/sqrt(fm*fm'*fn*fn')
        if cc>threshold
            cluster(cluster_no,no) = n;
            no = no+1;
            indication(n) = 1;
            n = n+1;
        else
            n = n+1;
        end
    end  
    cluster_no = cluster_no+1;
end
figure(1)
subplot(311)
plot(float_data);
subplot(312)
plot(filter_data);
subplot(313)
plot(pulse_data);
% for i = 1:cluster_no
%     for j = 1:no
%         cluster_data
%     end
% end
figure(2)

% cluster_data(i,:) = 
[p,q] = size(cluster);
cluster_data = zeros(p,length(raw_data));
for i=1:p;
%     cluster_data = zeros(1,length(raw_data));
    for j = 1:q
       if cluster(i,j)~=0
           start = pulse_infor(cluster(i,j),1);
           over  = pulse_infor(cluster(i,j),2);
           cluster_data(i,start:over) = pulse_data(start:over);
       else
%            start = pulse_infor(1:pulse_num-1,1);
%            over = pulse_infor(1:pulse_num-1,2);
           break;
       end
           
    end
    subplot(p,1,i)
    plot(cluster_data(i,:));
end
% figure(3)
% plot(cluster_data(2,:) + cluster_data(4,:));
toc