wavelet_check.m

小波压缩编码平台软件包

function [ort,Er,Vr]=wavelet_check(wvf,disp)
%[ort,mmnt]=wavelet_check(wavelet,disp)
%Version: 3.10, Date: 2005/26/12, author: Nikola Sprljan
%Lists properties of a wavelet and performs several tests
%
%Input:
% wvf - wavelet identification string, or wavelet data structure
% disp  - [optional, default = 1] specifies whether to display the output
%            if (disp ~= 0) then display
%            if (disp == 0) then do not display
%
%Output:
% ort - specifies orthogonality of the wavelet:
%        'b'-biorthogonal;
%        'o'-orthogonal;
%        'n'-none
% Er - error energy contribution ratio between the high-pass and low-pass
%      subband
% Vr - error energy ratio between the reconstructed even and odd 
%      samples. Even samples are coincide with the low-pass subsampling 
%      position, while the odd coincide with the high-pass.
%
%Note:
% The function displays (if disp = 1):
%  -filter taps for both analysis and synthesis;
%  -the number of vanishing moments;
%  -the result of the orthogonality test;
%  -DC and Nyquist gain, and L1 and L2 norms
%  -Time-Frequency localisation measure
%  -the result of the Perfect Recostruction (PR) property test for a array of
%   random size and content. Two tests: one for periodic and one for symmetric
%   extension on the matrix borders.
%
%Uses:
% load_wavelet.m (Wavelet Toolbox)
%
%Example:
% wavelet_check('CDF_9x7');

if nargin<2 %display by default
    disp=1;
end;

if ischar(wvf)
    wvf = load_wavelet(wvf);
end;
lLD=length(wvf.filt_AL);
lHD=length(wvf.filt_AH);
lLR=length(wvf.filt_SL);
lHR=length(wvf.filt_SH);
if disp
    fprintf('Analysis (decomposition) filters (LP taps, HP taps) = (%d,%d)\n',lLD,lHD);
    fprintf('LP taps = [');
    for i=1:lLD
        fprintf('%1.5f',wvf.filt_AL(i));
        if (i ~=lLD) fprintf(' '); else fprintf(']\n');end;
    end;
    fprintf('HP taps = [');
    for i=1:lHD
        fprintf('%1.5f',wvf.filt_AH(i));
        if (i ~=lHD) fprintf(' '); else fprintf(']\n');end;
    end;
    fprintf('\nSynthesis (reconstruction) filters (LP taps, HP taps) = (%d,%d)\n',lLR,lHR);
    fprintf('LP taps = [');
    for i=1:lLR
        fprintf('%1.5f',wvf.filt_SL(i));
        if (i ~=lLR) fprintf(' '); else fprintf(']\n');end;
    end;
    fprintf('HP taps = [');
    for i=1:lHR
        fprintf('%1.5f',wvf.filt_SH(i));
        if (i ~=lHR) fprintf(' '); else fprintf(']\n');end;
    end;
end;
%INSTRUCTION: change norm_tolerance depending on how precise the coefficients have to be
norm_tolerance = 10^(-6);
%simple orthogonality check
tmpLoF_D = wvf.filt_AL / sum(abs(wvf.filt_AL));
tmpHiF_D = wvf.filt_AH / sum(abs(wvf.filt_AH));
tmpHiF_R = wvf.filt_SH / sum(abs(wvf.filt_SH));
if (lLD==lHD) & (abs(tmpLoF_D)-fliplr(abs(tmpHiF_D)) < norm_tolerance)
    ort='o';
elseif (lLD==lHR) & (abs(tmpLoF_D)-fliplr(abs(tmpHiF_R)) < norm_tolerance)
    ort='b';
else
    ort='n';
end;
%Er
Er = sum(wvf.filt_SH.^2)/sum(wvf.filt_SL.^2);
%Vr
HiRodd = mod(wvf.filt_SH_delay,2) == 1;
LoRodd = mod(wvf.filt_SL_delay,2) == 1;
%even samples
emse_Lo = sum(wvf.filt_SL(~LoRodd).^2);
emse_Hi = sum(wvf.filt_SH(HiRodd).^2);
%odd samples
omse_Lo = sum(wvf.filt_SL(LoRodd).^2);
omse_Hi = sum(wvf.filt_SH(~HiRodd).^2);
Vr = (omse_Lo + omse_Hi) / (emse_Lo + emse_Hi);

%%%DISPLAY SECTION%%%
if disp
    fprintf('\nOrthogonality (test) = ');
    switch ort
        case 'o'
            fprintf('orthogonal (o)');
        case 'b'
            fprintf('biorthogonal (b)')
        case 'n'
            fprintf('none (n)');
    end;
    fprintf('\nHP/LP reconstructed error energy contribution = %f',Er);
    fprintf('\nOdd samples/even samples reconstructed error energy  = %f',Vr);
    fprintf('\nOrthonormality parameter = %f',orthnormpar(wvf.filt_AL,wvf.filt_AH));
    %vanishing moments for decomposition wavelet
    %INSTRUCTION: change moment_tolerance depending on how precise the computation has to be
    moment_tolerance = 10^(-1);
    plypnts=0:lHD-1;
    %plypnts = plypnts./sum(plypnts);
    mmntD=0;summnt=0;
    while abs(summnt) < moment_tolerance;
        polyf=plypnts.^mmntD;
        summnt=sum(wvf.filt_AH.*polyf);
        mmntD=mmntD+1;
    end;
    %vanishing moments for reconstruction wavelet
    plypnts=0:lHR-1;
    mmntR=0;summnt=0;
    while abs(summnt) < moment_tolerance;
        polyf=plypnts.^mmntR;
        summnt=sum(wvf.filt_SH.*polyf);
        mmntR=mmntR+1;
    end;
    mmnt=[mmntD-1 mmntR-1]; %number of the vanishing moments of the wavelet; 
    %mmnt(1) given by test performed on the analysis HP filter, and mmnt(2) 
    %on the synthesis HP    
    fprintf('\nVanishing moments (test) = (%d,%d)\n',mmnt(1),mmnt(2));
    %HiTmp=wvf.filt_AH;
    %HiTmp(1:2:lHD)=HiTmp(1:2:lHD)*-1;
    fprintf('\nAnalysis LP :\n');
    fprintf('  DC gain = %f \n',abs(sum(wvf.filt_AL)));
    fprintf('  Nyquist gain = %f \n',abs(sum(wvf.filt_AL(1:2:end)) - sum(wvf.filt_AL(2:2:end))));
    fprintf('  L1 norm = %f \n',sum(abs(wvf.filt_AL)));
    fprintf('  L2 norm = %f \n',sqrt(sum(wvf.filt_AL.^2)));
    [TFL,t,w]=tfl(wvf.filt_AL);
    fprintf('  Time-frequency localisation  (t, w) = %f (%f, %f) \n',TFL, t, w);

    fprintf('\nAnalysis HP :\n');
    fprintf('  DC gain = %f \n',abs(sum(wvf.filt_AH)));
    fprintf('  Nyquist gain = %f \n',abs(sum(wvf.filt_AH(1:2:end)) - sum(wvf.filt_AH(2:2:end))));
    fprintf('  L1 norm = %f \n',sum(abs(wvf.filt_AH)));
    fprintf('  L2 norm = %f \n',sqrt(sum(wvf.filt_AH.^2)));
    [TFL,t,w]=tfl(wvf.filt_SL);
    fprintf('  Time-frequency localisation (t, w) = %f (%f, %f) \n',TFL, t, w);

    test_signal_size=32;
    x=rand(1,test_signal_size);
    [a,d] = dwt_lifting1D(x,wvf);
    y = idwt_lifting1D(a,d,wvf);
    fprintf('\nPerfect reconstruction property check for matrix of size %d, Error = %f\n',...
        test_signal_size,max(max(abs(x-y))));
end;


function OP=orthnormpar(LoF_Dx,HiF_Dx)
%Orthonormality parameter, as in:
%M. Lightstone & E. Majani, 揕ow bit-rate design considerations for wavelet
% -based image coding,