北工大通信系统实验

1、信号处理与分析课程设计训练报告- ECG信号高频噪声消除学号：姓名：指导老师：张新峰完成日期：2016年6月5日绪论本实验是信号处理相关方法在心电图ECG（electocardiogram）信号处理中的一个具体应用。通过心电图进行疾病诊断是心血管疾病诊断的一种重要方法。该方法可靠简便，在临床中应用很广泛。近年来给予数字信号处理和模式识别的ECG信号自动分析成为一个研究这点。ECG信号在采集过程中收到工频信号、肌电信号以及基线漂移等干扰，由此带来的噪声对自动分析的结果影响很大。因此在对ECG信号分析之前，必须进行噪声消除研究。本实验内容给定ECG信号，合理设计滤波器，消除高频噪声的影响。一、实验

2、内容心电信号是微弱低频人体生理电信号，通常频率在0.05100Hz，幅值不超过4mv，通过安装在皮肤表面的电极来获取。由于人体是一个复杂的生 命系统，存在50H 工频干扰及基线漂移等其他生理电信号的干扰。噪声可能会影响到医生的临床诊断，因此，需对心电信号进行滤波，即必须做好前端数据采集的软硬件设计以保证心电数据的可靠和准确。（本次实验的心电数据是老师提供，为txt文档）传统医疗设备分别采用50Hz 带阻滤波器和RC 高通滤波器滤除工频干扰和基线漂移。但带阻滤波器电路复杂，其特性对元器件的精度敏感，而基线漂移本质上是一种缓慢变化的低频信号，采用RC 滤波器很难将高通滤波器的过渡带做得十分陡峭，基

3、线漂移补偿效果不理想。因此，模拟方法往往不太容易获得很好的特性。在本实验中，我们以C语言为基础，来实现高频噪声的消除。二、实验原理1. 低通巴特沃兹模拟滤波器的实现1）概述巴特沃兹滤波器的特点是通频带内的频率响应曲线最大限度平坦，没有起伏，而在阻频带则逐渐下降为零。 在振幅的对数对角频率的波得图上，从某一边界角频率开始，振幅随着角频率的增加而逐步减少，趋向负无穷大。 一阶巴特沃兹滤波器的衰减率为每倍频6分贝，每十倍频20分贝。二阶巴特沃兹滤波器的衰减率为每倍频12分贝、 三阶巴特沃兹滤波器的衰减率为每倍频18分贝、如此类推。巴特沃兹滤波器的振幅对角频率单调下降，并且也是唯一的无论阶数，振幅对角

4、频率曲线都保持同样的 形状的滤波器。只不过滤波器阶数越高，在阻频带振幅衰减速度越快。其他滤波器高阶的振幅对角频率图和低级数的振幅对角频率有不同的形状。2）设计巴特沃兹模拟低通滤波器1确定参数:，2求滤波器阶次N： = N= 4确定系统函数N为奇数时： N为偶数时： 双线性Z法把模拟滤波器转换成数字滤波器3. FFT的算法按时间抽选（DIT）的基-2 FFT算法。当输入序列点数为 （L为整数），先将序列x(n)按n的奇偶分成以下两组：则可以将DFT化为： k=0,1,2, 其中 分别为 点DFT。又根据系数的周期性可得到X(k)的完整表达：三、实现方法流程图1. 读取心电信号数据的实现float

5、 agTemp21024+2; /定义一个数组用来存放数据FILE *p; /定义指针p=fopen("119_heart1.txt","r"); / 打开心电图数据文件for(i=0;i<1024+2;i+) / 通过循环放入数组中fscanf(p,"%f",&agTemp2i);fclose(p);2. FFT变换的实现a. 倒位序处理，先将xn进行倒位序后，才能进行fft变换 void change() complex temp; unsigned short i=0,j=0,k=0; double t; for(

6、i=0;i<Num;i+) k=i;j=0; t=(log(Num)/log(2); while( (t-)>0 ) j=j<<1; j|=(k & 1); k=k>>1; if(j>i) temp=xi; xi=xj; xj=temp; b. 进行蝶形运算，将保存在全局复数数组xn中的时域信号转为频域信号:void fft() int i=0,j=0,k=0,l=0; complex up,down,product; change(); for(i=0;i< log(Num)/log(2) ;i+) /*一级蝶形运算*/ l=1<

7、<i; for(j=0;j<Num;j+= 2*l ) /*一组蝶形运算*/ for(k=0;k<l;k+) /*一个蝶形运算*/ mul(xj+k+l,WNum*k/2/l,&product); add(xj+k,product,&up); sub(xj+k,product,&down); xj+k=up; xj+k+l=down; c. 初始化, 复数类全局数组Wkvoid initW() int i; W=(complex *)malloc(sizeof(complex) * Num); for(i=0;i<Num;i+) Wi.real=

8、cos(2*PI/Num*i); Wi.img=-1*sin(2*PI/Num*i); 3. BTWZ滤波器和双线性Z转化为数字滤波器的实现设计巴特沃兹低通模拟滤波器时，我严格按照上学期的方法进行思考，思路是先做出2阶的滤波器，如果阶次高的话，我就把二阶滤波器做级联完成高阶滤波器的滤波。以下是二阶滤波器实现过程。完成BTWZ滤波器的系数计算void BTWZ(float fp, float fs,float ap, float as)int Ni,i;double N,P1000,Z3,z3,Wp;Wp=2*PI*fp;N=log(sqrt(pow(10,(as/10)-1)/(pow(10,

9、(ap/10)-1)/log(2);if(N=(double)Ni)Ni=(int)N;elseNi=(int)N+1;if(Ni%2=0)for(i=0;i<(Ni/2);i+)Pi=2*cos(2*i+Ni-1)*PI/(2*N);elseP0=0;for(i=0;i<(Ni-1)/2);i+)Pi+1=2*cos(2*(i+1)+Ni-1)*PI/(2*N);for(i=0;i<Ni;i+)printf("%fn",Pi);printf("%d",Ni);return P1000;通过z变换完成二阶滤波器从模拟到数字的转变void

10、 SHUANGXIANXINGZ(double Wp,double P,double in1000,double out1000)int i;out1000=0;in-1=in-2=out-1=out-2=0;a1=Wp*Wp;a2=2*a1;a3=a1;b1=1-P*Wp+a1；b2=2*a1-2;b3=a1+P*Wp+1;for(i=0;i<1000;i+)outi=ini*a1/b1+ini-1*a2/b1+ini-2*a3/b1-outi-1*b2/b1-outi-2*b3/b1;return out1000;4. 调用MATLAB画图的实现#ifdef USE_MATLAB/s

11、endIntToMatlab(agTemp1, 1024, "temp1");sendFloatToMatlab (agTemp3, 1024, "temp4");sendFloatToMatlab (Amp, 1024, "temp3");sendFloatToMatlab(HPs, 1024, "temp2");sendFloatToMatlab (Amp1, 1024, "temp1");sendFloatToMatlab(HPs1, 1024, "temp5");se

12、ndFloatToMatlab(agTemp4, 1024, "temp6");sendFloatToMatlab(Amp2, 1024, "temp7");sendFloatToMatlab(HPs2, 1024, "temp8");matlabEval("clf");matlabEval("subplot(6, 1, 1)");matlabEval("plot(0:length(temp2)-1, temp2, 'r', 0:length(temp2)-1, tem

13、p2, 'g');, axis(0 length(temp2)-1 min(temp2) max(temp2)");matlabEval("grid on");matlabEval("xlabel('sample');, ylabel('amplitude');");matlabEval("title('test');, legend('temp1', 'temp2');");matlabEval("subplot(6,

14、 1, 2)");matlabEval("plot(0:length(temp3)-1, temp3, 'r');, axis(0 length(temp3)-1 0 100)");matlabEval("set(gca,'XTickLabel',str2num(get(gca,'XTickLabel')/1024*360)");/重新设置x横坐标刻度值为原坐标数值/1024个点*频率范围matlabEval("grid on");matlabEval("xlabel

15、('sample');, ylabel('amplitude');");matlabEval("title('test2');");matlabEval("subplot(6, 1, 3)");/一列六个图中的第三个图表示原信号时域图形matlabEval("plot(0:length(temp4)-1, temp4, 'r');, axis(0 length(temp4)-1 min(temp4) max(temp4)");matlabEval("g

16、rid on");matlabEval("xlabel('sample');, ylabel('amplitude');");matlabEval("title('test2');");matlabEval("subplot(6, 1, 4)");matlabEval("plot(0:length(temp1)-1, temp1, 'r');, axis(0 length(temp1)-1 0 100)");matlabEval("

17、set(gca,'XTickLabel',str2num(get(gca,'XTickLabel')/1024*360)");/重新设置x横坐标刻度值为原坐标数值/1024个点*频率范围matlabEval("grid on");matlabEval("xlabel('sample');, ylabel('amplitude');");matlabEval("title('test2');");matlabEval("subplot(6

18、, 1, 5)");matlabEval("plot(0:length(temp5)-1, temp5, 'r');, axis(0 length(temp5)-1 min(temp5) max(temp5)");matlabEval("set(gca,'XTickLabel',str2num(get(gca,'XTickLabel')/1024*360)");/重新设置x横坐标刻度值为原坐标数值/1024个点*频率范围matlabEval("grid on");matlabE

19、val("xlabel('sample');, ylabel('amplitude');");matlabEval("title('test2');"); matlabEval("subplot(6, 1, 6)");matlabEval("plot(0:length(temp6)-1, temp6, 'r');, axis(0 length(temp6)-1 min(temp6) max(temp6)");matlabEval("set(g

20、ca,'XTickLabel',str2num(get(gca,'XTickLabel')/1024*360)");/重新设置x横坐标刻度值为原坐标数值/1024个点*频率范围matlabEval("grid on");matlabEval("xlabel('sample');, ylabel('amplitude');");matlabEval("title('test2');"); matlabEval("subplot(6, 2,

21、 1)");matlabEval("plot(0:length(temp7)-1, temp7, 'r');, axis(0 length(temp7)-1 min(temp7) 100)");matlabEval("set(gca,'XTickLabel',str2num(get(gca,'XTickLabel')/1024*360)");/重新设置x横坐标刻度值为原坐标数值/1024个点*频率范围matlabEval("grid on");matlabEval("

22、xlabel('sample');, ylabel('amplitude');");matlabEval("title('test2');"); matlabEval("subplot(6, 2, 2)");matlabEval("plot(0:length(temp8)-1, temp8, 'r');, axis(0 length(temp8)-1 min(temp8) max(temp8)");matlabEval("set(gca,'XT

23、ickLabel',str2num(get(gca,'XTickLabel')/1024*360)");/重新设置x横坐标刻度值为原坐标数值/1024个点*频率范围matlabEval("grid on");matlabEval("xlabel('sample');, ylabel('amplitude');");matlabEval("title('test2');");matlabEval("sqrt(8)");#endif将低通

24、滤波器转化为高通滤波器的实现void HIGH(float Wp, float Ws, float Ap, float As)float OMGp,OMGs,ys,LMDs; int Ni,i;OMGp=tan(Wp/2); /技术指标归一OMGs=tan(Ws/2); /ys=OMGs/OMGp;LMDs=1/ys;double N,P1000,Z3,z3,Wp;Wp=2*PI*fp;N=log(sqrt(pow(10,(As/10)-1)/(pow(10,(Ap/10)-1)/log(LMDs)/2;if(N=(double)Ni)Ni=(int)N;elseNi=(int)N+1;if(

25、Ni%2=0)for(i=0;i<(Ni/2);i+)Pi=2*cos(2*i+Ni-1)*PI/(2*N);elseP0=0;for(i=0;i<(Ni-1)/2);i+)Pi+1=2*cos(2*(i+1)+Ni-1)*PI/(2*N);for(i=0;i<Ni;i+)printf("%fn",Pi);printf("%d",Ni);return P1000;void SHUANGXIANXINGZ(double Wp,double P,double in1000,double out1000)int i;out1000=0;in-

26、1=in-2=out-1=out-2=0;a1=1;a2=-2;a3=1;b1=Wp*Wp-P*Wp+1;b2=2*Wp*Wp-2;b3=Wp*Wp+P*Wp+1;for(i=0;i<1000;i+)outi=ini*a1/b1+ini-1*a2/b1+ini-2*a3/b1-outi-1*b2/b1-outi-2*b3/b1;return out1000;四、实验结果及分析通过低通滤波后波形由于本实验数据本来就是低频信号，所以又做了一个高通滤波器来检验对错。通过高通滤波后波形由图可得，低频信号被滤掉五、实验感想在之前的学习中，我学习过FFT，巴特沃兹滤波器，双线性Z变换这些东西，但是并

27、不知道在实际应用中应该如何应用，通过张新峰老师的实验，我知道了课本上的东西是很重要的，而之前欠下的时间在未来也会加倍补回来。然后就是实践在我这个专业中时非常重要的，只有通过实践才能对学习的知识有更深刻的理解。比如通过这次试验，我就基本理解了设计滤波器的步骤和方法，如果再让我做一个类似的滤波器，我就会有思路了。然后谢谢张新峰老师的教导，作为对这个实验一无所知的学生到能够自己编出来东西，我自己的水平也得到了很大的提升。参考文献：C语言程序设计（第二版） 何钦铭 颜晖数字图像处理导论（第二版） 胡广书附录：源程序#include <stdio.h>#include <stdlib.

28、h>#include <malloc.h>#include <memory.h>#include <math.h>/#include "MSP.H"#include "matlab.h"#include "math.h"#define USE_MATLAB#define N 1024#define PI 3.1415926typedef struct double real; double imag;complex;complex xN,*W,yN,AN;float Wp,P1024;floa

29、t out1024;float input1024,output1024;int ln,l,NN;float a20,a1120,b20;void add(complex a,complex b,complex *c)/定义复数加法实部与实部相加，虚部与虚部相加 c->real=a.real+b.real; c->imag=a.imag+b.imag;void mul(complex a,complex b,complex *c)/定义复数乘法 c->real=a.real*b.real - a.imag*b.imag; c->imag=a.real*b.imag +

30、a.imag*b.real;void sub(complex a,complex b,complex *c)/定义复数减法 c->real=a.real-b.real; c->imag=a.imag-b.imag;void divi(complex a,complex b,complex *c)/定义复数除法 c->real=( a.real*b.real+a.imag*b.imag )/( b.real*b.real+b.imag*b.imag); c->imag=( a.imag*b.real-a.real*b.imag)/(b.real*b.real+b.imag

31、*b.imag);/*变址计算，将x(n)码位倒置*/void change() complex temp; unsigned short i=0,j=0,k=0; double t; for(i=0;i<N;i+) k=i;j=0; t=(log(N)/log(2); while( (t-)>0 ) j=j<<1; j|=(k & 1); k=k>>1; if(j>i) temp=xi; xi=xj; xj=temp; void fft() int i=0,j=0,k=0,l=0; complex up,down,product; chang

32、e(); for(i=0;i< log(N)/log(2) ;i+) /*一级蝶形运算*/ l=1<<i; for(j=0;j<N;j+= 2*l ) /*一组蝶形运算*/ for(k=0;k<l;k+) /*一个蝶形运算*/ mul(xj+k+l,WN*k/2/l,&product); add(xj+k,product,&up); sub(xj+k,product,&down); xj+k=up; xj+k+l=down; /*初始化变换核*/void initW() int i; W=(complex *)malloc(sizeof(c

33、omplex) * N); for(i=0;i<N;i+) Wi.real=cos(2*PI/N*i); Wi.imag=-1*sin(2*PI/N*i); void BTWZ(float fp, float fs,float ap, float as)int Ni,i;float n,OMGp,OMGs,ys,LMDs;Wp=2*PI*fp;OMGp=tan(PI*fp); /技术指标归一OMGs=tan(PI*fs); /ys=OMGs/OMGp;LMDs=1/ys;n=log(sqrt(pow(10,(as/10)-1)/(pow(10,(ap/10)-1)/log(LMDs)/2

34、;if(n=(float)Ni)Ni=(int)n;elseNi=(int)n+1;if(Ni%2=0)for(i=0;i<(Ni/2);i+)Pi=2*cos(2*i+Ni-1)*PI/(2*n);elseP0=0;for(i=0;i<(Ni-1)/2);i+)Pi+1=2*cos(2*(i+1)+Ni-1)*PI/(2*n);NN=Ni;void Z(float PP,float in1024)int i;float a1,a2,a3,b1,b2,b3;in-1=in-2=out-1=out-2=0;a1=1;a2=-2;a3=1;b1=Wp*Wp-PP*Wp+1;b2=2*W

35、p*Wp-2;b3=Wp*Wp+PP*Wp+1;for(i=0;i<1024;i+)outi=ini*a1/b1+ini-1*a2/b1+ini-2*a3/b1-outi-1*b2/b1-outi-2*b3/b1;void main()FILE *p; /定义文件型变量int agTemp11024;float agTemp21024+2;float agTemp31024,agTemp41024;float Amp1024,Amp11024,Amp21024;/Amp为滤波后信号幅频，Amp1为原信号幅频。float HPs1024,HPs11024,HPs21024;float *p

36、p3DFigure;int e,i,j,debug_flag,jj;p=fopen("105_heart1.txt","r");/将文件105_heart.txt中的数据读入到程序中for(i=0;i<1024+2;i+) /取1024个数作分析，注意前两个数据不要，将文件中的数据存放到agTemp2数组中。fscanf(p,"%f",&agTemp2i);fclose(p);for(i=0;i<1024;i+)/去除最开始的2个数据，保存在agTemp3内 agTemp3i=agTemp2i+2;for(i=0;

37、i<N;i+) xi.real=agTemp3i;xi.imag=0;/取前1024个数据，实部和虚部分别储存在 dataR 与 dataI 中initW();fft();/做FFTfor(i=0;i<N;i+) Ampi=sqrt(xi.real*xi.real+xi.imag*xi.imag); BTWZ(2,350,3,30);if(NN%2=0)Z(P0,agTemp3);for(i=1;i<NN/2;i+)Z(Pi,out);elseZ(P0,agTemp3);for(i=1;i<(NN-1)/2;i+)Z(Pi,out);for(i=0;i<N;i+)

38、 xi.real=outi;xi.imag=0;/取前1024个数据，实部和虚部分别储存在 dataR 与 dataI 中initW();fft();/做FFTfor(i=0;i<N;i+) Amp1i=sqrt(xi.real*xi.real+xi.imag*xi.imag); /画图#ifdef USE_MATLABstartMatlab("t");#endifpp3DFigure = (float *)calloc(128, sizeof(float *);/在内存的动态存储区中分配N个长度为size的连续空间/sizeof计算float*所占的字节数for (

39、i = 0; i < 128; i+)pp3DFigurei = (float *)calloc(64, sizeof(float);for (i = 0; i < 128; i+)for (j = 0; j < 64; j+)pp3DFigureij = (float)i + j;#ifdef USE_MATLABsendFloatToMatlab (agTemp3, 1024, "temp2");sendFloatToMatlab (Amp, 1024, "temp3");sendFloatToMatlab(out, 1024, &

40、quot;temp4");sendFloatToMatlab (Amp1, 1024, "temp1");matlabEval("clf");matlabEval("subplot(2, 2, 1)");matlabEval("plot(0:length(temp2)-1, temp2, 'r', 0:length(temp2)-1, temp2, 'g');, axis(0 length(temp2)-1 min(temp2) max(temp2)");matlabEva

41、l("grid on");matlabEval("xlabel('sample');, ylabel('amplitude');");matlabEval("title('test');, legend('temp1', 'temp2');");matlabEval("subplot(2, 2, 2)");matlabEval("plot(0:length(temp3)-1, temp3, 'r');, axi

42、s(0 length(temp3)-1 0 100)");matlabEval("set(gca,'XTickLabel',str2num(get(gca,'XTickLabel')/1024*360)");/重新设置x横坐标刻度值为原坐标数值/1024个点*频率范围matlabEval("grid on");matlabEval("xlabel('sample');, ylabel('amplitude');");matlabEval("title(

43、'test2');");matlabEval("subplot(2, 2, 3)");/一列六个图中的第三个图表示原信号时域图形matlabEval("plot(0:length(temp4)-1, temp4, 'r');, axis(0 length(temp4)-1 min(temp4) max(temp4)");matlabEval("grid on");matlabEval("xlabel('sample');, ylabel('amplitude');");matlabEval("title('test2');")