基于matlab的数字滤波器设计及心电信号分析.docx
《基于matlab的数字滤波器设计及心电信号分析.docx》由会员分享，可在线阅读，更多相关《基于matlab的数字滤波器设计及心电信号分析.docx（18页珍藏版）》请在冰点文库上搜索。
基于matlab的数字滤波器设计及心电信号分析
设计题目
基于matlab的数字滤波器设计及心电信号分析
设
计
技
术
参
数
阈值=0.6相对幅值的差
RR间期1mv
设
计
要
求
能够准确的提取与分析QRS波
绘制图表来说明检测结果
工
作
量
学会使用matlab软件
用数字滤波器测量心电信号
对心电信号有个初步的认识与掌握
搜集一组正常人体心电信号数据，编写matlab程序，并对它进行分析与检测
根据设计结果写出报告
工
作
计
划
查资料，搜集关于此次课设的相关内容
对题目进行分析，确定合理的解决方案
编程并且调试
写报告
参
考
资
料
周辉数字信号处理基础及Matlab实现中国林业出版社2005
肖伟刘忠Matlab程序设计与应用清华大学出版社2005
指导教师签字
孟辉
基层教学单位主任签字
李昕
1课程设计的目的
通过对常用数字滤波器的设计和实现，掌握数字信号处理的工作原理及设计方法；熟悉设计数字滤波器的方法，掌握利用数字滤波器对信号进行滤波的方法，掌握数字滤波器的计算机仿真，并能够对设计结果加以分析。
2课程设计的原理
2.1、用窗函数法设计FIR滤波器
根据过渡宽带及阻带衰减要求，选择窗函数的类型并估计窗口长度N，窗函数类型可根据最小阻带衰减As独立选择，因为窗口长度N对最小阻带衰减As没有影响，在确定窗函数类型以后，可根据过渡宽带小于给定指标确定所拟用的窗函数的窗口长度N，设待求滤波器的过渡宽带，它与窗口长度N近似成反比，窗函数类型确定后，其计算公式也确定了，不过这些公式是近似的，得出的窗口长度还要在计算中逐步修正，原则上在保证阻带衰减满足要求的情况下，尽量选择较小的N，在N和窗函数类型确定后，即可调用MATLAB中的窗函数求出窗函数。
2.2、用巴特沃斯法设计IIR低通滤波器
巴特沃斯低通滤波器的幅度平方函数中，N为滤波器的阶数。
在Ω=Ωc附近，随着Ω加大，幅度迅速下降。
幅度下降的速度与阶数N有关，N越大，通带愈平坦，过渡带愈窄，过渡带与阻带幅度下降的速度愈快，总的频响特性与理想低通滤波器的误差愈小。
3课程设计设计步骤及结果分析
3.1、心电信号采集
心电信号作为心脏电活动在人体体表的表现，信号一般比较微弱，幅度在10μV~5mV，频率为0.05~100Hz。
在心电信号的采集、放大、检测及记录过程中，有来自外界的各种干扰。
记录一段时间内的人体心电信号波形，要求长度不小于10秒，并对记录的信号进行数字化，保存为数据文件；这里，请同学们使用美国的MIT/BIH心电原始数据，由实验老师给出一定长度的的心电原始数据，数据保存在文件“a01.txt～a10.txt”中，在MATLAB中通过如下语句读取：
%从当前路径下的a01.txt文件读取心电原始数据到变量a01中，a01为二维数据，第一列%为心电信号时间，第二列为心电信号幅度。
3.2、源数据的导入
a=load（'F:
/心电信号数据/a16.txt'）;%选第十六组数据
3.3、绘出源心电信号的时域波形图和频谱图并进行分析
%时域波形图以及频谱图
a=load（'1.txt'）;
t=a（:
1）;b=a（:
2）;
figure
（1）;
subplot（211）;
plot（t,b）;
title（'心电信号的时域波形'）;
xlabel（'时间t/s'）;ylabel（'幅值/A'）;
n=1000;
m=abs（fft（b,n））;
fs=100;
f=fs/n*（0:
n-1）;
subplot（212）;
plot（f,m）;
title（'心电信号的频谱图'）;
xlabel（'频率f/Hz'）;ylabel（'幅值/db'）;
3.4、含噪心电信号合成
1、加入单频正弦干扰
y=0.2*sin（2*40*pi*t）;
y1=y+b;
figure
（2）;
subplot（211）;
plot（t,y1）;
title（'加单频噪声时域波形'）;
xlabel（'时间t/s'）;ylabel（'幅值/A'）;
k=abs（fft（y1,n））;
subplot（212）;
plot（f,k）;
title（'加单频噪声频谱图'）;
xlabel（'频率/Hz'）;ylabel（'幅值'）
选择添加单频正弦信号。
通过转置后添加到原始信号b上，添加后的信号用y1表示。
对于变量、函数的操作是矩阵操作，两个信号相加必须长度一致。
2、继续加入白噪声
y2=awgn（y1,5）;
figure（3）;
subplot（211）;
plot（t,y2）;
title（'加单频,白噪声时域波形'）;
xlabel（'时间t/s'）;ylabel（'幅值/A'）;
h=abs（fft（y2,n））;
subplot（212）;
plot（f,h）;
title（'加单频,白噪声频谱图'）;
xlabel（'频率/Hz'）;ylabel（'幅值'）;
通过调用函数y2=awgn（y1,5）对已经加了单频噪声的信号进行加噪，加噪后的信号用y2保存。
对y2信号的分析，如图所示分别给出了它的时域波形图和频谱图，通过对比，时域波形变化不大，但频谱图有了大幅度的变化。
3.5、时域波形差
s=y2-b;
figure（4）;
subplot（111）;
plot（t,s）;
title（'时域波形差'）;
3.6、数字滤波器设计及滤波
1、设计低通FIR滤波器
%低通FIR滤波器
N=50;
wc=0.6;
window=blackman（N）;
hn=fir1（N-1,wc,window）;
lv=filter（hn,1,y2）;
freqz（hn,1）;
2、用低通FIR滤波器对单频噪声进行滤波
%低通滤波器滤波后图形
subplot（211）;
plot（t,lv）;
title（'滤波后时域波形'）;
xlabel（'时间t/s'）;ylabel（'幅值/A'）;
LV=abs（fft（lv））;
subplot（212）;
plot（f,LV）;
title（'滤波后频谱图'）;
xlabel（'频率/Hz'）;ylabel（'幅值'）
3、设计带阻FIR滤波器
%带阻FIR滤波器
wlp=0.2*pi;wls=0.4*pi;wus=1.2*pi;wup=0.6*pi;
B=wls-wlp;
M=ceil（12*pi/B）;
wp=[（wls+wlp）/2/pi,（wus+wup）/2/pi];
hn=fir1（M,wp,'stop',blackman（M+1））;
liu=filter（hn,1,y2）;
freqz（hn,1）;
4、用带阻FIR滤波器对单频噪声滤波后图形
%带阻FIR滤波器滤波后图形
subplot（211）;
plot（t,liu）;
title（'滤波后时域波形'）;
xlabel（'时间t/s'）;ylabel（'幅值/A'）;
LIU=abs（fft（liu））;
subplot（212）;
plot（f,LIU）;
title（'滤波后频谱图'）;
xlabel（'频率/Hz'）;ylabel（'幅值'）;
5、设计低通IIR滤波器
%低通IIR滤波器
wp=0.4*pi;ws=0.35*pi;rp=1;rs=10;
[N,wc]=buttord（wp,ws,rp,rs,'s'）;
[B,A]=butter（N,wc,'s'）;
[Bz,Az]=impinvar（B,A）;
wen=filter（Bz,Az,y2）;
freqz（Bz,Az）;
6、用低通IIR滤波器对单频噪声滤波后图形
%低通IIR滤波器滤波后图形
subplot（211）;
plot（t,wen）;
title（'滤波后时域波形'）;
xlabel（'时间t/s'）;ylabel（'幅值/A'）;
WEN=abs（fft（wen））;
subplot（212）;
plot（f,WEN）;
title（'滤波后频谱图'）;
xlabel（'频率/Hz'）;ylabel（'幅值'）;
3.7、心电信号时域波形及频谱的观察和分析
1、低通FIR滤波器滤波效果
subplot（331）;
plot（t,b）;
title（'心电信号的时域波形'）;
xlabel（'时间t/s'）;ylabel（'幅值/A'）;
subplot（332）;
plot（f,m）;
title（'心电信号的频谱图'）;
xlabel（'频率f/Hz'）;ylabel（'幅值/db'）;
subplot（333）;
plot（t,y1）;
title（'加单频噪声时域波形'）;
xlabel（'时间t/s'）;ylabel（'幅值/A'）
subplot（334）;
plot（f,k）;
title（'加单频噪声频谱图'）;
xlabel（'频率/Hz'）;ylabel（'幅值'）;
subplot（335）;
plot（t,y2）;
title（'加单频,白噪声时域波形'）;
xlabel（'时间t/s'）;ylabel（'幅值/A'）;
subplot（336）;
plot（f,h）;
title（'加单频,白噪声频谱图'）;
xlabel（'频率/Hz'）;ylabel（'幅值'）;
subplot（337）;
plot（t,lv）;
title（'滤单频波后时域波形'）;
xlabel（'时间t/s'）;ylabel（'幅值/A'）;
subplot（338）;
plot（f,LV）;
title（'滤单频波后频谱图'）;
xlabel（'频率/Hz'）;ylabel（'幅值'）;
2、带阻FIR滤波器滤波效果
subplot（331）;
plot（t,b）;
title（'心电信号的时域波形'）;
xlabel（'时间t/s'）;ylabel（'幅值/A'）;
subplot（332）;
plot（f,m）;
title（'心电信号的频谱图'）;
xlabel（'频率f/Hz'）;ylabel（'幅值/db'）;
subplot（333）;
plot（t,y1）;
title（'加单频噪声时域波形'）;
xlabel（'时间t/s'）;ylabel（'幅值/A'）
subplot（334）;
plot（f,k）;
title（'加单频噪声频谱图'）;
xlabel（'频率/Hz'）;ylabel（'幅值'）;
subplot（335）;
plot（t,y2）;
title（'加单频,白噪声时域波形'）;
xlabel（'时间t/s'）;ylabel（'幅值/A'）;
subplot（336）;
plot（f,h）;
title（'加单频,白噪声频谱图'）;
xlabel（'频率/Hz'）;ylabel（'幅值'）;
subplot（337）;
plot（t,liu）;
title（'滤单频波后时域波形'）;
xlabel（'时间t/s'）;ylabel（'幅值/A'）;
subplot（338）;
plot（f,LIU）;
title（'滤单频波后频谱图'）;
xlabel（'频率/Hz'）;ylabel（'幅值'）;
3、低通IIR滤波器滤波效果
subplot（331）;
plot（t,b）;
title（'心电信号的时域波形'）;
xlabel（'时间t/s'）;ylabel（'幅值/A'）;
subplot（332）;
plot（f,m）;
title（'心电信号的频谱图'）;
xlabel（'频率f/Hz'）;ylabel（'幅值/db'）;
subplot（333）;
plot（t,y1）;
title（'加单频噪声时域波形'）;
xlabel（'时间t/s'）;ylabel（'幅值/A'）
subplot（334）;
plot（f,k）;
title（'加单频噪声频谱图'）;
xlabel（'频率/Hz'）;ylabel（'幅值'）;
subplot（335）;
plot（t,y2）;
title（'加单频,白噪声时域波形'）;
xlabel（'时间t/s'）;ylabel（'幅值/A'）;
subplot（336）;
plot（f,h）;
title（'加单频,白噪声频谱图'）;
xlabel（'频率/Hz'）;ylabel（'幅值'）;
subplot（337）;
plot（t,wen）;
title（'滤单频波后时域波形'）;
xlabel（'时间t/s'）;ylabel（'幅值/A'）;
subplot（338）;
plot（f,WEN）;
title（'滤单频波后频谱图'）;
xlabel（'频率/Hz'）;ylabel（'幅值'）;
4心得体会
通过这次的课程设计，我从中得到了许多经验和MATLAB软件设计的一些新思路；在学习MATLAB的过程中，我明白了“实践出真知”这句话的真谛，在书上看到的内容，如果不懂的话，我们可以把程序输入代码窗口，运行之后，看得出的结果，然后再加上书本上的讲解，就很好理解了，相反，如果我们只是想着这个程序或者这个函数有什么功能，肯定想很长时间都不明白。
通过这次课程设计，MATLAB的使用不仅仅只是看书就能够学会的，要自己在计算机上面动手操作，才能熟练的使用MATLAB软件，我还发现了MATLAB一个强大的功能是它为我们提供一套功能强大的绘图命令，这些命令可以根据输入的数据自动完成图形的绘制，可以完成对图形的加标号，加标题等操作。
总结一下这次课程设计，发现自己虽然在不仅在理论上没有掌握牢固，并且在实践的时候也遇到了问题，所以自己还是远远的不足，不管是在MATLAB的设计上，还是其他专业课上，在以后的一段学习时间里必须坚持自己思考，自己多动脑，多动手，这样才能脱离理论，让自己的学习更上一层楼。
升级会员 