数字信号处理选做实验--心电信号数字滤波系统设计

F0=F0';
t1=t1';
figure(8)
subplot(3,1,1);
plot(1000*t1,F0);
wp=0.7*2*pi;ws=0.25*2*pi;Rp=0.1;As=40;T=1;
[N,wc]=buttord(wp,ws,Rp,As,'s')
[B,A]=butter(N,wc,'s');
title('低通滤波后信号频谱_{2.5}')
4.2．Simulink仿真参数设置及波形对比：
六、实验结果分析
subplot(3,1,3);plot(y21);
title('低通滤波后_{2.5}')
M=length(F0);
m=1:M;
fs=1000;
f2=m*fs/M;
F01=fft(F0);
Y11=fft(y11);
Y21=fft(y21)
figure(9)
subplot(3,1,1)
plot(f2,abs(F01))
F1=fft(F);
Y1=fft(y1);
Y2=fft(y2)
figure(2)
subplot(3,1,1)
plot(f2,abs(F1))
axis([0,1000,0,200])
title('原始信号频谱_{59.997}')
subplot(3,1,2)
plot(f2,abs(Y1))
axis([0,1000,0,200])
[BH,AH]=lp2hp(B,A,wph);
[h,w]=freqs(BH,AH);
plot(w,20*log10(abs(h)));
axis([0,1,-80,5]);
grid on
N =
7
wc =
3.0327
B0 =
1.0e+003 *
0 0 0 0 0 0 0 2.3595
A0 =
1.0e+003 *
[b1,a1]=imp_invr(cs1,ds1,T)
[db1,mag1,pha1,w1]=freqz_m(b1,a1);
y2=filter(b1,a1,y1);
subplot(3,1,3);plot(y2);
title('低通滤波后')
M=length(F);
m=1:M;
fs=1000;
f2=m*fs/M;
本文是利用MATLAB软件对美国麻省理工学院提供的MIT-BIH数据库的122号心电信号病例进行分析，利用MATLAB软件及simulink平台对122号心电信号的病例进行读取、插值、高通滤波、低通滤波等的处理。将心电信号中的高频和低频的杂波进行滤除后对插值前后滤波前后的时域波形及频谱进行分析。同时也将滤波器的系统函数进行读取，分析，画出滤波的信号流程图，也分析各个系统及级联后的系统的幅频响应、相位响应来判断系统的稳定性，并用MATLAB软件将图形画出，以便于以后的对系统进行分析。
五、实验过程原始记录（数据、图标、计算等）
1.提取txt格式心电信号：
fid=fopen('122.txt');
C=textscan(fid,'%8c %f %*f','headerlines',2);
fclose(fid);
a=C{1};
y=C{2};
k=length(a)
for i=1:k
c(i)=strread(a(i,:),'%*s %f','delimiter',':');
wp=0.7*2*pi;ws=0.25*2*pi;Rp=0.1;As=40;T=1;
[N,wc]=buttord(wp,ws,Rp,As,'s')
[B,A]=butter(N,wc,'s');
[b,a]=imp_invr(B,A,T)
[db,mag,pha,w]=freqz_m(b,a);
y1=filter(b,a,F);
3.对处理前后的心电信号分别做频谱分析
利用Matlab软件对处理前后的心电信号编程显示其频谱，分析比较滤波前后的频谱，得出结论。
如果分析频谱，滤波效果不明显，则需变动滤波器参数指标，重新设计滤波器。通过频谱分析，多次试验确定最合适的滤波器。
4 .Simulink仿真.
根据前面的设计，进行基于Simulink的动态仿真设计。实现心电信号的分析和处理。给出系统的基于Simulink的动态建模和仿真的系统方框图，同时记录系统的各个输出点的波形和频谱图。
实验题目
心电信号数字滤波系统设计
实验时间
2013年12月16日
实验地点
A1118
实验成绩
实验性质
□验证性□设计性□综合性
教师评语：
□出勤率好□原理正确□方案合理
□实验结果正确□回答问题正确□报告规范
一、实验目的
本实验项目主要研究数字心电信号的初步分析及滤波器的应用。通过完成本பைடு நூலகம்验的设计，主要达到以下几个目的：
Columns 11 through 12
0.0926 5.3949
6模拟高通滤波器：
wp=0.7*2*pi;ws=0.25*2*pi;Rp=0.1;As=40;
[N,wc]=buttord(wp,ws,Rp,As,'s')
[B0,A0]=butter(N,wc,'s');
wph=2*pi*0.25;hk=freqs(B0,A0,wph);
三、使用仪器、器材
Matlab软件
四、实验步骤
1.对原始心电信号做线性插值
运算公式如下:
其中是第i个数据时间点，Ai是与之对应的数据，N是两数据之间需要的插值数，是需要插值的两点数据差，，
2.根据心电信号的频域特征，设计相应的低通和带通滤波器
一般正常人的心电信号频率在0.7～100HZ范围内，人体心电信号微弱，信噪比小，因此，在采集心电信号时，易受到仪器、人体活动等因素的影响，而且所采集的心电信号常伴有干扰。采集心电数据时，由于人的说话呼吸，常常会混有约为0.1Hz到0.25Hz频段的干扰，对于这些低频干扰，可以让信号通过一个高频滤波器，低截止频率设置为0.25,来滤除低频信号，对于高频信号干扰，可以让信号再通过一个低频滤波器，其中截止频率设置为99H
axis([0,1000,0,200])
title('原始信号频谱_{2.5}')
subplot(3,1,2)
plot(f2,abs(Y11))
axis([0,1000,0,200])
title('高通滤波后信号频谱_{2.5}')
subplot(3,1,3)
plot(f2,abs(Y21))
axis([0,1000,0,200])
m=1:M;
f2=m*fs/M;
Y2=fft(F);
subplot(2,2,4)
plot(f2,abs(Y2))
title('插值后信号频谱')
axis([0 1000 0 200])
5模拟低通滤波器：
wp=60*2*pi;ws=99*2*pi;Rp=1;As=40;
[N,wc]=buttord(wp,ws,Rp,As,'s')
fid = fopen('F.txt','wt');
fprintf(fid,'%g\n',F);
fclose(fid);
4.插值前后波形比较
subplot(2,2,1)
plot(x,y)
title('初始信号时域波形')
axis([0 2.5 -2 1])
subplot(2,2,2)
fs=1000;
end
x=c';
plot(x,y)
2.对原始心电信号进行线性插值
t=0.006:0.001:2.5;
F=interp1(x,y,t);
F=F';
t=t';
plot(t,F)
3.把数据读到txt中
fid = fopen('t.txt','wt');
fprintf(fid,'%g\n',t);
fclose(fid);
N=length(y)
n=1:N;
f1=n*fs/N;
Y1=fft(y);
plot(f1,abs(Y1))
title('初始信号频谱')
axis([0 1000 0 200])
subplot(2,2,3)
plot(t,F)
title('差值后信号时域波形')
axis([0 2.5 -2 1])
M=length(F);
[b,a]=imp_invr(B,A,T)
[db,mag,pha,w]=freqz_m(b,a);
y11=filter(b,a,F0);
subplot(3,1,2);plot(y11);
title('高通滤波后_{2.5}')
wp1=2*pi*60;ws1=2*pi*99;Rp1=0.1;As1=40;T1=1000;
OmegaP1=wp1/T1;OmegaS1=ws1/T1;
[cs1,ds1]=afd_butt(OmegaP1,OmegaS1,Rp1,As1);
[b1,a1]=imp_invr(cs1,ds1,T)
[db1,mag1,pha1,w1]=freqz_m(b1,a1);
y21=filter(b1,a1,y11);
subplot(3,1,2);plot(y1);
title('高通滤波后')
wp1=2*pi*60;ws1=2*pi*99;Rp1=0.1;As1=40;T1=1000;
OmegaP1=wp1/T1;OmegaS1=ws1/T1;
[cs1,ds1]=afd_butt(OmegaP1,OmegaS1,Rp1,As1);
Columns 1 through 10
0 0 0 0 0 0 0 0 0 0
Columns 11 through 12
0 5.3949
A =
1.0e+28 *
Columns 1 through 10
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0008
1、了解MATLAB软件的特点和使用方法，熟悉基于Simulink的动态建模和仿真的步骤和过程。
2、了解人体心电信号的时域特征和频谱特征。
3、进一步了解数字信号的分析方法。
4、通过应用具体的滤波器进一步加深对滤波器理解。
5、通过本课题的设计，培养学生运用所学知识分析和解决实际问题的能力。
二、实验原理
B,A]=butter(N,wc,'s')
k=0:511;fk=0:1000/512:1000;wk=2*pi*fk;
Hk=freqs(B,A,wk);
plot(fk,20*log10(abs(Hk)));
grid on
>> Untitled5
N =
11
wc =
409.2596
B =
1.0e+28 *
0.0010 0.0136 0.0929 0.4070 1.2343 2.5905 3.4964 2.3595
7滤波前后图形对比：
t=0.006:0.001:59.997;
F=interp1(x,y,t);
F=F';
t=t';
figure(1)
subplot(3,1,1);
plot(1000*t,F);
title('高通滤波后信号频谱_{59.997}')
subplot(3,1,3)
plot(f2,abs(Y2))
axis([0,1000,0,200])
title('低通滤波后信号频谱_{59.997}')
8.截取到2.5s对截取的部分进行滤波及频谱分析
t1=0.006:0.001:2.5;
F0=interp1(x,y,t1);