心电信号基于matlab心电信号特征提取+分析处理【含Matlab源码 289期】.zip

clear Fs=360;N=2048; 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)谱分析 mf1=fft(x1,N);%进行频谱变换（傅里叶变换） mf2=fft(x2,N); mag1=abs(mf1); mag2=abs(mf2); f1=(0:length(mf1)-1)*Fs/length(mf1); %进行频率变换 f2=(0:length(mf2)-1)*Fs/length(mf2); figure 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([0 360 -1.3 1.3]);grid; xlabel('频率(HZ)');ylabel('相位');title('心电信号x1相位谱'); figure plot(f1,angle(mf2)/pi);axis([0 360 -1.3 1.3]);grid; xlabel('频率(HZ)');ylabel('相位');title('心电信号x2相位谱'); wn1=x1; wn2=x2; P1=10*log10(abs(fft(wn1).^2)/N); 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; %阻带衰减 [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）维纳滤波器去除工频干扰： dt=1/1023.5; % 取樣時距(或週期),sec t=(0:dt:2)'; % 建立一個0-2秒的時間向量 y=sin(2*pi*50*t)/10; % 100 Hz的高頻訊號,振福1/10 x_noise=x1+y; % 基頻載上一組高頻的輸出。 figure plot(t,x_noise);grid;title('加入噪声50hz'); %维纳滤波 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('噪声信号在半对数坐标系下功率谱密度图像'); xlabel('频率','color','b');ylabel('幅度','color','b')