clc ;clear all; close all;
%--------------------------------------------------------------
sfreq=500;
M=dlmread('T103dat.txt');
%M=dlmread('ecg.txt');
[ECG]=Read500HzECGdat(M,-1,sfreq,1024);%500Hz Ecg数据的转化
sfreq=250;
ECG=ECG(405000:408000,2);%截取数据点405000,到408000的数,第二列
%ECG=ECG(50:800,1);
% plot(ECG)
s=ECG;
[C,L] = wavedec(s,11,'db8');%多级一维小波分解,返回小波信号x在N级分解,? N必须是一个严格的正整数(见WMAXLEV)。
%输出分解结构包含小波分解向量C和簿记向量L.对于[C,L] = wavedec(X,N,Lo_D,Hi_D),
%Lo_D是分解低通滤波器
% Hi_D是分解高通滤波器。
%该结构组织为:
%C = [app。COEF(N)|。DET。 系数(N)| ... | det。COEF。(1)]
%L(1)=应用程序的长度。COEF。(N)
% L(i)= det的长度。 (N-i + 2),对于i = 2,...,N + 1
% L(N + 2)=长度(X)。
% a10 = appcoef(C, L, 'db8',10);
% d1 = detcoef(C,L,1);
% d2 = detcoef(C,L,2);
% d3 = detcoef(C,L,3);
% d4 = detcoef(C,L,4);
A9 = wrcoef('a',C,L,'db8',9); % 从近似系数得到近似信号A
D1 = wrcoef('d',C,L,'db8',1); % 从细节系数得到细节信号D
D2 = wrcoef('d',C,L,'db8',2);
D3 = wrcoef('d',C,L,'db8',3);
D4 = wrcoef('d',C,L,'db8',4);
D5 = wrcoef('d',C,L,'db8',5); % 从细节系数得到细节信号D
D6 = wrcoef('d',C,L,'db8',6);
D7 = wrcoef('d',C,L,'db8',7);
D8 = wrcoef('d',C,L,'db8',8);
D9 = wrcoef('d',C,L,'db8',9);
%wrcoef从一维小波系数重建单个分支。
%wrcoef重建一维信号的系数,
%?给出一个小波分解结构(C和L)和
%要么是指定的小波('wname',参见WFILTERS获取更多信息)
% 或指定的重建滤波器(Lo_R和Hi_R)。
% ?
% X = wrcoef('type',C,L,'wname',N)计算的向量
% 基于小波的重构系数
% 分解结构[C,L](更多信息见WAVEDEC),
% ????在N级。“wname”是一个包含小波名称的字符串。
% ??
% ????参数“类型”决定是否近似
% ('type'='a')或细节('type'='d')系数
% 重建。
% 当'type'='a'时,允许N为0;除此以外,
% 严格的正数N是必需的。
% 级别N必须是一个整数,使得N <= length(L)-2。
%
% X = wrcoef('type',C,L,Lo_R,Hi_R,N)计算系数
% 如上所述,鉴于你指定的重建。
% ?
% X = wrcoef('type',C,L,'wname')和
% X = wrcoef('type',C,L,Lo_R,Hi_R)重建系数
% 最大水平N =长度(L)-2。
s1=D1+D2+D3+D4+D5+D6+D7;
s2=smooth(ECG,7,'sgolay');%平滑滤波
%index = 1:1024;
%x = leleccum(index);
%产生噪声信号
% init = 2055615866;
% randn('seed',init);
% nx =s+18*randn(size(s));
%获取消噪的阈值
[thr,sorh,keepapp] = ddencmp('den','wv',s1);
% ddencmp用于去噪或压缩的默认值。
% ????[THR,SORH,KEEPAPP,CRIT] = ddencmp(IN1,IN2,X)
% ????返回去噪或压缩的默认值,
% ????使用小波包或小波包,输入矢量
% ????或矩阵X,其可以是1-D或2-D信号。
% ????THR是门槛,SORH是软或硬
% ????阈值,KEEPAPP允许你保持近似
% ????系数和CRIT(仅用于小波包)
% ????是熵名(见WENTROPY)。
% ????IN1是“den”或“cmp”,IN2是“wv”或“wp”。
% ?
% ????对于小波(三个输出参数):
% ????[THR,SORH,KEEPAPP] = ddencmp(IN1,'wv',X)
% ????返回用于去噪的默认值(如果IN1 =“den”)
% ????或压缩(如果IN1 ='cmp')X.
% ????这些值可以用于WDENCMP。
% ?
% ????对于小波包(四个输出参数):
% ????[THR,SORH,KEEPAPP,CRIT] = ddencmp(IN1,'wp',X)
% ????返回用于去噪的默认值(如果IN1 =“den”)
% ????或压缩(如果IN1 ='cmp')X.
% ????这些值可以用于WPDENCMP。
%对信号进行消噪
xd = wdencmp('gbl',s1,'db4',2,thr,sorh,keepapp);
figure(3)
subplot(221);
plot(s);
title('原始信号');
subplot(222);
plot(s1);
title('含噪信号');
subplot(223);
plot(xd);
title('消噪后的信号');
subplot(224);
plot(s2);
title('平滑滤波信号');
% wdencmp使用小波去噪或压缩。
% ????wdencmp执行去噪或压缩过程
% ????信号或图像使用小波。
% ?
% ????[XC,CXC,LXC,PERF0,PERFL2] =
% ????wdencmp( 'GBL',X, 'wname',N,THR,SORH,KEEPAPP)
% ????返回输入的去噪声或压缩版本XC
% ????信号X(1-D或2-D)由小波系数获得
% ????使用全局正阈值THR进行阈值化。
% ????额外的输出参数[CXC,LXC]是
% ????XC的小波分解结构,
% ????PERFL2和PERF0是L ^ 2的恢复和压缩
% ????百分比分数。
% ????PERFL2 = 100 *(C的向量范数/ C的向量范数)^ 2
% ????其中[C,L]表示小波分解结构
% ????X.
% ????小波分解在N和N级执行
% ????'wname'是一个包含小波名称的字符串。
% ????SORH('s'或'h')用于软或硬阈值
% ????(请参阅WTHRESH了解更多详情)。
% ????如果KEEPAPP = 1,近似系数不能
% ????thresholded,否则是可能的。
% ?
% ????wdencmp( 'GBL',C,L,W,N,THR,SORH,KEEPAPP)
% ????具有相同的输出参数,使用相同的
% ????选项如上,但直接从中获得
% ????输入小波分解结构[C,L]的
% ????信号在N级被去噪或压缩,
% ????使用'wname'小波。
% ?
% ????对于一维情况和“lvd”选项:
% ????wdencmp('lvd',X,'wname',N,THR,SORH)或
% ????wdencmp('lvd',C,L,'wname',N,THR,SORH)
% ????有相同的输出参数,使用相同的
% ????选项如上,但允许水平相关
% ????向量THR中包含的阈值(THR必须是
% ????长度N)。另外,近似值保持不变。
% ?
% ????对于二维案例和“lvd”选项:
% ????wdencmp('lvd',X,'wname',N,THR,SORH)或
% ????wdencmp('lvd',C,L,'wname',N,THR,SORH)
% ????THR必须是包含该级别的N个矩阵
% ????依赖于三个方向的阈值
% ????水平,对角线和垂直。
% function[stp2_tim,stp2_value]=find_up_peaks(s,sfreq)
% %此函数的功能是找到zn中的所有的向上的峰值,
%
% %---------------------------------------------
% [sample,C]=size(s);
%
% %----------------求差分--------------------------
% znDEF=zeros(sample,C);
% for load=1:C
% for i=2:sample-1
% znDEF(i,load)=zn(i+1,load)-zn(i,load);
% end
% end
% %--------------------求向上的峰值对应的过零点--------------------
% stph_tim=zeros(1,C);
% stph_value=zeros(1,C);
% inter=ceil(0.5*sfreq);
% for load=1:C
% k=0;
% for i=1:sample-3 %i=inter:sample-inter
% if ( znDEF(i,load)>0 & znDEF(i+1,load)<0 ) | (znDEF(i,load)*znDEF(i+1,load)==0 & znDEF(i,load)>0 &znDEF(i+2,load)<0) | (znDEF(i,load)*znDEF(i+1,load)==0 & znDEF(i,load)>0 & znDEF(i+3,load)<0 & znDEF(i+2)>0)
% k=k+1;
% stph_tim(k,load)=i+1;
% stph_value(k,load)=zn(i+1,load);
% end
% end
% end
% stp2_tim=stph_tim;
% stp2_value=stph_value;
% end
% figure(3)
% plot(stp2_tim,stp2_value);
figure(1)
subplot(311)
plot(s)%除去基线漂移,对于心电信号一般是1Hz
title('原始信号')
xlim([0,3000])
ylim([-200,200])
xlabel('采样点数')
ylabel('幅度')
subplot(312)
plot(s1)
title('去除基线漂移后的信号')
xlim([0,3000])
% xlim X限制。
% ???? XL = xlim获取当前轴的x极限。
% ???? xlim([XMIN XMAX])设置x限制。
% ???? XLMODE = xlim('mode')获取x限制模式。
% ???? xlim(模式)设置x限制模式。
% ????????????????????????????? (模式可以是“自动”或“手动”)
% ???? xlim(AX,...)使用轴AX而不是当前轴。
% ?
% ???? xlim设置或获取轴的XLim或XLimMode属性。
% xDFT = fft(s);
% Freq = 0:50;
% subplot(211);
% plot(s); xlabel('Seconds'); ylabel('Amplitude');
% subplot(212);
% plot(Freq,abs(xDFT(1:length(xDFT)/2+1)))
% set(gca,'xtick',[4:4:64]);
% xlabel('Hz'); ylabel('Magnitude');
xlabel('采样点数')
ylabel('幅度')
subplot(313)
plot(s-s1)
title('基线')
xlabel('采样点数')
ylabel('幅度')
xlim([0,3000])
figure(2)
subplot(911)
plot(s)
xlim([0,3000])
xlabel('采样点数')
ylabel('幅度')
title('原始信号')
subplot(912)
plot(D1)
xlim([0,3000])
xlabel('采样点数')
ylabel('幅度')
title('D1')
subplot(913)
plot(D2)
xlim([0,3000])
xlabel('采样点数')
ylabel('幅度')
title('D2')
subplot(914)
plot(D3)
xlim([0,3000])
xlabel('采样点数')
ylabel('幅度')
title('D3')
subplot(915)
plot(D4)
xlim([0,3000])
xlabel('采样点数')
ylabel('幅度')
title('D4')
subplot(916)
plot(D5)
xlim([0,3000])
xlabel('采样点数')
ylabel('幅度')
title('D5')
subplot(917)
plot(D6)
xlim([0,3000])
xlabel('采样点数')
ylabel('幅度')
title('D6')
subplot(918)
plot(D7)
xlim([0,3000])
xlabel('采样点数')
ylabel('幅度')
title('D7')
subplot(919)
plot(D8)
xlim([0,3000])
xlabel('采样点数')
ylabel('幅度')
title('D8')
- 1
- 2
前往页