tic
clear;
clc;

fp = fopen('G:\gzb現地單元數據\PD071206180736609.dat','r');
count = 600000;
threshold1 = 0.05;
threshold2 = 2.48e-4;
threshold3 = 0.05;
N = 10;
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)-0.4*threshold1)*(filter_data(k)-0.4*threshold1);
             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
figure(1)
subplot(311)
plot(float_data);
subplot(312)
plot(filter_data);
subplot(313)
plot(pulse_data);
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.80;
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;
    end
    while n<pulse_num
        if isequal(indication(n),1)
            n = n+1;
        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
        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(2)
% for i=1:6;
%     subplot(3,2,i)
%     plot(all_value(i,:));
% end

toc