ECG基线漂移.zip

clear all; clc; A = xlsread('C:\Users\mxy\Desktop\波形数据.xlsx');%提取第一列元素data=load('E:\signal\1\1_1tang.txt','r') channelOne = A(1:2:end,1);%取第1列，奇数行的数据 channelTwo = A(2:2:end,1);%双通道，1通道为ECG信号，2通道为PPG信号 %×××××××××××××××ECG_PPG信号预处理×××××××××××××××% %***********************************************% %************ECG消除基线漂移*************% %[b,a]=cheby1(n,Rp,W'n,ftype')其中n为滤波器阶数；Rp为通频带阶数；ftype为可以设置的滤波器类型； %b、a分别是滤波器分子和分母的多项式系数.求得n=4，Rp=0.5，Wn=0.6/500(归一化截止频率),通带截止频率f1=0.3，阻带截止频率f2=0.5 %[channelOne,count]=fscanf(fid,'%f%f',18000);%读取向量的前18000个数据 %fclose(fid); fs=500; t=0:1/500:(36-1/500);%渐进，从0到（36-1/500），步幅1/500（采样频率500Hz），共18000个元素 %t=(1:fs)/fs=0:1/500:(36-1/500);%时间序列（时间轴）T = 1/Fs；L=18000；%总的采样点数；t = (0:L-1)*T L=18000;%信号长度，总的采样点数 x=zeros(1,18000); for i=1:L x(i)=channelOne(i); end x=smooth(x,30,'lowess'); fig=figure(1); fig.Position=get(0,'ScreenSize'); subplot(2,1,1)%把图形窗口分割成2行2列，t-x在第一行 plot(t,x);%以t为横坐标值，x为纵坐标值作图 zoom on xlabel('time');%横坐标标注：时间 ylabel('magntitude');%纵坐标标注：量级 title('(a)原始的心电信号'); grid on %打开坐标网格线 for i=1:length(x) x(i)=x(i)-1628;%将信号平移 end %----------------------------------------------------------% % for k=1:L-49; % SSF=cumsum(y,1,'reverse');%从后往前累加 % for kk=1:L-1;%L-1=k+50-2=j % SSF(k)=cumsum(y(kk));%窗宽50 % end % flag = SSF(k)==0 && SSF(k+1)>0; indmax=find(diff(sign(diff(x)))<0)+1;%极大值点 indmin=find(diff(sign(diff(x)))>0)+1;%极小值点 subplot(2,1,2); plot(t,x,t(indmin),x(indmin),'r*')%,t(indmax),x(indmax),'ro' title('极值点图');legend('ECG信号','极小值点');%加图例 grid on % if flag % else x(i)=0; % end % end %----------------------------------------------------------% %获取需要的时域特征点--------------------------------------% b=length(indmin);%设定b为极小值点的个数 B=ones(b,2);%创建一个新矩阵B，有b行2列 B(:,1)=t(indmin);%第1列为特征点的横坐标 B(:,2)=x(indmin);%第2列为特征点的纵坐标 BX=B(:,1);BY=B(:,2); D=0.25;%窗宽60*1/500HZ=0.12，也令其为阈值 d_ecg = diff(BX); % 计算两两极小值点之间横轴距离 j=1:400; ecg_x=ones();%预设一个数组,准备存放符合条件的特征点的横坐标 ecg_y=ones();%纵坐标 for i=1:length(d_ecg) if abs(d_ecg(i)) > D%距离的绝对值大于窗宽 ecg_x(j)=BX(i);%把符合条件的左端点放入big数组第1位 ecg_y(j)=BY(i); ecg_x(j+1)=BX(i+1);%把符合条件的右端点放入big数组第2位 ecg_y(j+1)=BY(i+1); j=j+1; %完成一次计算后进行自加，跳入下一数值进行运算 elseif BY(i)>BY(i+1) %如果距离小于窗宽 %当左端点的信号值>右端点信号值时 ecg_x(j)=BX(i+1);%取右端点存入数组 ecg_y(j)=BY(i+1); j=j+1; else ecg_x(j)=BX(i);%否则，取左端点存入数组 ecg_y(j)=BY(i); j=j+1; end end %这种方法得到的数组肯定会出现相邻的相同元素 ecg_x(diff(ecg_x)==0)=[];%去除掉相邻两个差值为0的点其中之一 ecg_y(diff(ecg_y)==0)=[]; %disp(big_x);%显示数组 %figure(2);plot(big_x,big_y); %对特征点做插值------------------------------------------% %三次多项式插值法 e_inter=interp1(t,x,ecg_x,'cubic','extrap');%对ecg_x中超出已知点集的插值点用指定插值方法计算函数值 figure(2);subplot(211); plot(t,x,ecg_x,e_inter,'r'); title('三次多项式插值');%即立方插值 %legend('a','b');%加图例 hold on %用插值结果对原信号进行拟合-------------------------------% c=length(ecg_x);%设定c为符合条件的极小值点的个数 C=ones(c,2);%创建一个新矩阵C，有c行2列 C(:,1)=ecg_x;%第1列为特征点的横坐标 C(:,2)=ecg_y;%第2列为特征点的纵坐标 %big_x=big_x(1:end);big_y=big_y(1:end);%将数组转化为向量 %figure(3) % plot(t,x,big_x,inter,'b.','markersize',15); % hold on; e10=polyfit(ecg_x,e_inter,10);%12次多项式拟合 e_poly=polyval(e10,t);%求多项式在t处的值 plot(t,e_poly,'k-');grid on legend('原ECG信号','特征点立方插值','拟合曲线');%加图例 e_fit=x'-e_poly;subplot(212) plot(t,x,t,e_fit);legend('原始ECG信号','插值拟合减基线');%加图例 title('减基线漂移后的ECG信号'); figure(3); hua_fft(x,fs,1); grid on hold on hua_fft(e_fit,fs,1); title('插值拟合前后的频域图');grid on axis([-0.5 3 0 0.5*10^6]);legend('插值拟合前','插值拟合后') %*************************************************************% %************PPG消除基线漂移**********************************% % fs=500; % t=0:1/500:(36-1/500); % for j=1:18000 % x(j)=channelTwo(j); % end %************构造拟合函数*************% %*************选取拟合点**************% %**根据数字滤波器性能指标设计巴特沃斯模拟滤波器滤除工频噪声**% %******借助设计好的模拟滤波器设计巴特沃斯数字滤波器*********% %******利用设计好的数字滤波器对ECG信号进行滤波处理**********% %××××××××××××××××××××××××××××××××××××% %×××××××××××根据条件筛选出符合要求的ECG信号××××××××××% %××××××××××××××××××××××××××××××××××××% %×××××××××××××××ECG信号特征提取×××××××××××××××% %*************小波变换提取ECG特征点**************% %×××××××××××××××PPG信号特征提取×××××××××××××××% %*************差分法提取波峰数**************% %*计算主波峰值后20个点的两点间斜率平均值得到主波斜率*% %*************20%主波宽度**************% %*************傅立叶变换获得脉率**************% %*************局部搜索极大值法获得A2/A1、A3/A1**************% %×××××××××××××××特征提取信号预处理×××××××××××××××% %******************归一化处理****************% %*************脉搏波信号主成份分析**************% %×××××××××构建自适应神经网络模糊推理系统××××××% %×××××××××××ANFIS检测绩效分析××××××××××%