医学信号处理实验报告心脑电信号认知Word文件下载.docx
《医学信号处理实验报告心脑电信号认知Word文件下载.docx》由会员分享，可在线阅读，更多相关《医学信号处理实验报告心脑电信号认知Word文件下载.docx（15页珍藏版）》请在冰点文库上搜索。
3）对于所有大于阈值的峰值点作为检测到的R波尖峰。
4）由生理基础可以知道，R波间隔是相对稳定的。
可以通过检测峰值点的间隔，去除那些较高的伪迹。
3.Pan-Tompkins法检测心电R波尖峰
1）将信号分别通过给定的低通滤波器、高通滤波器
2）对滤波后的信号求一阶导数
3）对求导之后的信号进行平方运算
4）将信号通过滑动窗口进行积分，这里选取窗口长度为30
5）应用阈值法检测经过前四步处理之后的心电信号R波尖峰
流程图如下图所示
其演示效果如下图所示
5、实验目的:
1）周期图法的改进方法，和分段平均对图像数据的影响。
2）能够利用两种方法处理心电波形并计算一些特征值。
6、实验内容：
（一）上机题3：
改进周期图法估计功率谱
1、接着上机题2做，任选一种窗函数，用分段、平均的思想改进周期图法，观察改进前后功率谱的差异；
2、给出一段文字总结周期图法的缺点，改进法的优点。
（二）上机题4：
心电R波检测和RR间隔估计
使用阈值法和Pan-Tompkins的检测方法验证信号使用数据：
ECG3.dat、ECG4.dat、ECG5.dat和ECG6.dat，采样率为200Hz（参考文件ECGS.m）。
计算每个数据的RR波间隔和心率的平均值。
比较两种方法结果的差异。
分别画出四个数据，并通过视觉直接观察法测量其参数来验证估计结果。
七、实验器材（设备、元器件）：
matlab2014b
八、实验步骤：
构思，设计算法，写程序，写报告。
9、实验数据及结果分析：
（一）1.程序
clearall
close=load（'
eegclose'
）;
open=load（'
eegopen'
close=close.eegclose（:
8*2）;
open=open.eegopen（:
n=length（close）;
fs=250;
x=（0:
n-1）'
/fs;
figure
（1）
subplot（2,2,1）
plot（x,close）;
xlabel（'
时间/s'
ylabel（'
幅度/mV'
title（'
闭眼信号'
）
subplot（2,2,2）
plot（x,open）;
睁眼眼信号'
len=500;
%序列长度等于4500-len
d_num=100;
%分段数为len/d_num
%构造
win=cell（1,4）;
rect=boxcar（n-len）;
win{1,1}=rect;
trian=triang（n-len）;
win{1,2}=trian;
hamm=hamming（n-len）;
win{1,3}=hamm;
bman=blackman（n-len）;
;
win{1,4}=bman;
close_w=cell（1,4）;
open_w=cell（1,4）;
fori=1:
length（win）
close_w{1,i}=zeros（n-len,1）;
open_w{1,i}=zeros（n-len,1）;
forj=1:
d_num:
len;
close_w{1,i}=10*log10（abs（fft（close（j:
n-len-1+j）.*win{1,i}））.^2/（n-len））+close_w{1,i};
open_w{1,i}=10*log10（abs（fft（open（j:
n-len-1+j）.*win{1,i}））.^2/（n-len））+open_w{1,i};
end
end
%功率谱平均
close_w{1,i}=close_w{1,i}/len;
open_w{1,i}=open_w{1,i}/len;
freq=（0:
（n-len）/2-1）/（n-len）*fs;
names={'
矩形窗'
'
三角窗'
海明窗'
凯撒窗'
};
%在功率谱中标注出δ（1-3Hz）、θ（4-7Hz）、α（8-13Hz）、β（14-30Hz）
locate_mid=[1.5,5.5,10.5,22,45];
%生成四种脑电波的定位矩阵
loc={[1,3][4,7][8,13][14,30]};
4
subplot（4,2,2*i-1）
plot（freq,close_w{1,i}（1:
（n-len）/2））
xlabel（'
频率/Hz'
db'
title（['
闭眼加'
names{1,i},'
的功率谱'
]）;
axis（[0100min（close_w{1,i}）max（close_w{1,i}）]）
set（gca,'
XTickmode'
manual'
Xtick'
locate_mid）
subplot（4,2,2*i）
plot（freq,open_w{1,i}（1:
ylabel（'
睁眼加'
axis（[0100min（close_w{1,i}）max（open_w{1,i}）]）
2.结果：
1）每段500点：
图1.4000段平均
图2.400段平均
图3.40段平均
图4.4段平均
2）每段2500点:
图5.20段平均
（二）1.程序
ecg=cell（1,4）;
peaks=cell（1,4）;
becg=cell（1,4）;
rs=zeros（1,4）;
hb=zeros（1,4）;
fori=3:
6
ecg{1,i-2}=load（['
ECG'
num2str（i）,'
.dat'
peaks{1,i-2}=ECGfindpeaks（0.4,200,ecg{1,i-2}）;
subplot（2,2,i-2）
plot（peaks{1,i-2}（:
1）,peaks{1,i-2}（:
2）,'
o'
（1:
length（ecg{1,i-2}））/200,ecg{1,i-2}）
mean（diff（peaks{1,i-2}（:
1）））
rs（i-2）=mean（diff（peaks{1,i-2}（:
hb（i-2）=60/（mean（diff（peaks{1,i-2}（:
1））））
num2str（i）]）;
电压/uv'
[N,Wn]=buttord（40/200,60/200,1,32）;
[b1,a1]=butter（N,Wn）;
[N,Wn]=buttord（10/200,30/200,1,32）;
[b2,a2]=butter（N,Wn,'
high'
figure
（2）
becg{1,i-2}=filtfilt（b1,a1,ecg{1,i-2}）;
becg{1,i-2}=filtfilt（b2,a2,becg{1,i-2}）;
becg{1,i-2}=abs（diff（becg{1,i-2}））;
becg{1,i-2}=[zeros（30-1,1）;
becg{1,i-2};
zeros（30-1,1）];
r=zeros（length（becg{1,i-2}）-57,1）;
（length（becg{1,i-2}）-57）
r（j）=sum（becg{1,i-2}（j:
（j-1+30）））;
becg{1,i-2}=r;
plot（r）
R=peaks1（becg{1,i-2},20,0.4）-20;
subplot（2,2,i-2）
plot（R/200,r（R）,'
length（r））/200,r）
rs（i-2）=mean（diff（R））/200
hb（i-2）=60/（mean（diff（R））/200）
function[z_r]=ECGfindpeaks（t,fs,ecgdata）
sze=length（ecgdata）;
n=fix（sze/（t*fs））;
z=zeros（n-1,2）;
z_l=1;
forj=1:
t*fs
fori=1:
n-1
[peks,loca]=findpeaks（ecgdata（（fs*t）*（i-1）+j:
（fs*t）*i+j））;
x=find（peks>
0.5*max（peks）+0.5*mean（peks））;
loca=loca（x）+fs*t*（i-1）+1;
loca_l=length（loca）;
%plot（（loca+j）/fs,peks（x）,'
time,necg）;
%holdon;
z（z_l:
z_l+loca_l-1,1）=（loca+j）/fs;
z_l+loca_l-1,2）=peks（x）;
z_l=z_l+loca_l;
ifj==1
rpeks=z;
c=1;
z_r=zeros（length（rpeks）,2）;
length（rpeks）
z_1=find（z（:
1）==rpeks（i,1））;
iflength（z_1）>
=t*fs-1
z_2=find（z（:
2）==rpeks（i,2））;
iflength（z_2）>
z_r（c,:
）=[rpeks（i,1）rpeks（i,2）];
c=c+1;
cc=1;
fori=2:
c
if（z_r（i,1）-z_r（cc,1））<
=t-0.1;
z_r（i,1）=0;
elsecc=i;
z_r=z_r（z_r（:
1）>
0,:
）;
functionR=peaks1（ecg,m,th）
ecg_f=[zeros（1,m）'
ecg;
zeros（1,m）'
];
ecg_max=max（ecg_f）;
ecg_min=min（ecg_f（ecg_f>
0））;
th=ecg_min+（ecg_max-ecg_min）*th;
ind_ecg=find（ecg_f>
th）;
R=[];
fori=1:
length（ind_ecg）
ifecg_f（ind_ecg（i））==max（ecg_f（（ind_ecg（i）-m）:
（ind_ecg（i）+m）））
R=[R,ind_ecg（i）];
if（length（R）>
1）
if（（R（end）-R（end-1））<
70）
R（end）=[];
1）阈值法
RR间期（单位：
秒）：
0.81590.63050.46111.0953
心率（单位：
次/分钟）：
73.541295.1626130.124354.7798
2）Pan-Tompkins法检测心电R波尖峰
0.81080.63050.45051.0953
73.997995.1650133.195754.7798
10、总结及心得体会：
（一）由几次实验结果可以看出：
改进的周期图法结果的差异主要见于分段的长度，而非平均的段数。
并且段越长，dB图像越平滑，也就是说频率分辨率越差，但合适的图像平滑效果正是我们想要的结果。
（二）
改进阈值法的结果很有趣味，所以我多做了一些尝试来使得其检验效果更好，但是经过改良的阈值法，就不能适用于Pan-Tompkins法。
我不得不在Pan-Tompkins法里使用最开始的方法。
1.阈值法：
阈值法可以很好的检测规律的波形，但是当还有基线漂移的时候，其开始变得不稳定。
所以目标和方法就在于抵抗基线漂移。
（虽然滤除低频分量可以很好的改善基线漂移，但那并不在考虑范围）
1）因为R波波峰一般是最高的。
故假设，误检的发生其概率是小于检测到R波峰的。
如果数据长为100*n（n为一次心动周期的近似数据长度），设计一个大窗长度为99*n，牺牲掉n点长的数据（牺牲更多数据并没有好处），以小窗长度n进行分段计算，得到n份R波尖峰数据，设置阈值采信在这n分数据中出现th*n以上的数据。
这样可以避免在较缓慢基线漂移中，因为小窗位置不当而造成的个别失检和误检。
2）在阈值的设定中，将阈值设为a*（max（g（n）））+（1-a）（mean（g（n）））,（g（n）为处理数据时小窗内的信号数据）。
将均值加入函数中使得基线的漂移得到体现。
3）但对于剧烈的基线漂移导致的伪迹比R波更像R波，改进检测算法无能为力。
只能用检测时间间隔的方法剔除掉，从第一个R波峰值开始，如果下一个值与其距离在0.4s以内，则去除该点。
若不是，则将该点作为开始点重复上面的操作。
但这样至少要保证第一个点是R波峰值，这很难做到。
但是基于第一条的假设下，我们可以选择一个在99条数据中出现次数最多的R波峰值作为是第一个点。
但是由于数据太短，没能验证想法。
理论上来说这种方法是具有积极意义的。
2.Pan-Tompkins法
滤波加积分的方式使得峰值的位置发生了稍稍的偏移，这是使其结果和阈值法出现区别的原因。
值得一提的是，由于Pan-Tompkins法滤波器的选择对波形影响很大，会在上升和下降的边沿出现毛刺。
加之波形密集能量集中，使得为阈值法写的函数出现很大误差。
.
十一、对本实验过程及方法、手段的改进建议：
无
报告评分：
指导教师签字：
升级会员 