%*************************************************************************%
% 程序功能：此程序對輸入信號進行一維二代整數提升小波變換，對變換后的系數進行閾值處理%
%          后經整數提升小波反變換重構原信號,同時提供了最高層分解的低頻系數的十六進制表示，以便于后續RAY-PERIOD壓縮算法的進行
% 編寫人：孫金鳳                                                           %
% 編寫時間：2006年3月10日--4月10日                                          %
%*************************************************************************%

% 初始化
clc;
clear;
close all;
% 信號的裝載
i=220*faultnoise();
%load leleccum;
%i=leleccum;
i=round(i);
lx=length(i);
%------------------------------------------------------
% 程序參數設置
wname='db4';%使用的小波
N=6;% 設置分解層次
%------------------------------------------------------
%計算所用小波的提升方案
LS=LiftScheme(wname,'int2int');

% 顯示濾波器相關信息
disp('使用的小波名稱為：');disp(wname);
disp('分解層數為：');disp(N);
disp('其小波濾波器組的勞倫多項式如下：');
[Hs,Gs,Ha,Ga]=wave2lp(wname)
disp('對應的提升方案為：');
displs(LS); 

tic  %開始計時
%計算提升小波分解得到的全部尺度和小波系數
xDec=MyLwt(i,LS,N); %注意：原位運算造成尺度與小波系數交織存儲
% 提取最高層的低頻系數
c=mylwtcoef('ca',xDec,LS,N,N);
la=length(c);
%-------------------------------------------------
% 設置計算閾值的參數
m=la;
alpha=1.5;
%-------------------------------------------------
dc=[];ds=[];
x=c;
thr=[];
% 提取各層高頻系數，計算局部閾值，并進行閾值處理
for j=N:-1:1
    d=mylwtcoef('cd',xDec,LS,N,j);
    % 計算閾值
    n = m/(N+2-j)^alpha; 
    n = min(round(n),length(d));
    if n == length(d)
        thr(j) = 0;
    else
        dd = sort(abs(d));
        thr(j) = dd(end-n);
    end
    nkeep(j) = n;%每層保留系數的個數
    % 進行閾值處理
    dt = thresh(d,'h',thr(j));
    x = myilwt(x,dt,LS);
    dc = [dt dc];%存儲處理之后的高頻系數
    ds = [d ds];%存儲處理之前的高頻系數
end
%將低頻系數轉化為16進制，以便于RAY-PERIOD壓縮算法的進行。
[abin,ind2]=DectoBin(c,32);%將低頻系數轉化為2進制
[row,col]=size(abin);
acomb=[];%初始化組合在一起的低頻系數
for j=1:row
    ahex(j,:)=BintoHex(abin(j,:),32);
    acomb=[acomb ahex(j,:)];
end
disp('最高層分解的低頻系數為：')
ahex;
acomb
ind2 %存儲負系數的位置
disp('低頻系數的長度為:')
la
%-----------------------------------------------------------
% 性能的估計
%-----------------------------------------------------------
disp('各層高頻的閾值為：');
thr
disp('各層高頻的保留系數個數為：');
nkeep
% 計算閾值處理后零值占所有數據量的百分比
disp('系數為0的個數百分比：');
perf0 = 100*(length(find(dc==0))/length(dc))
perf00=100*(length(find(dc==0))/length(xDec));
% 壓縮比
if perf0<=96.875;
    z=ones(1,length(dc));
    ind16=find(dc==0);
    z([ind16])=0;
    lz=length(z)
    disp('非零代表位圖為');
    hex=BintoHex(z,lz)
    ddcr=(100-perf00)+(3.125*lz+100*la)/length(xDec);
else
    ddcr=(100-perf00)*2+100*la/length(xDec);
end
disp('閾值壓縮比為')
ddcr
% 計算能量恢復的百分比
nc = norm(xDec);
if nc<eps
    perfl2 = 100;
else
    cxDec=[c,dc];
    perfl2 = 100*((norm(cxDec)/nc)^2);
    disp('保留能量百分比：');
    sprintf('%.10f',perfl2)
end
disp('相對誤差')
100*norm(dc-ds)/norm(xDec)
toc
%-----------------------------------------------------------
% 壓縮結果的圖像表示
%-----------------------------------------------------------
[Lo_D,Hi_D,Lo_R,Hi_R]=ls2filt(LS);
bswfun(Lo_D,Hi_D,Lo_R,Hi_R,'plot');%分析與綜合尺度及小波函數的圖像顯示
figure(2);
subplot(411);plot(i);title('原始信號');
dd=[c ds];
subplot(412);plot(dd);title('分解后的能量分布');
subplot(413);plot(x);title('重構信號');
y = i-x;
subplot(414);plot(y);title('差值信號');
xlabel1 = ['2-norm rec:',num2str(sprintf('%.10f',perfl2))];
xlabel2 = ['%     ---zero cfs:',num2str(perf0),'%'];
xlabel([xlabel1 xlabel2]);