生物医学信号处理实习报告模版.doc

《生物医学信号处理》实习报告一姓名：学号：实验室名称：实习名称：心电信号的预处理内容：1) 根据相关文献资料，综述信号消噪的基本理论和基本方法；2) 设计两种分别用于抑制不同噪声的滤波器；3) 基于W-H方程，设计最优滤波器；4) 运用前面设计的几种滤波器，对加有噪声的模拟ECG信号进行去噪；5) 总结不同滤波器的去噪效果；滤波器设计原理一：两种能消噪的滤波器的设计原理1. 巴特沃斯滤波器的设计原理其特点是通频带内的频率响应曲线最大限度平坦，没有起伏，而在阻频带则逐渐下降为零（对理想低通滤波的逼近：巴特沃思滤波器是以原点附近的最大平坦响应来逼近理想低通滤波器）如下图一表示，可以看出滤波器的幅频特性是随着滤波器的阶次N的增加而变得越来越好，在截止频率有：（1）衰减具有不变性通带、阻带均具有单调下降的特性图一低通巴特沃斯滤波器的幅频响应为： （2）其中，为滤波器的阶数，取正整数； 为滤波器的截止频率，为其的频率； 为滤波器的通带截止频率；2. 切比雪夫滤波器的设计原理切比雪夫滤波器有两类，一类是在同带内幅度频率响应呈现等波纹特性，而阻带内是单调的全极点滤波器；另一类是在通带内幅度频率响应呈现单调特性，而阻带内是等波纹特性，同时有零点和极点的滤波器，这类滤波器的零点位于s平面的虚轴上。

如图二，为第Ⅰ类切比雪夫滤波器的幅频特性曲线图二第Ⅰ类切比雪夫滤波器的幅度平方函数定义为（3）其中，为切比雪夫多项式（），为滤波器阶数；为通带截止角频率，此处是指被通带波纹所限制的最高角频率， ；为小于1的正数，表示通带内幅度波动的程度，越小，幅度波动越小；其特点为：1) 当时，在之间呈等波纹变化，越小，波动幅度越小；2) 无论为何值，所有的曲线都通过点，该点被定义为点；3) 当时，若为奇数，；若为偶数， （4）；4) 当时，曲线呈单调下降；5) 通带内的起伏使对应的相频特性呈现非线性；3. 维纳滤波器的设计原理维纳滤波器属于现代滤波器，传统的滤波器只能滤除信号和干扰频带没有重叠的情况，当信号和干扰频带有重叠的时候传统滤波器将无能为力，这时就需要用到现代滤波器，现代滤波器利用信号和干扰的统计特征（如自相关函数、功率谱等）导出一套最佳估计算法，然后用硬件或软件予以实现下面是对维纳滤波器的设计仿真：维纳滤波器是以均方误差最小（LMS）为准则的，它根据过去观测值和当前观测值来估计信号的当前值，因此它的解形式是系统的传递函数或单位脉冲响应图三设计维纳滤波器的过程就是寻求在最小均方误差下滤波器的单位脉冲响应或传递函数的表达式，其实质就是解维纳-霍夫（Wiener-Hopf）方程。

均方误差为：(5)维纳—霍夫方程最小均方误差是下的解为：(6)其中，为滤波器的系数向量，为含有噪声的混合信号的自相关矩阵，为混合信号和原始信号的互相关向量，为此先求出混合信号的自相关函数以及混合信号和原始信号的互相关函数，这两个函数，我们需要有样本得到滤波器实现的代码1. 巴特沃斯滤波器（Butter.m） [B,A]=butter(3,2*pi*1000,'s'); %脉冲响应不变法低通滤波器 [num1,den1]=impinvar(B,A,4000); [h1,w]=freqz(num1,den1); [B,A]=butter(3,2/0.00025,'s'); %双线性变换法低通滤波器 [num2,den2]=bilinear(B,A,4000); [h2,w]=freqz(num2,den2); f=w/pi*2000; plot(f,abs(h1),'GREEN',f,abs(h2),'RED'); legend('脉冲响应不变法','双线性变换法'); grid; xlabel('频率/Hz ') ylabel('幅值/dB') title('巴特沃斯低通滤波器');图四2. 切比雪夫滤波器（qiebixuefu1.m） clc,clear %确定数字滤波器指标 Fs=10000; %采样频率 fp=280;fs=450; wp=2*pi*fp/Fs; %通带截止频率 Ap=0.1; %通带内的最大衰减 ws=2*pi*fs/Fs; %阻带截止频率 As=40; %阻带内的最小衰减 %数字角频率转换成模拟角频率 %双线性变换法 wps=2*Fs*tan(wp/2); wss=2*Fs*tan(ws/2); %设计模拟滤波器 k1=sqrt((10^(0.1*As)-1)/(10^(0.1*Ap)-1)); qsp=wss/wps; N=ceil(acosh(k1)/acosh(qsp)); e=sqrt(10^(0.1*Ap)-1); Ap=10*log(1+e^2); %计算得到模拟低通H(s)分子和分母的系数 [z,p,k]=cheb1ap(N,Ap); [bp,ap]=zp2tf(z,p,k); [bs,as]=lp2lp(bp,ap,Fs); %复变量映射s-->z [b,a]=bilinear(bs,as,Fs); %画出切比雪夫1型低通滤波器的幅频响应曲线图 [H,w] =freqs(bp,ap,2000);%按n指定的频率点 给出频率响应 magH2 = (abs(H)).^2; %给出传递函数的幅度平方 %输出波形 plot(w,magH2);title('切比雪夫1型低通滤波器的幅频响应特性 '); xlabel ('w/wc'); ylabel('Chebyshev 1 | H(jw)|^2');图五3. 维纳滤波器（Winner.m） function[H,Emin]=Winner(Rs,Rw,N) L1=(length(Rs)+1)/2; Rss=zeros(1,L1); for k=1:L1 Rss(k)=Rs(k+L1-1); end L2=(length(Rw)+1)/2; Rww=zeros(1,L2); for k=1:L2 Rww(k)=Rw(k+L2-1); end Rx=zeros(1,N); for k=1:N Rx(k)=Rss(k)+Rww(k); end Rxx=zeros(N,N); for i=1:N for j=1:N if(i<=j) Rxx(i,j)=Rx(j-i+1); else Rxx(i,j)=Rx(i-j+1); end end end H=pinv(Rxx)*Rss'; Emin=Rss(1)-sum(H*Rss);测试程序：（WinnerTest.m） clc; t=-50:0.1:50; w=0.1*pi; s=cos(w.*t); w=randn(1,length(t)); x=s+w; Rs=xcorr(s);% 信号的自相关函数 Rw=xcorr(w); % 白噪声的自相关函数 [h,E]=Winner(Rs,Rw,length(t));% 维纳滤波 se=filter(h,1,x); plot(t,s); hold on plot(t,x,'g'); hold on plot(t,se,'r'); e=s-se; plot(t,e,'BLACK'); legend('原始信号','加噪后信号','维纳滤波后信号','误差');图六4.分别运用上面设计的滤波器和最优滤波器分别对含噪声的模拟ECG信号进行去噪，并运用SNR指标定量分析滤波器的去噪能力：1) 读入原始ECG模拟信号（代码见ecgsyn.m和derivsecgsyn.m），并存入ECG.txt中；2) 提取ECG.txt中的（500-2500个信号），保存至ECG（2000）.txt中，并将生成的心电信号中加入高斯噪声，使其生成含噪的心电信号，保存至ECGWithGauss(2000).txt中（代码见ECG(2000).m），如图七；图七3) 设计SNR（信噪比）的计算方法 function snr=SNR(I,In) %计算信噪比函数 %I:original signal %In:noisy signal(ie.original signal+noise signal) snr=0; Ps=sum(sum((I-mean(mean(I))).^2));%signal power Pn=sum(sum((I-In).^2));%noise power snr=10*log10(Ps/Pn) %其中I是纯信号，In是带噪声的信号，snr是信噪比4) 利用所设计的三种滤波器对所生成的ECG信号除噪；a. 巴特沃斯滤波器消噪（SNRButter.m） clear,clc X=importdata('ECG(2000).txt');%加噪心电数据读入 X1=X(:); X2=importdata('ECGWithGauss(2000).txt'); X3=X2(:); [B,A]=butter(3,2*pi*1000,'s'); [num1,den1]=impinvar(B,A,4000); [h1,w]=freqz(num1,den1); f=w/pi*2000; y=filter(num1,den1,X3); figure; subplot(3,1,1); plot(X1); %原ECG信号 title('原始信号'); subplot(3,1,2); plot(X3,'g') title('加入高斯噪声后的ECG信号');%加入高斯噪声之后的噪声 subplot(3,1,3); plot(y,'r'); title('除噪后的ECG信号'); snr1=SNR(X1,X3); sne2=SNR(X1,y);b. 切比雪夫1型滤波器(SNRQie.m) clc,clear %确定数字滤波器指标 Fs=1000; %采样频率 fp=5;%通带带截止频率 fs=200;%阻带截止频率 wp=2*pi*fp/Fs; %通带截止频率 Ap=0.1; %通带内的最大衰减 ws=2*pi*fs/Fs; %阻带截止频率 As=40; %阻带内的最小衰减 %数字角频率转换成模拟角频率 %双线性变换法 wps=2*Fs*tan(wp/2); wss=2*Fs*tan(ws/2); %设计模拟滤波器 k1=sqrt((10^(0.1*As)-1)/(10^(0.1*Ap)-1)); qsp=wss/wps; N。