杭电_数字信号处理课程设计_实验5

注：心电信号文件(.txt)在附属资源中。
一、实验目的
1、探究心电信号的初步分析。心电信号(频率-般在0.05Hz ~100Hz范围)是一种基本的人体生理信号，体表检测人体心电信号中常带有工频干扰(50HZ)、 基线漂移(频 率低于0.5Hz)和肌电干扰等各种噪声。
2、为了得到不失真的原始心电信号，需要滤波预处理。设计数字低通滤波器、高通滤波器、带阻滤波器，用MATLAB软件对含噪心电信号分别进行高通、带阻和低通滤波等处理，将心电信号中的低频基线漂移、50Hz 工频高频和高频杂波进行滤除。
3、通过观察对含噪心电图信号的滤波作用，获得数字滤波的感性知识。

二、实验要求及内容
实验题目：
给定一组干净心电信号数据，数据文件存于C盘Ecg.txt。采样频率 Fs = 500Hz。
1、编写程序读出心电信号，并在屏幕上打印出其波形。
2、产生模拟高斯白噪声信号，与干净心电混合,设计一个IIR低通滤波器和一个FIR 低通滤波器分别滤除心电信号中的白噪声干扰，调整白噪声信噪比大小，对滤波前后的心 电信号的频谱进行分析比较。其中数字低通滤波器指标要求，通带截止频率Wp=0.1π，阻带截止频率Ws=0.16π，阻带衰减不小于15 dB，通带衰减不大于1 dB。
要求:编写一个IIR低通滤波器和一个FIR低通滤波器仿真程序，在屏幕上打印出数字滤波器的频率区间[0, π]上的幅频响应特性由线(H(e^jw)) ;计算其对含噪心电信号的 低通滤波响应序列，并在屏幕上打印出干净心电信号波形，含工频干扰的心电信号波形以及IIR低通滤波和FIR低通后的信号波形，并进行比较;同时对滤波前后的心电信号的频谱进行分析比较，并在屏幕上打印出滤波前后的心电信号的频谱，观察其变化。
3、产生模拟工频信号，与干净心电混合，设计一个带阻滤波器(50Hz 陷波器)滤除心电信号中的电源线干扰，调整工频幅度大小，对滤波前后的心电信号的频谱进行分析比 较。其中带阻滤波器指标要求，通带下限频率Wp1=0.18π， 阻带下截止频率Ws1=0.192 π，阻带上截止频率Ws2=0.208π，通带上限频率Wp2=0.22π，阻带衰减不小于15 dB, 通带衰减不大于1 dB。
要求:编写IIR带阻滤波器仿真程序,在屏幕上打印出数字滤波器的频率区间[0, π]上 的幅频响应特性由线(H(e^jw ));计算其对含工频干扰的心电信号的带阻滤波响应序列，并在屏幕上打印出干净心电信号波形，含工频干扰的心电信号波形以及滤波后的信号波形，并进行比较;同时对滤波前后的心电信号的频谱进行分析比较，并在屏幕上打印出滤波前后的心电信号的频谱，观察其变化。
4、产生模拟基线漂移信号，与干净心电信号混合，设计一个高通滤波器滤除心电信号中的基线低频干扰，调整基线的幅度大小,对滤波前后的心电信号的频谱进行分析比较。 其中，高通滤波器指标要求，通带截止频率Wp=0.0028π，阻带 截止频率Ws=0.0012π， 阻带衰减不小于15 dB，通带衰减不大于1 dB。
要求:编写IIR高通滤波器(或FIR高通滤波器)仿真程序，在屏幕上打印出数字滤 波器的频率区间[0,π]上的幅频响应特性由线(H(e^jw);计算其对含基线低频干扰的心电信号的高通滤波响应序列，并在屏幕上打印出干净心电信号波形，含基线低频干扰的心 电信号波形以及滤波后的信号波形，并进行比较;同时对滤波前后的心电信号的频谱进行 分析比较，并在屏幕上打印出滤波前后的心电信号的频谱，观察其变化。

实验代码：
一、

clear,clc;
val = importdata('Ecg.txt');
signal = val(1,1:1800);
fs = 500;
figure(1);
subplot(1,1,1);
plot(signal);
title('干净的ECG信号');
xlabel('采样点');
ylabel('幅值(dB)');
grid on;

二、

clear,clc;
val = importdata('Ecg.txt');
signal = val(1,1:1800);
Fs = 500;
figure(2);
subplot(4,1,1);
plot(signal);
title('干净的ECG信号');
xlabel('采样点');ylabel('幅值(dB)');
grid on;
XK1=fft(signal,1800);
magXK1=abs(XK1); %幅频特性
figure(3);
subplot(4,1,1);
k1=0:length(magXK1)-1;
stem(k1,magXK1,'.'); %信号幅频特性曲线
xlabel('k');
ylabel('|X(k)|');
title('干净的ECG信号频谱');

% 加入高斯白噪声信号
signal1 = awgn(signal,10,'measured');
figure(2);
subplot(4,1,2);
plot(signal1);
title('含噪声的ECG信号');
xlabel('采样点');ylabel('幅值(dB)');
grid on;
XK2=fft(signal1,1800);
magXK2=abs(XK2); %幅频特性
figure(3);
subplot(4,1,2);
k2=0:length(magXK2)-1;
stem(k2,magXK2,'.'); %信号幅频特性曲线
xlabel('k');
ylabel('|X(k)|');
title('含噪声的ECG信号频谱');

% IIR低通滤波器设计
wp=0.1*pi;
ws=0.16*pi;
Fp=2*Fs*tan(wp/2);
Fc=2*Fs*tan(ws/2);
Rp=1;
Rs=15;
[N,Wn] = buttord(Fp,Fc,Rp,Rs,'s');
[Z,P,K] = buttap(N);
[Bap,Aap] = zp2tf(Z,P,K);
[b,a] = lp2lp(Bap,Aap,Wn);
[bz,az] = bilinear(b,a,Fs);
[H1,W]=freqz(bz,az);
y1=filter(bz,az,signal1);
figure(1);
subplot(2,1,1);
plot(W/pi,20*log10(abs(H1)));
xlabel('\omega/\pi');
ylabel('幅度 (dB)');
title('IIR低通频谱');
grid on;
figure(2);
subplot(4,1,3);
plot(y1);
title('IIR滤波');
xlabel('采样点');ylabel('幅值(dB)');
grid on;
XK3=fft(y1,1800);
magXK3=abs(XK3); %幅频特性
figure(3);
subplot(4,1,3);
k3=0:length(magXK3)-1;
stem(k3,magXK3,'.'); %信号幅频特性曲线
xlabel('k');
ylabel('|X(k)|');
title('IIR低通滤波后的ECG信号频谱');

% FIR低通滤波器设计
wp=0.1*pi;
ws=0.16*pi;
wdelta=ws-wp;
N=ceil(8*pi/wdelta);
wn=(wp+ws)/2;
h=fir1(N-1,wn/pi,hanning(N));
[H2,W2]=freqz(h,1,512);
y2=conv(signal1, h);
figure(1);
subplot(2,1,2);
plot(W/pi,20*log10(abs(H2)));
xlabel('\omega/\pi');
ylabel('幅度 (dB)');
title('FIR低通频谱');
grid on;
figure(2);
subplot(4,1,4);
plot(y2);
title('FIR滤波')
xlabel('采样点');ylabel('幅值(dB)');
grid on;
XK4=fft(y2,1800);
magXK4=abs(XK4); %幅频特性
figure(3);
subplot(4,1,4);
k4=0:length(magXK4)-1;
stem(k4,magXK4,'.'); %信号幅频特性曲线
xlabel('k');
ylabel('|X(k)|');
title('FIR低通滤波后的ECG信号频谱');

三、

clear,clc;
val = importdata('Ecg.txt');
signal = val(1,1:1800);
Fs = 500;
figure(2);
subplot(3,1,1);
plot(signal);
title('干净的ECG信号');
xlabel('采样点');ylabel('幅值(dB)');
grid on;
XK1=fft(signal,1800);
magXK1=abs(XK1); %幅频特性
figure(3);
subplot(3,1,1);
k1=0:length(magXK1)-1;
stem(k1,magXK1,'.'); %信号幅频特性曲线
xlabel('k');
ylabel('|X(k)|');
title('干净的ECG信号频谱');

% 加入工频信号
av=100;
f0=50;
t=1:length(signal);
noise2=av*cos(2*pi*f0*t/Fs);
signal2=noise2+signal;
figure(2);
subplot(3,1,2);
plot(signal2);
title('含干扰的ECG信号');
xlabel('采样点');ylabel('幅值(dB)');
grid on;
XK2=fft(signal2,1800);
magXK2=abs(XK2); %幅频特性
figure(3);
subplot(3,1,2);
k2=0:length(magXK2)-1;
stem(k2,magXK2,'.'); %信号幅频特性曲线
xlabel('k');
ylabel('|X(k)|');
title('含干扰的ECG信号频谱');

% IIR带阻滤波器设计
wp = [0.18,0.22];
ws = [0.192,0.208];
Rp = 1;
Rs = 15;
[N,Wn] = buttord(wp,ws,Rp,Rs,'s');
[b,a] = butter(N,Wn,'stop');
n=0:0.001:pi;
[H,W] = freqz(b,a,n);
[h,t]=impz(b,a);
y=conv(signal2,h); % 进行带阻滤波
figure(1);
subplot(1,1,1);
plot(W/pi,20*log10(abs(H)));
xlabel('\omega/\pi');
ylabel('幅度 (dB)');
title('IIR带阻频谱');
grid on;
figure(2);
subplot(3,1,3);
plot(y);
title('IIR带阻滤波');
xlabel('采样点');ylabel('幅值(dB)');
grid on;
XK3=fft(y,1800);
magXK3=abs(XK3); %幅频特性
figure(3);
subplot(3,1,3);
k3=0:length(magXK3)-1;
stem(k3,magXK3,'.'); %信号幅频特性曲线
xlabel('k');
ylabel('|X(k)|');
title('IIR带阻滤波后的信号频谱');

四、

clear,clc;
val = importdata('Ecg.txt');
signal = val(1,1:1800);
Fs = 500;
figure(2);
subplot(3,1,1);
plot(signal);
title('干净的ECG信号');
xlabel('采样点');ylabel('幅值(dB)');
grid on;
XK1=fft(signal,1800);
magXK1=abs(XK1); %幅频特性
figure(3);
subplot(3,1,1);
k1=0:length(magXK1)-1;
stem(k1,magXK1,'.'); %信号幅频特性曲线
xlabel('k');
ylabel('|X(k)|');
title('干净的ECG信号频谱');

% 加入基线漂移信号
n1=length(signal)/3;
x1=zeros(1,n1);
t=1:length(signal)-n1;
x2=(length(signal)-n1)/2000*(t-1)+1;
noise3=[x1,x2];
signal3=noise3+signal;
figure(2);
subplot(3,1,2);
plot(signal3);
title('含基线漂移的ECG信号');
xlabel('采样点');ylabel('幅值(dB)');
grid on;
XK2=fft(signal3,1800);
magXK2=abs(XK2); %幅频特性
figure(3);
subplot(3,1,2);
k2=0:length(magXK2)-1;
stem(k2,magXK2,'.'); %信号幅频特性曲线
xlabel('k');
ylabel('|X(k)|');
title('含基线漂移的ECG信号频谱');

% IIR高通滤波器设计
wp = 0.0028*pi;
ws = 0.0012*pi;
Rp = 1;
Rs = 15;
wp1=2*Fs*tan(wp/2);
ws1=2*Fs*tan(ws/2);
[N,Wn]=buttord(wp1,ws1,Rp,Rs,'s');
[b,a]=butter(N,Wn,'high','s');
[Z,P,K]=buttap(N);
[bz,az]=bilinear(b,a,Fs);
[H,W]=freqz(bz,az);
y=filter(bz,az,signal3); % 进行高通滤波
figure(1);
subplot(1,1,1);
plot(W/pi,20*log10(abs(H)));
xlabel('\omega/\pi');
ylabel('幅度 (dB)');
title('IIR高通频谱');
grid on;
figure(2);
subplot(3,1,3);
plot(y);
title('IIR高通滤波');
xlabel('采样点');ylabel('幅值(dB)');
grid on;
XK3=fft(y,1800);
magXK3=abs(XK3); %幅频特性
figure(3);
subplot(3,1,3);
k3=0:length(magXK3)-1;
stem(k3,magXK3,'.'); %信号幅频特性曲线
xlabel('k');
ylabel('|X(k)|');
title('IIR高通滤波后的信号频谱');

三、实验结果与分析：
一、
1.心电信号波形如下图所示

经检验，该结果符合题目要求。
二、
1.IIR低通滤波器和FIR低通滤波器频谱图：

经检验，该结果符合题目要求。
2.干净心电信号波形，含工频干扰的心电信号波形以及IIR低通滤波和FIR低通后的信号波形如下图所示：

经检验，该结果符合题目要求。
3.滤波前后的心电信号的频谱图如下所示：

经检验，该结果符合题目要求。

三、
1.IIR带阻滤波器频谱图：

2.干净心电信号波形，含工频干扰的心电信号波形以及滤波后的信号波形如下图所示：

3.滤波前后的心电信号的频谱图如下所示：

四、
1.IIR高通滤波器频谱图：

2.干净心电信号波形，含基线低频干扰的心电信号波形以及滤波后的信号波形如下图所示：

3.滤波前后的心电信号的频谱图如下所示：

四、实验总结：
通过本次实验，将前面几次数字信号处理课程设计课程中学习到的滤波方法应用到了综合应用中，为了得到不失真的干净的心电信号，需要进行滤波预处理。设计数字低通滤波器、高通滤波器、带阻滤波器，用MATLAB软件对含噪心电信号分别进行高通、带阻和低通滤波等处理，将心电信号中的低频基线漂移、50Hz 工频高频和高频杂波进行滤除。使我了解并掌握了数字滤波器的计算机仿真方法，并通过观察对实际心电图信 号的滤波作用，进一步获得了对数字滤波器的感性认知。.