图片上的红线并不是心率-时间曲线，而是处理后的心电数据，用于识别心跳

识别率接近99.8%

BPF是9-11Hz带通滤波器
信号取自自己，睡觉时记录了10分钟，8位量化，没有进行数字滤波


程序代码：
%pan-tompkins算法与最大值搜索联用

%--------------------------需要配置的参数-----------------------------------

%最小阈值
thresh_min = 4;
%滑动平均窗口宽度
windowsize = 300;
%进一步搜索最大值的范围
findmax_l_edge=-400;
findmax_r_edge=-100;
%新阈值计算系数
%新阈值 = thresh_refresh_k*上一个峰值;
thresh_refresh_k = 0.4;

%--------------------------运行前要输入的数据-------------------------------

%原始心电数据  要校正到2uV/LSB
%ecg = ecg_readcsv('holter.csv',1);
%采样率
fs = 1024;
%需要处理的数据范围
x = ecg.value(1:fs*200)*500;
%x = ecg.value*500;
%滤波器对象及延时
%load BPF.mat
%filter_delay = 200;
filter_delay = 0;

%---------------------------程序-------------------------------------------
disp('Stage 1 Filtering');

disp('    Step 1/5');
%x->原始数据,y->处理后数据
y = x;
disp('          OK');

disp('    Step 2/5');
%滤波
%x = filter(HPF,x);
%x = filter(LPF,x);
y = filter(BPF,y);
disp('          OK');

disp('    Step 3/5');
%微分
y = diff(y);
disp('          OK');

disp('    Step 4/5');
%平方每一项
for i=1:length(y)
    y(i) = y(i) * y(i);
   %if(y(i)<0)
   %    y(i)=-y(i);
   %end
end
disp('          OK');

disp('    Step 5/5');
%滑动平均
y = filter(ones(1,windowsize)/windowsize,1,y);
disp('          OK');

%------------------------Pan-Tompkins算法----------------------------------
peak_value=0;%峰值大小
peak_addr=0;%峰值位置

thresh = thresh_min;
thresh_updated = 0;

last_r_pos = 0;

clear rwave_stack;
rwave_stack = zeros(1,ceil(length(y)/(fs*60/300)));
rwave_stack_pos = 1;

disp('Stage 2 Pan-Tompkins');
last_percent = 0;
for i=1:length(y)
    %报告进度
    percent = i/length(y)*100;
    if( last_percent +1 == floor(percent))
         disp(strcat(num2str(percent),'%'));
         last_percent = floor(percent);
    end
    
    %当fs*2个样本内（心率最低30）有一个有效的峰值，就保存下来
    %查找峰值的范围
    if( last_r_pos == 0 )
        %之前没找到R波，从i处开始找峰值
        peak_search_l = i;
    else
        %之前找到了R波，在下一个峰值处寻找
        %心率最多300，间隔时间=去掉小数(fs*60/300)
        peak_search_l = peak_addr + floor(fs*60/300);
    end
   
    %搜索完毕
    if(peak_search_l > ( length(y) -40 ))
        break;
    end
   
    %15秒找不到峰值，复位所有参数
    peak_search_r = peak_search_l + 15*fs;
   
    if(peak_search_r > length(y) - 40)
        peak_search_r = length(y) - 40;
    end
   
    if(peak_search_r < peak_search_l)
        peak_search_r = peak_search_l;
    end
       
    for j = peak_search_l:peak_search_r
        %抓住一个极大值
        if( y(j)<y(j+10) && y(j+10)<y(j+20) && y(j+20)>y(j+30) && y(j+30)>y(j+40) )
            %精确定位极大值
            peak_value = max( y(j:j+40) );
            for k = j:j+40
                if(y(k)==peak_value)
                    peak_addr = k;
                end
            end
            %如果100ms后（R波宽度的一半）有更大的，则丢弃现在这个
            %先不写
           
            %有效峰值？
            if( peak_value >= thresh )
                %保存这个R波
                rwave_stack(rwave_stack_pos) = peak_addr - filter_delay ;
                rwave_stack_pos = rwave_stack_pos +1;    
                last_r_pos = peak_addr - filter_delay ;
                %刷新阈值
                thresh = thresh_refresh_k*peak_value;
                if(thresh < thresh_min )
                    thresh = thresh_min;
                end
                thresh_updated =1;
                break;
            else
                thresh_updated =0;
            end
        end
    end
   
    %没找到R波
    if(thresh_updated == 0 )
        %复位检测参数
        peak=0;
        peak_addr=0;
        thresh=thresh_min;
        thresh_updated=0;
    end
   
end

%作图
figure;
%plotyy(x,'b',y,'r');
plot(x,'b');

%----------------------在原始信号上搜索峰值，进一步定位R波-------------------
disp('Stage 3 Find peaks');
last_percent = 0;
for i=1:rwave_stack_pos-1
    %报告进度
    percent = i/(rwave_stack_pos-1)*100;
    if( percent == floor(percent))
         disp(strcat(num2str(percent),'%'));
         last_percent = floor(percent);
    end
    peak_value = max( x(rwave_stack(i)+findmax_l_edge:rwave_stack(i)+findmax_r_edge) );
    for j=rwave_stack(i)+findmax_l_edge:rwave_stack(i)+findmax_r_edge
        if( peak_value == x(j))
            rwave_stack(i) = j;
            break;
        end
    end
end

%-----------------------在图上标出R波位置-----------------------------------
disp('Marking Rwave');
for i=1:rwave_stack_pos-1
hold on;
plot(rwave_stack(i),x(rwave_stack(i)),'go');
end

%--------------------------统计数据----------------------------------------
rr_interval = diff(rwave_stack);
disp('Stage 4 Statistic');
bpm_average = 0;
for i=1:rwave_stack_pos-2
     bpm_average = bpm_average + rr_interval(i);
end
bpm_average = bpm_average/(rwave_stack_pos-2);
bpm_average = bpm_average * 60 /fs;
bpm_max = max(rr_interval(1:rwave_stack_pos-2)) * 60 /fs;
bpm_min = min(rr_interval(1:rwave_stack_pos-2)) * 60 /fs;

%心率-时间图
hold on;
plot( (1: length(x)/(rwave_stack_pos-2): length(x)), rr_interval(1:rwave_stack_pos-2)*60/fs, 'r');
disp('The red line is bpm vs time');

%-----------------------------报告数据-------------------------------------
disp(strcat('Max beats per second:',num2str(bpm_max)));
disp(strcat('Min beats per second:',num2str(bpm_min)));
disp(strcat('Average beats per second:',num2str(bpm_average)));