[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);
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型低通滤波器的幅频响应曲线图
《生物医学信号处理》实习报告
学生姓名：
学号：
实验室名称：
项目名称：心电信号的预处理
项目内容：
1)阅读相关文献资料，理解噪声的相关知识，据此综述抑制信号中干扰噪声信号的基本理论和基本方法；
2)设计两种分别用于抑制不同噪声的滤波器；
3)基于W-H方程，设计最优滤波器；
4)运用前面设计的几种滤波器，对加有噪声的模拟ECG信号进行去噪；
(5)基线漂移和呼吸时ECG幅值的变化
基线漂移和呼吸时ECG幅值的变化一般由人体呼吸!电极移动等低频干扰所引起,频率小于5Hz;其变化可视为一个加在心电信号上的与呼吸频率同频率的正弦分量,在0.015一0.3Hz处基线变化变化幅度的为ECG峰峰值的15%"。
上面的电极接触噪声与人为运动所产生的噪声是人为因素造成的，当然也可以通过人为因素来避免。然而工频干扰、肌电干扰(EMG)与基线漂移和呼吸时ECG幅值的变化就不是人为因素所能消除的了。
图二
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;%阻带内的最小衰减
%数字角频率转换成模拟角频率
[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');
legend('原始信号','加噪后信号','维纳滤波后信号','误差');
图四
4.分别运用上面设计的滤波器和最优滤波器对含高斯噪声心律失常ECG数据消噪，并运用SNR指标定量分析滤波器的去噪能力
1）设计SNR（信噪比）的计算方法
functionsnr=SNR(I,In)
%计算信噪比函数
%I:original signal
fork=1:L2
Rww(k)=Rw(k&#os(1,N);
fork=1:N
Rx(k)=Rss(k)+Rww(k);
end
Rxx=zeros(N,N);
fori=1:N
forj=1:N
if(i<=j)
Rxx(i,j)=Rx(j-i+1);
else
Rxx(i,j)=Rx(i-j+1);
图三
3.维纳滤波器（Winner.m）
function[H,Emin]=Winner(Rss,Rww,N)
L1=(length(Rs)+1)/2;
Rss=zeros(1,L1);
fork=1:L1
Rss(k)=Rs(k+L1-1);
end
L2=(length(Rw)+1)/2;
Rww=zeros(1,L2);
%In:noisy signal(ie.original signal+noise signal)
snr=0;
Ps=sum(abs(I).^2);%signal power
Pn=sum(abs(In-I).^2);%noise power
snr=10*log10(Ps/Pn);
%其中I是纯信号，In是带噪声的信号，snr是信噪比
下面是对维纳滤波器的设计仿真：
维纳滤波器是以均方误差最小（LMS）为准则的，它根据过去观测值和当前观测值来估计信号的当前值，因此它的解形式是系统的传递函数或单位脉冲响应
图一
设计维纳滤波器的过程就是寻求在最小均方误差下滤波器的单位脉冲响应或传递函数的表达式，其实质就是解维纳-霍夫（Wiener-Hopf）方程。
第Ⅰ类切比雪夫滤波器的幅度平方函数定义为
（3）
其中， 为切比雪夫多项式（ ）， 为滤波器阶数；
为通带截止角频率，此处是指被通带波纹所限制的最高角频率，；
为小于1的正数，表示通带内幅度波动的程度， 越小，幅度波动越小；
其特点为：
1)当 时， 在 之间呈等波纹变化， 越小，波动幅度越小；
2)无论 为何值，所有的曲线都通过 点，该点被定义为 点；
2)巴特沃兹低通滤波器(snrbutter.m)
clc;
clear;
Test_Signal_1=load('118e00m.mat');
%导入数据
Signal_1=Test_Signal_1.val(1,:);
Fs=2*max(abs(Signal_1));
%采样频率确定
Signal_1=Signal_1/(Fs/2);
均方误差为： (5)
维纳—霍夫方程最小均方误差是下的解为： (6)
其中， 为滤波器的系数向量， 为含有噪声的混合信号的自相关矩阵， 为混合信号和原始信号的互相关向量，为此先求出混合信号的自相关函数以及混合信号和原始信号的互相关函数，这两个函数，我们需要有样本得到。
编写的源程序：
1.巴特沃斯滤波器(Butter.m)
为了滤除掉上述三种噪声，我按照实验要求设计了三种不同的滤波器。分别是巴特沃斯滤波器与切比雪夫滤波器。为了对比他们的滤波效果，又设计了一个维纳滤波器。最后运用SNR指标定量分析了不同滤波器的去噪能力。以下是3种滤波器的原理：
1.巴特沃斯滤波器的设计原理
其特点是通频带内的频率响应曲线最大限度平坦，没有起伏，而在阻频带则逐渐下降为零（对理想低通滤波的逼近：巴特沃思滤波器是以原点附近的最大平坦响应来逼近理想低通滤波器）。而滤波器的幅频特性是随着滤波器的阶次N的增加而变得越来越好，在截止频率 有：
5)总结不同滤波器的去噪效果。
原理（写出具体的计算公式）
人体心电信号是一种弱电信号,信噪比低"一般正常的心电信号频率范围为0.05一100Hz,幅度为10产v一smV"而这个心动周期信号带宽主要集中在0~58士1白Hz,P波的带宽为0一55士19Hz,T波的带宽为0~11士ZHz"而90%的心电信号(ECG)频谱能量集中在0.25/35Hz之间"采集心电信号时,一会受到各种噪声的干扰,噪声来源通常有下面几种:
(1)工频干扰
50Hz工频干扰是由人体的分布电容所引起,工频干扰的模型由50Hz的正弦信号及其谐波组成"由于本实验中所采用的是MIT一BIH数据库的心电信号,所有的工频干扰都为60Hz"幅值通常与ECG峰峰值相当或更强"。
(2)电极接触噪声
电极接触噪声是瞬时干扰,来源于电极与肌肤的不良接触,即病人与检测系统的连接不好"其连接不好可能是瞬时的,如病人的运动和振动导致松动;也可能是检测系统不断的开关!放大器输入端连接不好等"电极接触噪声可抽象为快速!随机变化的阶跃信号,它按指数形式衰减到基线值,包含工频成分"这种瞬态过渡过程可发生一次或多次!其特征值包括初始瞬态的幅值和工频成分的幅值!衰减的时间常数;其持续时间一般为15左右,幅值可达记录仪的最大值"。
3)当 时，若 为奇数， ；若 为偶数， （4）；
4)当 时，曲线呈单调下降；
5)通带内的起伏使对应的相频特性呈现非线性；
3.维纳滤波器的设计原理
维纳滤波器属于现代滤波器，传统的滤波器只能滤除信号和干扰频带没有重叠的情况，当信号和干扰频带有重叠的时候传统滤波器将无能为力，这时就需要用到现代滤波器，现代滤波器利用信号和干扰的统计特征（如自相关函数、功率谱等）导出一套最佳估计算法，然后用硬件或软件予以实现。
wp=fp/(Fs/2);
ws=fs/(Fs/2);
%数字转模拟
wp = tan(pi*wp/2);
ws = tan(pi*ws/2);
[N,wc]=buttord(wp,ws,Rp,Rs,'s');
[B,A]=butter(N,wc,'s');
[b,a]=bilinear(B,A,.5);%双线性变换法
（1）
衰减具有 不变性。通带、阻带均具有单调下降的特性。
低通巴特沃斯滤波器的幅频响应为：
（2）
其中， 为滤波器的阶数，取正整数；
为滤波器的截止频率，为其 的频率；
为滤波器的通带截止频率；
2.切比雪夫滤波器的设计原理
切比雪夫滤波器有两类，一类是在同带内幅度频率响应呈现等波纹特性，而阻带内是单调的全极点滤波器；另一类是在通带内幅度频率响应呈现单调特性，而阻带内是等波纹特性，同时有零点和极点的滤波器，这类滤波器的零点位于s平面的虚轴上。