心电信号的特征提取、分析与处理*

数据来源:MIT-BIH数据库(可从以下数据中任选两组进行实验)
给出4组不同病例的心电信号数据,分别命名为“100-2-3”,“105-2-3”,“109-2-3”,“111-2-3”,每组数据以“.mat”形式存储。(在文章最后面附带四组数据库的压缩包)
每组数据长度N=2048,采样率fs=360Hz(硬件采集心电信号时,每秒采集360个点。注意设计FIR,IIR时可能用到该采样率。). 即2048点对应时间约为5.69s()

##内容

(1)谱分析: 取两段心电信号数据(不同病例),对该数据进行频谱分析(幅度谱、相位谱、功率谱);
(2)相关分析:分别计算两段心电信号的均值、方差、自相关函数与互相关函数;分析两段信号的相干函数曲线
对于相关函数进行循环相关函数与线性相关函数的对比;
(3)数字滤波器设计:
取一段心电信号,添加白噪声。分别作出原始信号、加噪信号的图;作出原始信号、加噪信号的自相关曲线,分析估计心电信号的准周期;取一段心电信号,添加高频噪声(1k-2khz),如高频正弦信号,频率自己选择。结合(1)中得出的结论,即ECG的主要能量分布结果,设计数字滤波器(IIR或FIR),去除高频噪声。(也可直接设计数字滤波器去除基线漂移)要求:对滤波器的截止频率的设置要给出说明;
(4)维纳滤波器去除工频干扰:
取一段心电信号,添加频率为50Hz的高斯白噪声(工频干扰)。设计维纳滤波器,分析滤波结果。
##matlab代码实现

clear
Fs=360;N=2048;
load('109-2-3.mat');
x1=y;
load('111-2-3.mat');
x2=y;
figure
subplot(2,1,1);plot(x1);grid;title('原始心电信号x1.109-2-3.mat');
subplot(2,1,2);plot(x2);title('原始心电信号x2.109-2-3.mat');grid
%(1)谱分析
mf1=fft(x1,N);%进行频谱变换(傅里叶变换)
mf2=fft(x2,N);
mag1=abs(mf1);
mag2=abs(mf2);
f1=(0:length(mf1)-1)*Fs/length(mf1);  %进行频率变换
f2=(0:length(mf2)-1)*Fs/length(mf2);
figure
plot(f1,mag1);axis([0,370,0,190]);grid;      %画出频谱图
xlabel('频率(HZ)');ylabel('幅值');title('心电信号x1幅度谱');
figure
plot(f2,mag2);axis([0,370,0,190]);grid;      %画出频谱图
xlabel('频率(HZ)');ylabel('幅值');title('心电信号x2幅度图');
figure
plot(f1,angle(mf1)/pi);axis([0 360 -1.3 1.3]);grid;
xlabel('频率(HZ)');ylabel('相位');title('心电信号x1相位谱');
figure
plot(f1,angle(mf2)/pi);axis([0 360 -1.3 1.3]);grid;
xlabel('频率(HZ)');ylabel('相位');title('心电信号x2相位谱'); wn1=x1;
wn2=x2;
P1=10*log10(abs(fft(wn1).^2)/N);
P2=10*log10(abs(fft(wn2).^2)/N);
f1=(0:length(P1)-1)/length(P1);
f2=(0:length(P2)-1)/length(P2);
figure
plot(f1,P1);grid
xlabel('归一化频率');ylabel('功率(dB)');title('心电信号x1的功率谱');
figure
plot(f2,P2);grid
xlabel('归一化频率');ylabel('功率(dB)');title('心电信号x2的功率谱');
%(2)相关分析
avr1=mean(x1);avr2=mean(x2);
fprintf('信号x1的均值= %f\n',avr1);fprintf('信号x2的均值= %f\n',avr2);
var1=var(x1);var2=var(x2);
fprintf('信号x1的方差= %f\n',var1);fprintf('信号x2的方差= %f\n',var2);
rx1=xcorr(x1,'biased');
rx2=xcorr(x2,'biased');
rx1x2=xcorr(x1,x2);
figure
plot(rx1);grid;title('心电信号x1的自相关函数');
figure
plot(rx2);grid;title('心电信号x2的自相关函数');
figure
plot(rx1x2);grid;title('心电信号x1,x2的互相关函数');
%(3)数字滤波器设计
SNR=10*log(100/8); % 2%是能量比
x11=awgn(x1,SNR);
figure
subplot(2,1,1), plot(x1);title('原信号x1');% 加入噪声后有毛刺,但2%的噪声有点小,毛刺不明显。
subplot(2,1,2), plot(x11);title('加高斯白噪信号');dt=1/1023.5; % 取樣時距(或週期),sec
t=(0:dt:2)'; % 建立一個0-2秒的時間向量
y_high=sin(2*pi*1000*t)/10; % 100 Hz的高頻訊號,振福1/10
y_out=x1+y_high; % 基頻載上一組高頻的輸出。
figure
plot(t,y_out);grid;title('加入高频噪声1000hz');%——————IIR零相移数字滤波器纠正基线漂移——————-
Wp=1.4*2/Fs;     %通带截止频率
Ws=0.6*2/Fs;     %阻带截止频率
devel=0.005;    %通带纹波
Rp=20*log10((1+devel)/(1-devel));   %通带纹波系数
Rs=20;                          %阻带衰减
[N, Wn]=ellipord(Wp,Ws,Rp,Rs,'s');   %求椭圆滤波器的阶次
[b, a]=ellip(N,Rp,Rs,Wn,'high');       %求椭圆滤波器的系数
[hw,w]=freqz(b,a,512);
result =filter(b,a,x1);   figure
freqz(b,a);title('IIR数字滤波器幅频曲线');
figure;subplot(2,1,1); plot(x1);
xlabel('t(s)');ylabel('幅值');title('原始信号');grid
subplot(2,1,2); plot(result);
xlabel('t(s)');ylabel('幅值');title('线性滤波后信号');grid  figure
N=512;
subplot(2,1,1);plot(abs(fft(x1))*2/N);
xlabel('频率(Hz)');ylabel('幅值');title('原始信号x1频谱');grid;
subplot(2,1,2);plot(abs(fft(result))*2/N);
xlabel('频率(Hz)');ylabel('幅值');title('线性滤波去掉基线漂移频谱');grid;
%(5)维纳滤波器去除工频干扰:
dt=1/1023.5; % 取樣時距(或週期),sec
t=(0:dt:2)'; % 建立一個0-2秒的時間向量
y=sin(2*pi*50*t)/10; % 100 Hz的高頻訊號,振福1/10
x_noise=x1+y; % 基頻載上一組高頻的輸出。
figure
plot(t,x_noise);grid;title('加入噪声50hz');%维纳滤波
Mlag=100;
N=100;%维纳滤波器长度
Rxn=xcorr(x_noise,Mlag,'biased');
Rxnx=xcorr(x1,x_noise,Mlag,'biased');%产生输入信号与原始信号的互相关函数
rxnx=zeros(N,1);
rxnx(:)=Rxnx(101:101+N-1);
Rxx=zeros(N,N);%产生输入信号自相关矩阵
Rxx=diag(Rxn(101)*ones(1,N));
for i=2:N
c=Rxn(101+i)*ones(1,N+1-i);
Rxx=Rxx+diag(c,i-1)+diag(c,-i+1);
end
Rxx;
h=zeros(N,1);
h=inv(Rxx)*rxnx;       %计算维纳滤波器的h(n)
yn=filter(h,1,x_noise);        %将输入信号通过维纳滤波器
figure
plot(yn);title('经过维纳滤波器后信号');
ynmean=mean(yn)    %计算经过维纳滤波器后信号均值
ynms=mean(yn.^2)   %计算经过维纳滤波器后信号均方值
ynvar=var(yn,1)    %计算经过维纳滤波器后信号方差
Ryn=xcorr(yn,Mlag,'biased'); %   计算经过维纳滤波器后信号自相关函数
Y=fft(yn);     %计算经过维纳滤波器后信号序列的快速离散傅里叶变换
Y1=fft(x_noise);
Py=Y.*conj(Y)/600; %   计算信号频谱
Py1=Y.*conj(Y)/600;
figure
subplot(211)
semilogy(t,Py)  %绘制在半对数坐标系下频谱图像
title('经过维纳滤波器后信号在半对数坐标系下频谱图像')
xlabel('频率(HZ)','color','b');ylabel('幅度','color','b')
subplot(212)
semilogy(t,Py1) %绘制在半对数坐标系下频谱图像
title('噪声信号在半对数坐标系下频谱图像')
xlabel('频率(HZ)','color','b');ylabel('幅度','color','b')
pyn=periodogram(yn);%计算经过维纳滤波器后信号的功率谱密度
pyn1=periodogram(x_noise);
figure
subplot(211)
semilogy(pyn)%绘制在半对数坐标系下功率谱密度图像
title('经过维纳滤波器后信号在半对数坐标系下功率谱密度图像');
xlabel('频率','color','b');ylabel('幅度','color','b')
subplot(212)
semilogy(pyn);title('噪声信号在半对数坐标系下功率谱密度图像');
xlabel('频率','color','b');ylabel('幅度','color','b')

附录(里面有4组不同病例的心电信号数据,心电信号的特征提取、分析与处理的PPT及处理后报告分析过程)

链接:https://pan.baidu.com/s/1E70b0avxqRU6lIE18Cwe3Q
提取码:oux5

心电信号的特征提取、分析与处理相关推荐

  1. 【C4】基于深度学习的心电信号分析

    ★★★ 本文源自AI Studio社区精品项目,[点击此处]查看更多精品内容 >>> 基于深度学习的心电信号分析 一.项目背景 近年来,随着人工智能和算法的发展,以机器学习和深度学习 ...

  2. arma模型谱估计matlab_基于机器学习的心律失常分类(四)——心电信号特征提取[MATLAB]...

    目前比较常用的特征提取是提取心电信号的各波形间期长度.波峰高度等,本文是使用ARMA模型对心电信号进行处理,使用其系数来作为特征. 一.心拍划分 大多伴有异常波形的心律失常信号,通常都会具体表现在单个 ...

  3. 心音与心电信号分析之一--6.26--心音信号数字滤波

    这才是我要学习的重点哇... 传统的低通,高通滤波,维纳滤波,卡尔曼滤波,自适应滤波属于线性滤波,在消除噪声的同时,往往会造成信号的奇异点改变或信号边缘模糊不清等不利因素,而同态滤波主要解决的是信号与 ...

  4. 【ECG理论篇】(3)AI实现心律失常判别:心电信号的波形识别与特征提取

    心电图中的各个波形都包含了非常多的信息,例如RR间期可以反映心动周期的时限:相邻心动周期的 RR 间期的比值可以反映室性早搏:R 波和 S 波幅值的比值和 R 波和 S 波之间的时限可以反映房性早搏等 ...

  5. 【特征提取】心电信号PTT+HRV+PRV含Matlab源码

    1 简介 ​提取心电信号PTT+HRV+PRV 2 部分代码 clc;clear all;close all;x=load('day2_0917.txt');%% ECG signaly=x(1:95 ...

  6. Matlab心电信号的PQRST模拟-实验报告

    心电信号处理算法设计-实验要求 data4 是一段实际采样得到的心电数据, 采样频率为 100Hz, 波形如下图所示.设计相应的算法, 计算心率, 单位为: 次/分钟.可能会用到的知识为数字滤波器的设 ...

  7. 心电信号去噪(part4)--经验模态分解法(EMD)

    本系列上一篇(数学形态学)指路:https://blog.csdn.net/m0_37422217/article/details/90744326 注:这里是以小型手持心电图机为研究对象的(单导联) ...

  8. 心电电路算法滤波_心电信号噪声的数字滤波研究

    心电信号噪声的数字滤波研究 现代医学表明,心电信号( ECG )含有临床诊断心血管 疾病的大量信息, ECG 的检测与分析在临床诊断中具有重要 价值,是了解心脏的功能与状况.辅助诊断心血管疾病.评 估 ...

  9. 【心电信号】基于matlab心率检测【含Matlab源码 1993期】

    一.心电信号简介 0 引言 心电信号是人类最早研究的生物信号之一, 相比其他生物信号更易于检测, 且具有直观的规律.心电图的准确分析对心脏病的及早治疗有重大的意义.人体是一个复杂精密的系统, 有许多不 ...

  10. 根据心电信号计算心率的matlab代码

    先说说现有的资料: [1][2]弄了FIR1还有hamming窗,搞得老麻烦了,扔一边先. --------------------------------------------------自己另外 ...

最新文章

  1. ArcGis License 启动失败
  2. 壕!甲骨文创始人 8000 万美元买豪宅后打算拆掉
  3. 超级实用sap table
  4. OpenCV演示FloodFill()函数的实例(附完整代码)
  5. 为决战5G时代,小米出手一点不含糊,接连招揽牛人,这次是曾学忠
  6. SAP 前端技术的演化史简介
  7. python3线程池爬虫_python3爬虫中多线程的优势总结
  8. 自制人脸数据,利用keras库训练人脸识别模型
  9. 人物传记——周小川、李稻葵(央行智囊团)
  10. ASP.NET MVC 实现页落网资源分享网站+充值管理+后台管理(15)之前台网站页面
  11. 上位机软件需求说明书100元
  12. java利用反射映射两个不同对象的属性值
  13. 初识Kinect之一
  14. 双向链表:P1996约瑟夫问题的解决方法
  15. INS 、AHRS、VRU、IMU的区别与联系
  16. HCTF 2018-warmup
  17. myEclipse8.5注册码
  18. 演绎另类黑客马拉松,机智云中国第二届智能硬件36小时开发大赛完美收官
  19. 明日边缘:愈演愈烈的POS机网络犯罪
  20. Retrofit2学习项目_2

热门文章

  1. matlab 求平面方程,MATLAB求空间平面方程
  2. JAVA EE常见英文单词(一)
  3. Windows 上 curl 的安装
  4. 手术导航系统原理简介、主要工作及应用
  5. arcgis重分类读不出值的解决办法?
  6. 学习DSP28335--CCS软件打开例程时一直报错问题以及编译问题解决
  7. DSP入门小白学习日记第四篇
  8. eclipse php 教程,Eclipse PHPEclipse配置方法文字教程
  9. 关于adb no serial number的解决方案
  10. Spyder使用教程