北工大通信系統(tǒng)實驗

1、信號處理與分析課程設(shè)計訓(xùn)練報告- ECG信號高頻噪聲消除學(xué)號：姓名：指導(dǎo)老師：張新峰完成日期：2016年6月5日緒論本實驗是信號處理相關(guān)方法在心電圖ECG（electocardiogram）信號處理中的一個具體應(yīng)用。通過心電圖進(jìn)行疾病診斷是心血管疾病診斷的一種重要方法。該方法可靠簡便，在臨床中應(yīng)用很廣泛。近年來給予數(shù)字信號處理和模式識別的ECG信號自動分析成為一個研究這點。ECG信號在采集過程中收到工頻信號、肌電信號以及基線漂移等干擾，由此帶來的噪聲對自動分析的結(jié)果影響很大。因此在對ECG信號分析之前，必須進(jìn)行噪聲消除研究。本實驗內(nèi)容給定ECG信號，合理設(shè)計濾波器，消除高頻噪聲的影響。一、實驗

2、內(nèi)容心電信號是微弱低頻人體生理電信號，通常頻率在0.05100Hz，幅值不超過4mv，通過安裝在皮膚表面的電極來獲取。由于人體是一個復(fù)雜的生 命系統(tǒng)，存在50H 工頻干擾及基線漂移等其他生理電信號的干擾。噪聲可能會影響到醫(yī)生的臨床診斷，因此，需對心電信號進(jìn)行濾波，即必須做好前端數(shù)據(jù)采集的軟硬件設(shè)計以保證心電數(shù)據(jù)的可靠和準(zhǔn)確。（本次實驗的心電數(shù)據(jù)是老師提供，為txt文檔）傳統(tǒng)醫(yī)療設(shè)備分別采用50Hz 帶阻濾波器和RC 高通濾波器濾除工頻干擾和基線漂移。但帶阻濾波器電路復(fù)雜，其特性對元器件的精度敏感，而基線漂移本質(zhì)上是一種緩慢變化的低頻信號，采用RC 濾波器很難將高通濾波器的過渡帶做得十分陡峭，基

3、線漂移補(bǔ)償效果不理想。因此，模擬方法往往不太容易獲得很好的特性。在本實驗中，我們以C語言為基礎(chǔ)，來實現(xiàn)高頻噪聲的消除。二、實驗原理1. 低通巴特沃茲模擬濾波器的實現(xiàn)1）概述巴特沃茲濾波器的特點是通頻帶內(nèi)的頻率響應(yīng)曲線最大限度平坦，沒有起伏，而在阻頻帶則逐漸下降為零。 在振幅的對數(shù)對角頻率的波得圖上，從某一邊界角頻率開始，振幅隨著角頻率的增加而逐步減少，趨向負(fù)無窮大。 一階巴特沃茲濾波器的衰減率為每倍頻6分貝，每十倍頻20分貝。二階巴特沃茲濾波器的衰減率為每倍頻12分貝、 三階巴特沃茲濾波器的衰減率為每倍頻18分貝、如此類推。巴特沃茲濾波器的振幅對角頻率單調(diào)下降，并且也是唯一的無論階數(shù)，振幅對角

4、頻率曲線都保持同樣的 形狀的濾波器。只不過濾波器階數(shù)越高，在阻頻帶振幅衰減速度越快。其他濾波器高階的振幅對角頻率圖和低級數(shù)的振幅對角頻率有不同的形狀。2）設(shè)計巴特沃茲模擬低通濾波器1確定參數(shù):，2求濾波器階次N： = N= 4確定系統(tǒng)函數(shù)N為奇數(shù)時： N為偶數(shù)時： 雙線性Z法把模擬濾波器轉(zhuǎn)換成數(shù)字濾波器3. FFT的算法按時間抽選（DIT）的基-2 FFT算法。當(dāng)輸入序列點數(shù)為 （L為整數(shù)），先將序列x(n)按n的奇偶分成以下兩組：則可以將DFT化為： k=0,1,2, 其中 分別為 點DFT。又根據(jù)系數(shù)的周期性可得到X(k)的完整表達(dá)：三、實現(xiàn)方法流程圖1. 讀取心電信號數(shù)據(jù)的實現(xiàn)float

5、 agTemp21024+2; /定義一個數(shù)組用來存放數(shù)據(jù)FILE *p; /定義指針p=fopen("119_heart1.txt","r"); / 打開心電圖數(shù)據(jù)文件for(i=0;i<1024+2;i+) / 通過循環(huán)放入數(shù)組中fscanf(p,"%f",&agTemp2i);fclose(p);2. FFT變換的實現(xiàn)a. 倒位序處理，先將xn進(jìn)行倒位序后，才能進(jìn)行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. 進(jìn)行蝶形運(yùn)算，將保存在全局復(fù)數(shù)數(shù)組xn中的時域信號轉(zhuǎn)為頻域信號: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+) /*一級蝶形運(yùn)算*/ l=1<

7、<i; for(j=0;j<Num;j+= 2*l ) /*一組蝶形運(yùn)算*/ for(k=0;k<l;k+) /*一個蝶形運(yùn)算*/ 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. 初始化, 復(fù)數(shù)類全局?jǐn)?shù)組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轉(zhuǎn)化為數(shù)字濾波器的實現(xiàn)設(shè)計巴特沃茲低通模擬濾波器時，我嚴(yán)格按照上學(xué)期的方法進(jìn)行思考，思路是先做出2階的濾波器，如果階次高的話，我就把二階濾波器做級聯(lián)完成高階濾波器的濾波。以下是二階濾波器實現(xiàn)過程。完成BTWZ濾波器的系數(shù)計算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變換完成二階濾波器從模擬到數(shù)字的轉(zhuǎn)變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. 調(diào)用MATLAB畫圖的實現(xiàn)#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)");/重新設(shè)置x橫坐標(biāo)刻度值為原坐標(biāo)數(shù)值/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)");/重新設(shè)置x橫坐標(biāo)刻度值為原坐標(biāo)數(shù)值/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)");/重新設(shè)置x橫坐標(biāo)刻度值為原坐標(biāo)數(shù)值/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)");/重新設(shè)置x橫坐標(biāo)刻度值為原坐標(biāo)數(shù)值/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)");/重新設(shè)置x橫坐標(biāo)刻度值為原坐標(biāo)數(shù)值/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)");/重新設(shè)置x橫坐標(biāo)刻度值為原坐標(biāo)數(shù)值/1024個點*頻率范圍matlabEval("grid on");matlabEval("xlabel('sample');, ylabel('amplitude');");matlabEval("title('test2');");matlabEval("sqrt(8)");#endif將低通

24、濾波器轉(zhuǎn)化為高通濾波器的實現(xiàn)void HIGH(float Wp, float Ws, float Ap, float As)float OMGp,OMGs,ys,LMDs; int Ni,i;OMGp=tan(Wp/2); /技術(shù)指標(biāo)歸一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;四、實驗結(jié)果及分析通過低通濾波后波形由于本實驗數(shù)據(jù)本來就是低頻信號，所以又做了一個高通濾波器來檢驗對錯。通過高通濾波后波形由圖可得，低頻信號被濾掉五、實驗感想在之前的學(xué)習(xí)中，我學(xué)習(xí)過FFT，巴特沃茲濾波器，雙線性Z變換這些東西，但是并

27、不知道在實際應(yīng)用中應(yīng)該如何應(yīng)用，通過張新峰老師的實驗，我知道了課本上的東西是很重要的，而之前欠下的時間在未來也會加倍補(bǔ)回來。然后就是實踐在我這個專業(yè)中時非常重要的，只有通過實踐才能對學(xué)習(xí)的知識有更深刻的理解。比如通過這次試驗，我就基本理解了設(shè)計濾波器的步驟和方法，如果再讓我做一個類似的濾波器，我就會有思路了。然后謝謝張新峰老師的教導(dǎo)，作為對這個實驗一無所知的學(xué)生到能夠自己編出來東西，我自己的水平也得到了很大的提升。參考文獻(xiàn)：C語言程序設(shè)計（第二版） 何欽銘 顏暉數(shù)字圖像處理導(dǎo)論（第二版） 胡廣書附錄：源程序#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)/定義復(fù)數(shù)加法實部與實部相加，虛部與虛部相加 c->real=a.real+b.real; c->imag=a.imag+b.imag;void mul(complex a,complex b,complex *c)/定義復(fù)數(shù)乘法 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)/定義復(fù)數(shù)減法 c->real=a.real-b.real; c->imag=a.imag-b.imag;void divi(complex a,complex b,complex *c)/定義復(fù)數(shù)除法 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+) /*一級蝶形運(yùn)算*/ l=1<<i; for(j=0;j<N;j+= 2*l ) /*一組蝶形運(yùn)算*/ for(k=0;k<l;k+) /*一個蝶形運(yùn)算*/ 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); /技術(shù)指標(biāo)歸一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中的數(shù)據(jù)讀入到程序中for(i=0;i<1024+2;i+) /取1024個數(shù)作分析，注意前兩個數(shù)據(jù)不要，將文件中的數(shù)據(jù)存放到agTemp2數(shù)組中。fscanf(p,"%f",&agTemp2i);fclose(p);for(i=0;i<1024;i+)/去除最開始的2個數(shù)據(jù)，保存在agTemp3內(nèi) agTemp3i=agTemp2i+2;for(i=0;

37、i<N;i+) xi.real=agTemp3i;xi.imag=0;/取前1024個數(shù)據(jù)，實部和虛部分別儲存在 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個數(shù)據(jù)，實部和虛部分別儲存在 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èi)存的動態(tài)存儲區(qū)中分配N個長度為size的連續(xù)空間/sizeof計算float*所占的字節(jié)數(shù)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)");/重新設(shè)置x橫坐標(biāo)刻度值為原坐標(biāo)數(shù)值/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');")