心电信号分析与处理

心电信号原始数据

资料与题目

心脏内部产生的一系列非常协调的电刺激脉冲，使得心脏肌肉细胞有节奏的舒张和收缩，这些信号传递到人体表面的不同部位形成不同的电位差。通过仪器设备可以从体表检测到这些微弱的电位差信号，称之为心电信号。

正常的心电信号频率范围为 0.05Hz－100Hz，其能量集中在低频段，其中 99％的能量集中在0Hz－35Hz。在其采集过程中容易受到各种干扰，主要分为以下：①工频和工频的谐波频率干扰，工频频率在我国为50Hz；②肌颤噪声和采样电路参考电压引入的电源纹波等高频噪声，频率通常在100Hz以上；③呼吸基线漂移和采样引入的直流分量，频率一般分布在0-0.7Hz。

上图所示是某电子设备采样得到的4096点心电信号原始数据时域波形，该电子设备的采样频率为500Hz，每采样点数据使用16bits存储，分为2个8bits，其中低位在前，高位在后。文件名：serial_data_raw.dat，文件中至少包含4096点原始心电信号数据，读取参考代码见附件1。

现需要设计处理算法，消除该原始数据中存在的干扰，具体要求如下：

1 分析该原始数据的频谱；

2 设计一个数字直流陷波器，消除呼吸基线漂移和采样引入的直流分量，0频处衰减不低于30dB；

直流陷波器参考设计为： ，其中a为接近1的正实数；

3 用加窗法设计一个FIR数字低通滤波器，消除工频和各种高频干扰，保留频率35Hz以内的心电信号，阻带衰减不低于40dB（不能使用FDATOOL或者filterdesigner），过渡带不超过10Hz；

4 经过2和3滤波后的时域心电信号波形，需要绘制在320*240的屏幕上，因此需要把滤波后的时域波形的数据动态范围限制在120左右（即时域数据最大值和最小值相差在120左右）；

5 分析并比较每次滤波前后数据的频谱，横坐标需要处理为Hz；

6 估算心电信号的心率的频率。

附件一

读取数据的程序代码

dataLen = 4096;

stm32data=zeros(1,dataLen);

tmpOut = zeros(1,dataLen*2);

fid=fopen('serial_data_raw.dat','r');

tmpOut = fread(fid,dataLen*2,'uchar');

for i=1:dataLen

    stm32data(i)=(tmpOut((i-1)*2+1)*1+tmpOut((i-1)*2+2)*256);   

end

数据样式

附件二

工程中使用高精度传感器采集动态信号，采集的原始数据会因为环境变化、量化字长和参考电压等因素，包含较强的缓变直流分量，如果不予以消除，会导致在降噪等处理中出现运算饱和溢出。例如在采集心电信号时，因为参考电压受环境温度变化会产生一定的温漂，以及人的呼吸活动和电极滑动也导致基线漂移。这些干扰的频率很低，通常在几Hz以内，但和心电信号的有效频谱非常接近，因此需要过渡带较窄的IIR直流陷波器来消除干扰，其系统函数H(z)为。

数字系统H(z)的极点为z = a，零点为z =1，其中a为接近1的正实数。因为零点对应的幅角，

所以如图4-5-1所示，数字系统在零频处的增益显著衰减，其衰减的程度与过渡带的宽度和极点a的数值有关。采用系统频率响应的几何确定和计算机仿真，都可以得出结论：a越大衰减变小，同时过渡带变窄，反之则衰减增大和过渡带增宽。因此只需要根据实际需要，通过计算机仿真，就可以求得合适的a值。

代码

第零步——读取数据

clear;clc;
close all;
%读取心电信号
dataLen = 4096;
data=zeros(1,dataLen);
tmpOut = zeros(1,dataLen*2);
fid=fopen('serial_data_raw.dat','r');
tmpOut = fread(fid,dataLen*2,'uchar');
for i=1:dataLen
    data(i)=(tmpOut((i-1)*2+1)*1+tmpOut((i-1)*2+2)*256);    
end

第一步——原始心电信号时域与频域

%1
fs=500;%心电信号采样频率
t = 0:1/fs:(dataLen-1)/fs;%原始信号时间轴
%绘制图象
figure(1);
subplot(2,1,1);
plot(t,data);
title('原始心电信号时域图');
xlabel('时间/s');
ylabel('幅度');

X1=fftshift(fft(data));
fx = (-dataLen/2:dataLen/2-1)*(fs/dataLen);%横坐标/Hz
db = 20 * log10(abs(X1));%计算分贝
figure(1);
subplot(2,1,2);
plot(fx,db);
ylim([0 120]);
title('原始心电信号频谱图');
xlabel('频率/Hz');
ylabel('幅度');

第二步——直流陷波器

a=0.991;
B=[1,-1];%分子系数矢量
A=[1,-a];%分母系数矢量
w=-0.5*pi:0.1:0.5*pi;
H=freqz(B,A,w);
%绘制图象
y=filter(B,A,data);
figure(2);
subplot(2,1,1);
plot(t,y);
title('陷波后心电信号时域图');
xlabel('时间/s');
ylabel('幅度');
axis([1.2 8 -60 60]);

Y=fftshift(fft(y));
fy = (-dataLen/2:dataLen/2-1)*(fs/dataLen);%横坐标/Hz
db2 = 20 * log10(abs(Y));%计算分贝
subplot(2,1,2);
plot(fy,db2);
ylim([0 120]);
title('陷波后心电信号频谱图');
xlabel('频率/Hz');
ylabel('幅度');

直流限波器滤除0Hz附近的低频噪声，即呼吸基线漂移和采样引入的直流分量

第三步——加窗法

%窗函数低通滤波器
fc = 35; % 截止频率
Astop = 40; % 阻带衰减9
ftrans = 10; % 过渡带宽度
N=6.6/(ftrans/fs*2);
% 使用窗函数设计FIR低通滤波器
b = fir1(N, fc/fs, 'low', hann(N+1));
y2=filter(b,1,y);
Y2=fftshift(fft(y2));
db3=20 * log10(abs(Y2));%计算分贝
fy2 = (-dataLen/2:dataLen/2-1)*(fs/dataLen);%横坐标/Hz
figure(3);
subplot(2,1,1);
plot(t,y2);
title('加窗法后心电信号时域图');
xlabel('时间/s');
ylabel('幅度');

axis([1.2 8 -60 60]);
subplot(2,1,2);
plot(fy2,db3);
ylim([0 120]);
title('加窗法后心电信号频谱');
xlabel('频率/Hz');
ylabel('幅度');

采用汉宁窗设计滤波器，设置过渡带宽为6.6*pi/N，要求过渡带宽为10Hz，由此可以计算求出滤波器阶数N=165，高频噪声基本滤除。加窗法滤除50Hz及以上的高频噪声，即50Hz的工频噪声以及100Hz肌颤噪声和采样电路参考电压引入的电源纹波等高频噪声。

第四步——放缩

figure(4);
plot(t,y2);
title('放缩后心电信号时域图');
xlabel('时间/s');
ylabel('幅度');
axis([2 3 -60 60]);
y3=y2(1000:2000);
figure;
plot(y3);

第五步——检索心率

k=0;
j=0;
while(j<1000)
    j=j+1;
    if y3(j)>20
        if k==0
            t1=j;
            j=j+200;
            k=1;
        end
        if k==1
            t2=j;
            j=j+200;
        end
    end
end
freq=1/((t2-t1)/fs);
fprintf("两点间隔为");
disp(t2-t1);
fprintf("心电信号频率为");
disp(freq);

峰值间隔为500，采样频率为500Hz，所以间隔为1s，折合心率为60次/min即1Hz