newtestmywfun.m

此为基于二代提升小波的阈值压缩算法源码

%*************************************************************************%
% 程序功能：此程序对输入信号进行一维二代整数提升小波变换，对变换后的系数进行阈值处理%
%          后经整数提升小波反变换重构原信号,同时提供了最高层分解的低频系数的十六进制表示，以便于后续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]);