基于MATLAB实现的对心电信号进行滤波处理后，求得瞬时心率和心率变异系数，求得各个自律神经活动的评价指标+使用说明文档

clear all;clc %%%%%%%%%% read data %%%%%%%%%% load('F:/MP100/2min-3/9-21/kokatu 20120921 ansei arai ansei') ecg=data(:,1);ecg=ecg'; %% 'ecg' means 'ECG' samp=2000; %% 'samp' means 'sampling frequency' %%%%%%%%%% calculate R-R interval / heart rate / heart rate variability %%%%%%%%%% b=fir1(50,[10/(samp/2),50/(samp/2)]); %% design a 50th-order FIR bandpass filter with passband 10/(samp/2)≒肋≒50/(samp/2) ecg_filt=filtfilt(b,1,ecg); %% zero-phase filtering; 'b' provides numerator coefficients of the filter,'1'provides the denominator coefficients of the filter time_sec=(1:length(ecg_filt))/samp; time_min=time_sec/60; figure subplot(3,1,1);%% %% check the difference between original and filtered resipiration signals plot(time_sec,ecg,'b',time_sec,ecg_filt,'r','linewidth',2);title('Comparison of Different ECG Signals');xlabel('Time (s)');ylabel('Potential (V)');axis([60,70,-1,2]);legend('Original','Filtered');grid; subplot(3,1,2); plot(time_sec,ecg,'b',time_sec,ecg_filt,'r','linewidth',2);xlabel('Time (s)');ylabel('Potential (V)');axis([180,190,-1,2]);legend('Original','Filtered');grid; subplot(3,1,3); plot(time_sec,ecg,'b',time_sec,ecg_filt,'r','linewidth',2);xlabel('Time (s)');ylabel('Potential (V)');axis([300,310,-1,2]);legend('Original','Filtered');grid; ecg_filt_diff_ave=mean(abs(diff(ecg_filt))); find_r1=find(ecg_filt>ecg_filt_diff_ave*10); %% 10 is thought to be standard,10 can be changed to other values find_r2=find(ecg_filt<mean(ecg_filt(find_r1))*0.7); %% 0.7 is though to be standard, 0.7 can be changed to other values ecg_filt(find_r2)=0; %% mark points not belong to R-wave with 0 ecg_filt_abv0=find(ecg_filt>0); %% find the start points of R-wave % ecg_filt_abv0(1)=[]; r_start1=ecg_filt(ecg_filt_abv0).*ecg_filt(ecg_filt_abv0-1); r_start2=find(r_start1==0); r_start3=ecg_filt_abv0(r_start2); for i=1:length(r_start3) %% find the peak points of R-wave [r_max,n]=max(ecg_filt(r_start3(i):r_start3(i)+30)); %% 'r_max' is the maximum value; 'n' is the sequence number of the maximun value r_point1(i,1)=r_start3(i)+n-1; r_time1(i,1)=r_point1(i,1)/samp; r_value1(i,1)=r_max; end r_value_ave=mean(r_value1)*2/3; r_time=[]; r_value=[]; for i=1:length(r_value1); if r_value1(i)>r_value_ave r_time=[r_time;r_time1(i)]; r_value=[r_value;r_value1(i)]; end end figure %% check the peak points of R-wave subplot(6,1,1); plot(time_min,ecg,'k',r_time/60,r_value,'r.');title('Peak Points of R-wave');ylabel('Potential (V)');axis([0,1,-1,2]);set(gca,'XTick',[0,1],'XTickLabel',{'0';'1'});grid; subplot(6,1,2); plot(time_min,ecg,'k',r_time/60,r_value,'r.');ylabel('Potential (V)');axis([1,2,-1,2]);set(gca,'XTick',[1,2],'XTickLabel',{'1';'2'});grid; subplot(6,1,3); plot(time_min,ecg,'k',r_time/60,r_value,'r.');ylabel('Potential (V)');axis([2,3,-1,2]);set(gca,'XTick',[2,3],'XTickLabel',{'2';'3'});grid; subplot(6,1,4); plot(time_min,ecg,'k',r_time/60,r_value,'r.');ylabel('Potential (V)');axis([3,4,-1,2]);set(gca,'XTick',[3,4],'XTickLabel',{'3';'4'});grid; subplot(6,1,5); plot(time_min,ecg,'k',r_time/60,r_value,'r.');ylabel('Potential (V)');axis([4,5,-1,2]);set(gca,'XTick',[4,5],'XTickLabel',{'4';'5'});grid; subplot(6,1,6); plot(time_min,ecg,'k',r_time/60,r_value,'r.');xlabel('Time (min)');ylabel('Potential (V)');axis([5,6,-1,2]);set(gca,'XTick',[5,6],'XTickLabel',{'5';'6'});grid; start_con=0:120:240; rr_inter=diff(r_time); rr_inter_ave=mean(rr_inter); r_time(1)=[]; %% assign the R-R interval to the latter R-wave r_value(1)=[]; if rr_inter_ave>1.0 %% define the maximum threshold (when HR=60, rr_inter_ave=1.0) thre_max=1.7; %% when HR<60, the maximun threshold=1.7 else thre_max=1.2; %% the standard maximun threshold=1.2 end rr_overmax=find(rr_inter>thre_max); error_t=r_time(rr_overmax); %% find the time that error occurs error_n=length(error_t); %% count the number of error points error_t=[-1;error_t]; rr_inter(rr_overmax)=[]; %% omit values larger than the maximum threshold r_time(rr_overmax)=[]; r_value(rr_overmax)=[]; hr=60./rr_inter; %% calculate the instant HR based on R-R interval figure plot(r_time/60,hr,'.-',error_t,40,'r*');title('Instant Heart Rate (*error)');set(gca,'XTick',[0:2:max(time_min)],'XTickLabel',{'0';'2';'4';'6'});xlabel('Time (min)');ylabel('Heart Rate (bpm)');axis([0,max(time_min),30,100]); grid; time_inter_min=0:60:max(time_sec); peak_num_min=histc(r_time,time_inter_min); time_inter_con=0:120:max(time_sec); peak_num_ave_con=histc(r_time,time_inter_con)/2; peak_num_min(end)=[]; peak_num_ave_con(end)=[]; hr_ave_min=[]; for i=0:max(time_min)-1 %% calculate the average HR of each minute hr_ave_min1=mean(hr(find(r_time>i*60 & r_time<=(i+1)*60))); hr_ave_min=[hr_ave_min;hr_ave_min1]; end hr_ave_con=[]; for i=1:length(start_con) %% calculate the average HR of each condition hr_ave_con1=mean(hr(find(r_time>(i-1)*120 & r_time<=i*120))); hr_ave_con=[hr_ave_con;hr_ave_con1]; end hr_ave_min2=[hr_ave_min,peak_num_min]; hr_ave_con2=[hr_ave_con,peak_num_ave_con]; start_min=0:60:300; figure subplot(1,2,1); bar(hr_ave_min2,0.5);title('Average Heart Rate of Each Minute');legend('Instant HR','Peak Nuber');set(gca,'XTickLabel',{'1st';'2nd';'3rd';'4th';'5th';'6th'},'ygrid','on');xlabel('Time (min)');ylabel('Average Heart Rate (bpm)');axis([0,length(start_min)+1,0,100]);hold on subplot(1,2,2); bar(hr_ave_con2,0.5);title('Average Heart Rate of Each Condition');legend('Instant HR','Peak Nuber');set(gca,'XTickLabel',{'Cond1';'Cond2';'Cond3'},'ygrid','on');xlabel('Condition');ylabel('Heart Rate (bpm)');axis([0,length(start_con)+1,0,100]);hold on cvrr_con=[]; for i=1:length(start_con) %% calculate the coefficient of variation of R-R interval rr_inter1=rr_inter(find(r_time>(i-1)*120 & r_time<=i*120)); rr_inter_detr=detrend(rr_inter1)+mean(rr_inter1); cvrr_con1=(std(rr_inter_detr)/mean(rr_inter_detr))*100; cvrr_con=[cvrr_con;cvrr_con1]; end figure bar(cvrr_con,0.5,'b');title('Coefficient of Variation of R-R Interval');set(gca,'XTickLabel',{'Cond1';'Cond2';'Cond3'},'ygrid','on');xlabel('Condition');ylabel('Coefficient of Variation (%)');axis([0,length(start_con)+1,0,10]); %%%%%%%%%% frequency analysis %%%%%%%%%% resamp=10; r_time_spl=[1/resamp:1/resamp:max(time_sec)]'; rr_inter_spl=spline(r_time,rr_inter,r_time_spl); % figure % plot(r_time/60,rr_inter,'r*-',r_time_spl/60,rr_inter_spl,'g.','linewidth',2);title('Interpolated R-R interval');xlabel('Time (min)');ylabel('R-R interval (s)');axis([0,max(time_min),0,1]);grid; for i=1:length(start_con); %% frequency analysis of each condition P0=zeros(513,1); %% generate a 501-by-1 matrix of zeros for j=0:10; %% frequency analysis every 11 points rr_inter_spl_con=rr_inter_spl((start_con(i)+j)*resamp+1:(start_con(i)+j+110)*resamp); %% pick out the desired data window=hamming(length(rr_inter_spl_con)); %% generate a symmetric Hamming window noverlap=[]; nfft=2^10; fs=resamp; [Pxx,f]=pwelch(rr_inter_spl_con,window,noverlap,nfft,fs); %% estimate the power spectral density Pxx of the input signal using Welch's method LF=Pxx(find(f>0.04 & f<=0.15)); %% find the low-frequency components HF=Pxx(find(f>0.15 & f<0.4)); %% find the high-frequency components st1(j+1,:)=trapz(LF)/trapz(HF); %% integration ratio for estimating (sympathetic activity) st2(j+1,:)=trapz(HF); %% integration for estimating (parasympathetic activity) st3(j+1,:)=trapz(HF)/(trapz(LF)+trapz(HF)); %% integration ratio for estimating (parasympathetic activity) P0=P0+Pxx; end str1(i,:)=mean(st1); %% average of 11 points str2(i,:)=mean(st2); str3(i,:)=mean(st3); psd(:,i)=P0/11; end s