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')
评论0