【⼼电信号】基于matlab⼼电信号特征提取+分析处理【含Matlab源码289期】⼀、获取代码⽅式
获取代码⽅式1：
完整代码已上传我的资源：
获取代码⽅式2：
通过订阅紫极神光博客付费专栏，凭⽀付凭证，私信博主，可获得此代码。
备注：订阅紫极神光博客付费专栏，可免费获得1份代码（有效期为订阅⽇起，三天内有效）；
⼆、⼼电信号简介
0 引⾔
⼼电信号是⼈类最早研究的⽣物信号之⼀, 相⽐其他⽣物信号更易于检测, 且具有直观的规律。⼼电图的准确分析对⼼脏病的及早有重⼤的意义。⼈体是⼀个复杂精密的系统, 有许多不可抗的外界因素, 得到纯净的⼼电信号⾮常困难。可以采⽤神经⽹络算法去除⼼电信号的噪声, 但这种⽅法存在训练难度⼤、耗时长的缺点。⼩波变换在处理⾮线性、⾮平稳且奇异点较多的信号时具有⼀定的优越性, 近年来许多学者使⽤其对⼼电信号进⾏研究。
1 ⼼电信号简介
⼼电信号由以下⼏个波段组成, ⼀个典型的⼼电图如图1所⽰。
图1 典型⼼电图
(1) P波：反映⼼房肌在除极过程中的电位变化过程;
(2) P-R间期：反映的是激动从窦房结通过房室交界区到⼼室肌开始除极的时限;
(3) QRS波：反映⼼室肌除极过程的电位变化;
(4) T波：代表⼼室肌复极过程中所引起的电位变化;
(5) S-T段：从QRS波终点到达T波起点间的⼀段⽔平线[2];
(6) Q-T间期:⼼室从除极到复极的时间[3];
(7) U波：代表动作电位的后电位。
由于⼼电信号⼗分微弱, 且低频, 极易受到⼲扰, 不同的⼲扰源的噪声虽是随机的, 但来⾃同⼀个⼲扰源的噪声往往具有同⼀类特征。分析⼲扰的来源, 针对不同的来源使⽤合适的处理⽅法, 是数据采集重点考虑的⼀个问题。常见⼲扰有3种: (1) ⼯频⼲扰; (2) 基线漂移; (3) 肌电⼲扰。其中已经证明⼩波变换在抑制⼼电信号的⼯频⼲扰⽅⾯具有较⼤优势。具体噪声频带如表1所⽰。
表1 ⼼电信号以及主要噪声频带
三、部分源代码
clear
Fs=360;N=2048;
load('9.mat');
load('9.mat');
x1=y;
load('11.mat');
x2=y;
figure
subplot(2,1,1);plot(x1);grid;title('原始⼼电信号x1.9.mat');
subplot(2,1,2);plot(x2);title('原始⼼电信号x2.9.mat');grid
%(1)谱分析
plot(f1,mag1);axis([0,370,0,190]);grid;%画出频谱图
xlabel('频率(HZ)');ylabel('幅值');title('⼼电信号x1幅度谱');
figure
plot(f2,mag2);axis([0,370,0,190]);grid;%画出频谱图
xlabel('频率(HZ)');ylabel('幅值');title('⼼电信号x2幅度图');
figure
plot(f1,angle(mf1)/pi);axis([0360-1.31.3]);grid;
xlabel('频率(HZ)');ylabel('相位');title('⼼电信号x1相位谱');
figure
plot(f1,angle(mf2)/pi);axis([0360-1.31.3]);grid;
xlabel('频率(HZ)');ylabel('相位');title('⼼电信号x2相位谱');
wn1=x1;
wn2=x2;
P1=10*log10(abs(fft(wn1).^2)/N);matlab傅里叶变换的幅度谱和相位谱
P2=10*log10(abs(fft(wn2).^2)/N);
f1=(0:length(P1)-1)/length(P1);
f2=(0:length(P2)-1)/length(P2);
figure
plot(f1,P1);grid
xlabel('归⼀化频率');ylabel('功率(dB)');title('⼼电信号x1的功率谱');
figure
plot(f2,P2);grid
xlabel('归⼀化频率');ylabel('功率(dB)');title('⼼电信号x2的功率谱');
%（2）相关分析
avr1=mean(x1);avr2=mean(x2);
fprintf('信号x1的均值= %f\n',avr1);fprintf('信号x2的均值= %f\n',avr2);
var1=var(x1);var2=var(x2);
fprintf('信号x1的⽅差= %f\n',var1);fprintf('信号x2的⽅差= %f\n',var2);
rx1=xcorr(x1,'biased');
rx2=xcorr(x2,'biased');
rx1x2=xcorr(x1,x2);
figure
plot(rx1);grid;title('⼼电信号x1的⾃相关函数');
figure
plot(rx2);grid;title('⼼电信号x2的⾃相关函数');
figure
plot(rx1x2);grid;title('⼼电信号x1,x2的互相关函数');
%(3)数字滤波器设计
SNR=10*log(100/8);%2%是能量⽐
x11=awgn(x1,SNR);
figure
subplot(2,1,1),plot(x1);title('原信号x1');%加⼊噪声后有⽑刺，但2%的噪声有点⼩,⽑刺不明显。subplot(2,1,2),plot(x11);title('加⾼斯⽩噪信号');
dt=1/1023.5;%取樣時距(或週期),sec
t=(0:dt:2)';%建⽴⼀個0-2秒的時間向量
y_high=sin(2*pi*1000*t)/10;%100 Hz的⾼頻訊號，振福1/10
y_out=x1+y_high;%基頻載上⼀組⾼頻的輸出。
figure
plot(t,y_out);grid;title('加⼊⾼频噪声1000hz');
%——————IIR零相移数字滤波器纠正基线漂移——————-
Wp=1.4*2/Fs;%通带截⽌频率
Ws=0.6*2/Fs;%阻带截⽌频率
devel=0.005;%通带纹波
Rp=20*log10((1+devel)/(1-devel));%通带纹波系数
Rs=20;%阻带衰减
Rs=20;%阻带衰减
[N, Wn]=ellipord(Wp,Ws,Rp,Rs,'s');%求椭圆滤波器的阶次
[b, a]=ellip(N,Rp,Rs,Wn,'high');%求椭圆滤波器的系数
[hw,w]=freqz(b,a,512);
result =filter(b,a,x1);
figure
freqz(b,a);title('IIR数字滤波器幅频曲线');
figure;subplot(2,1,1);plot(x1);
xlabel('t(s)');ylabel('幅值');title('原始信号');grid
subplot(2,1,2);plot(result);
xlabel('t(s)');ylabel('幅值');title('线性滤波后信号');grid
figure
N=512;
subplot(2,1,1);plot(abs(fft(x1))*2/N);
xlabel('频率(Hz)');ylabel('幅值');title('原始信号x1频谱');grid;
subplot(2,1,2);plot(abs(fft(result))*2/N);
xlabel('频率(Hz)');ylabel('幅值');title('线性滤波去掉基线漂移频谱');grid;
%（5）维纳滤波器去除⼯频⼲扰：
%维纳滤波
Mlag=100;
N=100;%维纳滤波器长度
Rxn=xcorr(x_noise,Mlag,'biased');
Rxnx=xcorr(x1,x_noise,Mlag,'biased');%产⽣输⼊信号与原始信号的互相关函数rxnx=zeros(N,1);
rxnx(:)=Rxnx(101:101+N-1);
Rxx=zeros(N,N);%产⽣输⼊信号⾃相关矩阵
Rxx=diag(Rxn(101)*ones(1,N));
for i=2:N
c=Rxn(101+i)*ones(1,N+1-i);
Rxx=Rxx+diag(c,i-1)+diag(c,-i+1);
end
Rxx;
h=zeros(N,1);
h=inv(Rxx)*rxnx;%计算维纳滤波器的h(n)
yn=filter(h,1,x_noise);%将输⼊信号通过维纳滤波器
figure
plot(yn);title('经过维纳滤波器后信号');
ynmean=mean(yn)%计算经过维纳滤波器后信号均值
ynms=mean(yn.^2)%计算经过维纳滤波器后信号均⽅值
ynvar=var(yn,1)%计算经过维纳滤波器后信号⽅差
Ryn=xcorr(yn,Mlag,'biased');%计算经过维纳滤波器后信号⾃相关函数
Y=fft(yn);%计算经过维纳滤波器后信号序列的快速离散傅⾥叶变换
Y1=fft(x_noise);
Py=Y.*conj(Y)/600;%计算信号频谱
Py1=Y.*conj(Y)/600;
figure
subplot(211)
semilogy(t,Py)%绘制在半对数坐标系下频谱图像
title('经过维纳滤波器后信号在半对数坐标系下频谱图像')
xlabel('频率（HZ）','color','b');ylabel('幅度','color','b')
subplot(212)
semilogy(t,Py1)%绘制在半对数坐标系下频谱图像
title('噪声信号在半对数坐标系下频谱图像')
xlabel('频率（HZ）','color','b');ylabel('幅度','color','b')
pyn=periodogram(yn);%计算经过维纳滤波器后信号的功率谱密度
pyn1=periodogram(x_noise);
figure
subplot(211)
semilogy(pyn)%绘制在半对数坐标系下功率谱密度图像
title('经过维纳滤波器后信号在半对数坐标系下功率谱密度图像');
xlabel('频率','color','b');ylabel('幅度','color','b')
subplot(212)
semilogy(pyn);title('噪声信号在半对数坐标系下功率谱密度图像');
semilogy(pyn);title('噪声信号在半对数坐标系下功率谱密度图像');
xlabel('频率','color','b');ylabel('幅度','color','b')
四、运⾏结果
五、matlab版本及参考⽂献
1 matlab版本
2014a
2 参考⽂献
[1] 沈再阳.精通MATLAB信号处理[M].清华⼤学出版社，2015.
[2]⾼宝建，彭进业，王琳，潘建寿.信号与系统——使⽤MATLAB分析与实现[M].清华⼤学出版社，2020.
[3]王⽂光，魏少明，任欣.信号处理与系统分析的MATLAB实现[M].电⼦⼯业出版社，2018.
[4]焦运良，邢计元，靳尧凯.基于⼩波变换的⼼电信号阈值去噪算法研究[J].信息技术与⽹络安全. 2019,38(05)