clear all;
close all;
clc;
sy=load('F:\伸腕1.txt');
indx=100:1100;
s=sy(indx);
subplot(3,1,1);
plot(s);grid on;
axis([1 1000 -1.5 1.5]);
title('原始信號');
%t=0:1000;
%st=sin(0.03*t);
%subplot(4,1,1);
%plot(st);
%axis([1 1000 -1.5 1.5]);
%title('原始信號');
%load noissin;
%s=noissin;
%subplot(4,1,2);
%plot(s);
%title('染噪信號');




%使用sym8小波，得到相應的提升方案
issym8=liftwave('sym8');
%添加els到提升方案
els={'p',[-0.125 0.125],0}
isnew=addlift(issym8,els);

xDec = lwt(s,isnew,3);
% 提取低頻和高頻系數
ca1 = lwtcoef('ca',xDec,isnew,3,1);
cd1 = lwtcoef('cd',xDec,isnew,3,1);
ca2 = lwtcoef('ca',xDec,isnew,3,2);
cd2 = lwtcoef('cd',xDec,isnew,3,2);
ca3 = lwtcoef('ca',xDec,isnew,3,3);
cd3 = lwtcoef('cd',xDec,isnew,3,3);

% the first level ...
w = sort(abs(cd1));
n=length(cd1);
if rem(n, 2) == 1
    m= w((n+1) / 2);
else
    m= (w(n/2) + w(n/2+1)) / 2;
end
%閾值選取
e1 = abs(m) ./ 0.6745;
Thr1 = e1 * (2*log(n)) .^(1/2);
%閾值對高頻的處理
for i=1:length(cd1)
    if(abs(cd1(i))>=Thr1)
        cd1(i)=cd1(i);
    else 
        cd1(i)=0;
    end
end;   
    

% the second level ...
w = sort(abs(cd2));
n=length(cd2);
if rem(n, 2) == 1
    m= w((n+1) / 2);
else
    m= (w(n/2) + w(n/2+1)) / 2;
end
%閾值選取
e1 = abs(m) ./ 0.6745;
Thr1 = e1 * (2*log(n)) .^(1/2);
%閾值對高頻的處理
for i=1:length(cd2)
    if(abs(cd2(i))>=Thr1)
        cd2(i)=cd2(i);
    else
        cd2(i)=0;
    end
end; 
    
w = sort(abs(cd3));
n=length(cd3);
if rem(n, 2) == 1
    m= w((n+1) / 2);
else
    m= (w(n/2) + w(n/2+1)) / 2;
end
%閾值選取
e1 = abs(m) ./ 0.6745;
Thr1 = e1 * (2*log(n)) .^(1/2);
%閾值對高頻的處理
for i=1:length(cd3)
    if(abs(cd3(i))>=Thr1)
        cd3(i)=cd3(i);
    else
        cd3(i)=0;
    end
end; 

    

s2 = ilwt(ca3,cd3,isnew);
s1 = ilwt(s2,cd2,isnew);
s0 = ilwt(s1,cd1,isnew);

% 均方根誤差：RMSE
RMSE1 = (sum((s - s0).^2) / length(s)).^(1/2);
RMSE2 = (sum((s - S).^2) / length(s)).^(1/2);
% 信噪比:SNR
SNR1 = 10 * log(sum(s.^2)/sum((s - s0).^2));
SNR2 = 10 * log(sum(s.^2)/sum((s - S).^2));

% 峰值誤差：PE
PE1 = max(abs(s - s0));
PE2 = max(abs(s - S));
subplot(3,1,1);
plot(s);
axis([1 1000 -1.5 1.5]);
title('the original signal');grid on;
xlabel('t/ms');
ylabel('s/mv');
subplot(3,1,2);
plot(S);grid on;
axis([1 1000 -1.5 1.5]);
title('the signal after being de-noised by the first generation wavelet');
xlabel('t/ms');
ylabel('s/mv');
subplot(3,1,3);
plot(s0); grid on;
axis([1 1000 -1.5 1.5]);
title('the signal after being de-noised by the lifting wavelet');
xlabel('t/ms');
ylabel('s/mv');
%
% %使用默認閾值消噪
% [c,l]=wavedec(s,1,'sym8');
% [thr,sorh,keepapp]=ddencmp('den','wv',s);
% s3=wdencmp('gbl',c,l,'sym8',1,thr,sorh,keepapp);
% figure(1);
% subplot(2,1,1),plot(s),title('原始信號');
% ylabel('幅值');
% %subplot(4,1,2),plot(s2),title('強制消噪后的信號');
% xlabel('樣本序號n');
% ylabel('幅值');
% %subplot(4,1,3),plot(s3),title('默認閾值消噪后的信號');
% xlabel('樣本序號n');
% ylabel('幅值');

% Subfunction
% Calculate median.

% %去噪后的提升反變換
% s4=ilwt(ca1,cd1,isnew);
% 
% subplot(2,1,2),plot(s4),title('軟閾值消噪后的信號');
%xlabel('t/ms');
%ylabel('s/mv');