%% ================================讀取ECG訊號=============================%%
clc;clear all;
addpath(genpath('D:\PLMS\ECGDataProcess\Mit-Bih-ECG-Signal'));
 [TIME,M,Fs,siginfo]=rdmat('100m');

%% ===============================小波閾值去噪============================= %%
%應用db5作為小波函式進行三層分解
%利用無偏似然估計閾值
%對100.dat的單導聯資料進行去噪
E1
 
=M(:,1);   % M為第二篇中對100.dat檔案處理後得到的資料矩陣，M(:,1)指MLII導聯資料，M(:,2)指V5導聯資料
E1=E1';
n1=size(E1);
figure;
plot(E1)
% s1=E1(1:2000);
%小波分解
[C1 L1]=wavedec(E1,6,'db5');
%從C中提取尺度3下的近似小波係數
cA6_1=appcoef(C1,L1,'db5',6);

%從訊號C中提取尺度1,2,3下的細節小波係數
cD1_1=detcoef(C1,L1,1);
cD2_1=detcoef(C1,L1,2);
cD3_1=detcoef(C1,L1,3);
cD4_1
 
=detcoef(C1,L1,4);
cD5_1=detcoef(C1,L1,5);
cD6_1=detcoef(C1,L1,6);
figure;
subplot(711);plot(cA6_1);ylabel('近似訊號a6');   %畫出各層小波係數
title('小波分解示意圖');
subplot(712);plot(cD6_1);ylabel('細節訊號d6');
subplot(713);plot(cD5_1);ylabel('細節訊號d5');
subplot(714);plot(cD4_1);ylabel('細節訊號d4');
subplot(715);plot(cD3_1);ylabel('
 
細節訊號d3');
subplot(716);plot(cD2_1);ylabel('細節訊號d2');
subplot(717);plot(cD1_1);ylabel('細節訊號d1');
xlabel('取樣點數'); 

%使用stein的無偏似然估計原理進行選擇各層的閾值
%cD1_1,cD2_1,cD3_1為各層小波係數
%rigrsure為無偏似然估計的閾值型別
thr=thselect(E1,'rigrsure'); 
ysoft6=wthresh(cD6_1,'s',thr);
ysoft5=wthresh(cD5_1,'s',thr);
ysoft4=wthresh(cD4_1,'s',thr);
ysoft3=wthresh(cD3_1,'s',thr);
ysoft2=wthresh(cD2_1,'s',thr);
ysoft1 = wthresh(cD1_1,'s',thr);
cA6_1 = zeros(1,346);
C = [cA6_1,ysoft6,ysoft5,ysoft4,ysoft3,ysoft2,ysoft1];
XC1 = waverec(C,L1,'db5');
% thr2_1=thselect(cD2_1,'rigrsure');
% thr3_1=thselect(cD3_1,'rigrsure');
% thr4_1=thselect(cD4_1,'rigrsure');
% thr5_1=thselect(cD5_1,'rigrsure');
% thr6_1=thselect(cD6_1,'rigrsure');
%各層的閾值
TR1=[2,thr6_1,thr5_1,thr4_1,thr3_1,thr2_1];
%'s'為軟閾值，'h'為硬閾值
SORH1='s';
% [thr,sorh,keepapp] = ddencmp('den','wv',E1);
% [XC1,CXC1,LXC1,PERF0_1,PERF2_1] = wdencmp('gbl',E1,'db5',6,thr,sorh,keepapp);
%----------去噪----------
% XC為去噪後訊號
% [CXC,LXC]為小波分解結構
% PERF0和PERF2是恢復和壓縮的範數百分比
% 'lvd'為允許設定各層的閾值
% 'gbl'為固定閾值
% 3為閾值的長度

%  [XC1,CXC1,LXC1,PERF0_1,PERF2_1]=wdencmp('lvd',E1,'db5',6,TR1,SORH1);

%-----------利用waverec重構----------------
cD1_1 = zeros(1,10804);

cA6_1 = zeros(1,346);
C = [cA6_1,cD6_1,cD5_1,cD4_1,cD3_1,cD2_1,cD1_1];
XC1 = waverec(C,L1,'db5');

%畫圖
figure;
subplot(2,1,1);plot(TIME,E1);title('ECG Signal 100.mat 原訊號 V5（10秒）');
xlabel('Time/s');ylabel('Voltage/mV');
xlim([0 10]);
subplot(2,1,2);plot(TIME,XC1);title('去噪後訊號 V5（10秒）');
xlabel('Time/s');ylabel('Voltage/mV');
xlim([0 10]);
%----------去噪效果衡量----------
%SNR越大效果越好，MSE越小越好
%選取訊號的長度
N1=n1(2);
x1=E1;
y1=XC1;
F1=0;
MM1=0;
for ii=1:N1
    m1(ii)=(x1(ii)-y1(ii))^2;
    t1(ii)=y1(ii)^2;
    f1(ii)=t1(ii)/m1(ii);
    F1=F1+f1(ii);
    MM1=MM1+m1(ii);
end;
SNR1=10*log10(F1);
MSE1=MM1/N1;
SM1=SNR1/MSE1;

%列印結果
disp('****************訊號2****************')
disp('**                                    **')

disp(['**','    SNR1=',num2str(SNR1),'    MSE1=',num2str(MSE1),'  **']);
disp('**                                    **')
disp('****************去噪效果****************')