fp = fopen('G:\gzb現地單元數據\PD071213101631125.dat','r');
count = 600000;
threshold1 = 0.8; %0.05
threshold2 = 0.0248;%2.48e-4;
threshold3 = -0.5;%-0.05;
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;
% for i=1000:length(raw_data)
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; 
         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
             continue;
         else
%              pulse_data(i:j) = filter_data(i:j);
             end_post = j;
             begin_post = i;
             if begin_post>last_end_post+512
             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
%       for j = begin_post:-1:512
%          temp_power = 0.0; 
%          for k = j-N:j+N
%             temp_power = temp_power + filter_data(k)*filter_data(k);
%          end
%          if temp_power > threshold2
%              continue;
%          else
% %              pulse_data(i:j) = filter_data(i:j);
%              begin_post = j;
%              break;    
%          end
%       end
%       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;
    else
        i = i + 1;
    end
end

all_value = zeros(pulse_num-1,512);
for i=1:pulse_num-1
   array_pulse = filter_data(pulse_infor(i,1):pulse_infor(i,2));
   if length(array_pulse)>512
   array_pulse = array_pulse(1:512); 
   else
   temp = zeros(512 - length(array_pulse),1);
   array_pulse(1:512) = [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 = 1;
cluster_no = 20;
cluster = zeros(cluster_no,pulse_num);
while m<pulse_num-1
    if isequal(indication(m),1)
        m = m+1;
    else
        plusem = all_value(m,:);
        cluster(1,1) = m;
%         indication(m) = 1;
    end
    while n<pulse_num-1
        if isequal(indication(n),1)
            n = n+1;
        else
            plusen = all_value(n,:);
        end
         ffm = fft(plusem,512);
         ffn = fft(plusen,512);
         ffm = ffm(1:256);
         ffn = ffn(1:256);
         fm = ffm.*conj(ffm);
         fn = ffn.*conj(ffn);
         for i=1:256;
            fm(i) = sqrt(fm(i));
         end
         for i=1:256;
            fn(i) = sqrt(fn(i));
         end
        cc = fm*fn'/sqrt(fm*fm'*fn*fn')
        if cc>threshold
            cluster(1) = m;
            cluster_no(2) = n;
        else
            n = n+1;
        end
    end
    
end

% xx = all_value(1,:);
% yy = all_value(5,:);
% relation0 = abs(xx)*abs(yy')/sqrt(xx*xx'*yy*yy')
% xx = xx - mean(xx);
% yy = yy - mean(yy);
% relation = abs(xx)*abs(yy')/sqrt(xx*xx'*yy*yy')
% relation = all_value(5,:)*all_value(3,:)'/sqrt(all_value(5,:)*all_value(5,:)'*all_value(3,:)*all_value(3,:)')
% [m,n] = size(all_value);
clc
for kkk = 2:7
 xx = all_value(kkk,:);
 yy = all_value(1,:);
 ffxx = fft(xx,512);
 ffyy = fft(yy,512);
 ffxx = ffxx(1:256);
 ffyy = ffyy(1:256);
 fxx = ffxx.*conj(ffxx);
 fyy = ffyy.*conj(ffyy);
 

 for i=1:256;
     fxx(i) = sqrt(fxx(i));
 end
 for i=1:256;
     fyy(i) = sqrt(fyy(i));
 end
 rel = fxx*fyy'/sqrt(fxx*fxx'*fyy*fyy')
end
 figure(1)
subplot(311)
plot(float_data);
subplot(312)
plot(filter_data);
subplot(313)
plot(pulse_data);
figure(2)
for i=1:6;
    subplot(3,2,i)
    plot(all_value(i,:));
end
figure(3)
xx = all_value(3,:);
 ffxx = fft(xx,512);

 fxx = ffxx.*conj(ffxx);

 

 for i=1:512;
fxx(i) = sqrt(fxx(i));
 end
plot(fxx)