tic
clear;
clc;

fp = fopen('G:\gzb現(xiàn)地單元數(shù)據(jù)\PD071206180738453.dat','r');
count = 600000;
threshold1 = 0.05; 
threshold2 = 0.0006;
threshold3 = 0.05;
N = 20;
raw_data = fread(fp,count,'int16');
fclose(fp);
float_data = (raw_data - 2048)*10.0/4096.0;
fs = 30000000.0;
f_filter = 100000.0;
pulse_data = zeros(1,length(raw_data));
[b,a] = cheby1(6,0.5,f_filter/fs,'high');
filter_data = filter(b,a,float_data);

subplot(311)
plot(float_data);
subplot(312)
plot(filter_data);
subplot(313)

j = 0;
i0 = 0;
for i=20:length(raw_data)
    if abs(filter_data(i))> threshold1
       i0 = mod(i,10);
       i0 = i + 10 - i0;
       for j = i0:10:length(raw_data)
         
         c = 0.0;
         temp_power = 0.0; 
         c = filter_data(j-2) + filter_data(j+2) - 2*filter_data(j);
         for k = j-N:j+N
            if abs(filter_data(k))> 0.2*threshold1
              temp_power =  temp_power + (abs(filter_data(k)) - 0.4*threshold1)*(abs(filter_data(k)) - 0.4*threshold1);
            end
         end
         
          if (temp_power>threshold2)&&(abs(c)<threshold3)  
%          if (temp_power>threshold2)
             continue;
         else
%            if j>i+8 
%              c
%              temp_power
             pulse_data(i:j) = filter_data(i:j);
%            end
           i = j + 1;
           break;    
         end
      end    
    end
end

plot(pulse_data);
toc