read txt

clc;
clear;

%ekg_raw=load('2.txt');
ekg_raw=load('20210110133815000_20210110134234000.txt');
%x=ekg_raw(:,1);
%y1=ekg_raw(:,2);




[row,col]=size(ekg_raw);

ekg=ekg_raw(1,:);
for i=2:row
    if ekg_raw(i,1) == ekg_raw(i-1,1);
        %break
        continue
    else
        ekg=[ekg;ekg_raw(i,:)];
    end
end

%ekg=ekg_raw;


%读取电压数据
fs=500;%分析频率 /Hz
%ecg=interp1(ekg(:,1),ekg(:,2),ekg(1,1):1/fs:ekg(end,1))';       


X1=ekg(:,1);
Y1=ekg(:,2);
S1=ekg(1,1);
E1=ekg(end,1);

P1=S1:1/fs:E1;
P2=fs:E1;


%读取电压数据
%ecg=interp1(ekg(:,1),ekg(:,2),ekg(1,1):1/fs:ekg(end,1))';       
ecg=interp1(X1,Y1,P1)';



%-----------------------------------------------------%
% 利用pan_tomkin算法找到R点
[map,r,delay]=pan_tompkin2021(ecg,fs,1);
    map=map';
    r=r';
    l=length(r);
    
for i=2:l;
    t(i-1,1)=(r(i)-r(i-1))/fs*1000; %求出R-R间的时间值,单位为ms
end

x=r(2:l);%  index of R waves
y=interp1(x,t,r(2):1:r(l),'spline'); %利用插值法求出以原ecg信号的采样率fs的拟合函数

figure %整个时间段内R-R时间差
plot([r(2):1:r(l)]/fs,y')
xlabel('时间/s');ylabel('R-R时间差/ms')
hold on
scatter(x/fs,t(1:l-1));

rawdata=y';
    [row,col]=size(rawdata);
    Nfft=2^(nextpow2(length(rawdata)));  % FFT数量

data=rawdata-mean(rawdata);%去直流电平
    data_fft=fft(data,Nfft);
    mag=abs(data_fft);
Pxx=mag.^2/Nfft/fs;% 功率谱密度
    f=(0:Nfft/2)'*fs/Nfft;%频率轴

figure
f_cut=0.5;%绘图截止频率/Hz，超过0.5的基本为0
plot(f(1:f_cut*fs),Pxx(1:f_cut*fs));
    xlabel('频率/Hz');
    ylabel('PSD/(ms^2/Hz)');

disp ('时域分析')
RR=mean(t);
SDNN=sqrt(var(t));
rMSSD=rms(diff(t));
    disp(['RR=',num2str(RR),' ms']);
    disp(['SDNN=',num2str(SDNN),' ms']);
    disp(['rMSSD=',num2str(rMSSD),' ms']);

disp('频域分析')
f1=0.4;%/Hz
f2=0.15;%/Hz
f3=0.04;%/Hz
TP=fs/Nfft*trapz(Pxx(1:floor(f1/(fs/Nfft))));%/ms^2
HF=fs/Nfft*trapz(Pxx(floor(f2/(fs/Nfft)):floor(f1/(fs/Nfft))));%/ms^2
LF=fs/Nfft*trapz(Pxx(floor(f3/(fs/Nfft)):floor(f2/(fs/Nfft))));%/ms^2
    disp(['TP=',num2str(TP),' ms^2']);
    disp(['HF=',num2str(HF),' ms^2']);
    disp(['LF=',num2str(LF),' ms^2']);
    disp(['LF/HF=',num2str(LF/HF)]);