MIT数据库的心电信号预处理matlab程序

文件包含了matlab程序及从MIT数据库下载的文件，对心电信号进行了低通、带陷、线性滤波处理，去除了干扰信号及进行了基线漂移纠正。

fid=fopen（‘109a.txt‘）;
C=textscan（fid‘%8c %f %*f‘‘headerlines‘2）;
fclose（fid）;
a=C{1};
M=C{2};
 k=length（a）
for i=1:k
   c（i）=strread（a（i:）‘%*s %f‘‘delimiter‘‘:‘）;
end
TIME=c‘;
plot（TIMEM）
clc;
%------------------------------低通滤波器滤除肌电信号------------------------------
Fs=1500;                        %采样频率
fp=80;fs=100;                    %通带截止频率，阻带截止频率
rp=1.4;rs=1.6;                    %通带、阻带衰减
wp=2*pi*fp;ws=2*pi*fs;   
[nwn]=buttord（wpwsrprs‘s‘）;     %‘s‘是确定巴特沃斯模拟滤波器阶次和3dB
[zPk]=buttap（n）;   %设计归一化巴特沃斯模拟低通滤波器，z为极点，p为零点和k为增益
[bpap]=zp2tf（zPk）  %转换为Ha（p）bp为分子系数，ap为分母系数
[bsas]=lp2lp（bpapwp） %Ha（p）转换为低通Ha（s）并去归一化，bs为分子系数，as为分母系数

[hsws]=freqs（bsas）;         %模拟滤波器的幅频响应
[bzaz]=bilinear（bsasFs）;     %对模拟滤波器双线性变换
[h1w1]=freqz（bzaz）;         %数字滤波器的幅频响应
m=filter（bzazM（:1））;

figure
freqz（bzaz）;title（‘巴特沃斯低通滤波器幅频曲线‘）;
      
figure
subplot（211）;
plot（TIMEM（:1））;
xlabel（‘t（s）‘）;ylabel（‘mv‘）;title（‘原始心电信号波形‘）;grid;

subplot（212）;
plot（TIMEm）;
xlabel（‘t（s）‘）;ylabel（‘mv‘）;title（‘低通滤波后的时域图形‘）;grid;
   
N=512
n=0:N-1;
mf=fft（M（:1）N）;               %进行频谱变换（傅里叶变换）
mag=abs（mf）;
f=（0:length（mf）-1）*Fs/length（mf）;  %进行频率变换

figure
subplot（211）
plot（fmag）;axis（[01500150]）;grid;      %画出频谱图
xlabel（‘频率（HZ）‘）;ylabel（‘幅值‘）;title（‘心电信号频谱图‘）;

mfa=fft（mN）;                    %进行频谱变换（傅里叶变换）
maga=abs（mfa）;
fa=（0:length（mfa）-1）*Fs/length（mfa）;  %进行频率变换
subplot（212）
plot（famaga）;axis（[01500150]）;grid;  %画出频谱图
xlabel（‘频率（HZ）‘）;ylabel（‘幅值‘）;title（‘低通滤波后心电信号频谱图‘）;
    
wn=M（:1）;
P=10*log10（abs（fft（wn）.^2）/N）;
f=（0:length（P）-1）/length（P）;
figure
plot（fP）;grid
xlabel（‘归一化频率‘）;ylabel（‘功率（dB）‘）;title（‘心电信号的功率谱‘）;
%-----------------带陷滤波器抑制工频干扰-------------------  
%50Hz陷波器：由一个低通滤波器加上一个高通滤波器组成  
%而高通滤波器由一个全通滤波器减去一个低通滤波器构成  
Me=100;               %滤波器阶数  
L=100;                %窗口长度  
beta=100;             %衰减系数  
Fs=1500;  
wc1=49/Fs*pi;     %wc1为高通滤波器截止频率，对应51Hz  
wc2=51/Fs*pi     ;%wc2为低通滤波器截止频率，对应49Hz  
h=ideal_lp（0.132*piMe）-ideal_lp（wc1Me）+ideal_lp（wc2Me）; %h为陷波器冲击响应  
w=kaiser（Lbeta）;  
y=h.*rot90（w）;         %y为50Hz陷波器冲击响应序列  
m2=filter（y1m）;  
  
figure  
subplot（211）;plot（abs（h））;axis（[0 100 0 0.2]）;  
xlabel（‘频率（Hz）‘）;ylabel（‘幅度（mv）‘）;title（‘陷波器幅度谱‘）;grid;  
N=512;  
P=10*log10（abs（fft（y）.^2）/N）;  
f=（0:length（P）-1）;  
subplot（212）;plot（fP）;  
xlabel（‘频率（Hz）‘）;ylabel（‘功率（dB）‘）;title（‘陷波器功率谱‘）;grid;  
     
figure  
subplot （211）; plot（TIMEm）;  
xlabel（‘t（s）‘）;ylabel（‘幅值‘）;title（‘原始信号‘）;grid;  
subplot（212）;plot（TIMEm2）;  
xlabel（‘t（s）‘）;ylabel（‘幅值‘）;title（‘带阻滤波后信号‘）;grid;  
    
figure  
N=512  
subplot（211）;plot（abs（fft（m））*2/N）;axis（[0 100 0 1]）;  
xlabel（‘t（s）‘）;ylabel（‘幅值‘）;title（‘原始信号频谱‘）;grid;  
subplot（212）;plot（abs（fft（m2））*2/N）;axis（[0 100 0 1]）;  
xlabel（‘t（s）‘）;ylabel（‘幅值‘）;title（‘带阻滤波后信号频谱‘）;grid;    
%------------------IIR零相移数字滤波器纠正基线漂移-------------------  
Wp=1.4*2/Fs;     %通带截止频率   
Ws=0.

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件      118864  2017-12-26 19:55  109a.txt
     文件        4359  2018-01-03 16:19  Untitled3.m