perform_wavelet_transform.m

是一本非常不错的书

%
% Part of WaveLab Version 802
% Built Sunday, October 3, 1999 8:52:27 AM
% This is Copyrighted Material
% For Copying permissions see COPYING.m
% Comments? e-mail [e-mail]
%

function y = UpDyadHi(x,qmf)
% UpDyadHi -- Hi-Pass Upsampling operator; periodized
%  Usage
%    u = UpDyadHi(d,f)
%  Inputs
%    d    1-d signal at coarser scale
%    f    filter
%  Outputs
%    u    1-d signal at finer scale
%
%  See Also
%    DownDyadLo, DownDyadHi, UpDyadLo, IWT_PO, aconv
%

y = aconv( MirrorFilt(qmf), rshift( UpSample(x) ) );

%
% Copyright (c) 1993. Iain M. Johnstone
%


%
% Part of WaveLab Version 802
% Built Sunday, October 3, 1999 8:52:27 AM
% This is Copyrighted Material
% For Copying permissions see COPYING.m
% Comments? e-mail [e-mail]
%

function y = rshift(x)
% rshift -- Circular right shift of 1-d signal
%  Usage
%    r = rshift(x)
%  Inputs
%    x   1-d signal
%  Outputs
%    r   1-d signal
%        r(i) = x(i-1) except r(1) = x(n)
%

n = length(x);
y = [ x(n) x( 1: (n-1) )];

%
% Copyright (c) 1993. Iain M. Johnstone
%


%
% Part of WaveLab Version 802
% Built Sunday, October 3, 1999 8:52:27 AM
% This is Copyrighted Material
% For Copying permissions see COPYING.m
% Comments? e-mail [e-mail]
%




function [qmf,dqmf] = MakeBSFilter(Type,Par)
% MakeBSFilter -- Generate Biorthonormal QMF Filter Pairs
%  Usage
%    [qmf,dqmf] = MakeBSFilter(Type,Par)
%  Inputs
%    Type   string, one of:
%             'Triangle'
%             'Interpolating' 'Deslauriers' (the two are same)
%             'Average-Interpolating'
%             'CDF' (spline biorthogonal filters in Daubechies's book)
%             'Villasenor' (Villasenor's 5 best biorthogonal filters)
%    Par    integer list, e.g. if Type ='Deslauriers', Par=3 specifies
%           Deslauriers-Dubuc filter, polynomial degree 3
%  Outputs
%    qmf    quadrature mirror filter  (odd length, symmetric)
%    dqmf   dual quadrature mirror filter  (odd length, symmetric)
%
%  See Also
%    FWT_PBS, IWT_PBS, FWT2_PBS, IWT2_PBS
%
%  References
%    I. Daubechies, "Ten Lectures on Wavelets."
%
%    G. Deslauriers and S. Dubuc, "Symmetric Iterative Interpolating Processes."
%
%    D. Donoho, "Smooth Wavelet Decompositions with Blocky Coefficient Kernels."
%
%    J. Villasenor, B. Belzer and J. Liao, "Wavelet Filter Evaluation for
%    Image Compression."
%

if nargin < 2,
    Par = 0;
end

sqr2 = sqrt(2);

if strcmp(Type,'Triangle'),
    qmf = [0 1 0];
    dqmf = [.5 1 .5];

elseif strcmp(Type,'Interpolating') | strcmp(Type,'Deslauriers'),
    qmf  = [0 1 0];
    dqmf = MakeDDFilter(Par)';
    dqmf =  dqmf(1:(length(dqmf)-1));

elseif strcmp(Type,'Average-Interpolating'),
    qmf  = [0 .5 .5] ;
    dqmf = [0 ; MakeAIFilter(Par)]';

elseif strcmp(Type,'CDF'),
    if Par(1)==1,
        dqmf = [0 .5 .5] .* sqr2;
        if Par(2) == 1,
            qmf = [0 .5 .5] .* sqr2;
        elseif Par(2) == 3,
            qmf = [0 -1 1 8 8 1 -1] .* sqr2 / 16;
        elseif Par(2) == 5,
            qmf = [0 3 -3 -22 22 128 128 22 -22 -3 3].*sqr2/256;
        end
    elseif Par(1)==2,
        dqmf = [.25 .5 .25] .* sqr2;
        if Par(2)==2,
            qmf = [-.125 .25 .75 .25 -.125] .* sqr2;
        elseif Par(2)==4,
            qmf = [3 -6 -16 38 90 38 -16 -6 3] .* (sqr2/128);
        elseif Par(2)==6,
            qmf = [-5 10 34 -78 -123 324 700 324 -123 -78 34 10 -5 ] .* (sqr2/1024);
        elseif Par(2)==8,
            qmf = [35 -70 -300 670 1228 -3126 -3796 10718 22050 ...
                10718 -3796 -3126 1228 670 -300 -70 35 ] .* (sqr2/32768);
        end
    elseif Par(1)==3,
        dqmf = [0 .125 .375 .375 .125] .* sqr2;
        if Par(2) == 1,
            qmf = [0 -.25 .75 .75 -.25] .* sqr2;
        elseif Par(2) == 3,
            qmf = [0 3 -9 -7 45 45 -7 -9 3] .* sqr2/64;
        elseif Par(2) == 5,
            qmf = [0 -5 15 19 -97 -26 350 350 -26 -97 19 15 -5] .* sqr2/512;
        elseif Par(2) == 7,
            qmf = [0 35 -105 -195 865 363 -3489 -307 11025 11025 -307 -3489 363 865 -195 -105 35] .* sqr2/16384;
        elseif Par(2) == 9,
            qmf = [0 -63 189 469 -1911 -1308 9188 1140 -29676 190 87318 87318 190 -29676 ...
                1140 9188 -1308 -1911 469 189 -63] .* sqr2/131072;
        end
    elseif Par(1)==4,
        dqmf = [.026748757411 -.016864118443 -.078223266529 .266864118443 .602949018236 ...
            .266864118443 -.078223266529 -.016864118443 .026748757411] .*sqr2;
        if Par(2) == 4,
            qmf = [0 -.045635881557 -.028771763114 .295635881557 .557543526229 ...
                .295635881557 -.028771763114 -.045635881557 0] .*sqr2;
        end
    end

elseif strcmp(Type,'Villasenor'),
    if Par == 1,
        % The "7-9 filters"
        qmf = [.037828455506995 -.023849465019380 -.11062440441842 .37740285561265];
        qmf = [qmf .85269867900940 reverse(qmf)];
        dqmf = [-.064538882628938 -.040689417609558 .41809227322221];
        dqmf = [dqmf .78848561640566 reverse(dqmf)];
    elseif Par == 2,
        qmf  = [-.008473 .003759 .047282 -.033475 -.068867 .383269 .767245 .383269 -.068867...
            -.033475 .047282 .003759 -.008473];
        dqmf = [0.014182  0.006292 -0.108737 -0.069163 0.448109 .832848 .448109 -.069163 -.108737 .006292 .014182];
    elseif Par == 3,
        qmf  = [0 -.129078 .047699 .788486 .788486 .047699 -.129078];
        dqmf = [0 .018914 .006989 -.067237 .133389 .615051 .615051 .133389 -.067237 .006989 .018914];
    elseif Par == 4,
        qmf  = [-1 2 6 2 -1] / (4*sqr2);
        dqmf = [1 2 1] / (2*sqr2);
    elseif Par == 5,
        qmf  = [0 1 1]/sqr2;
        dqmf = [0 -1 1 8 8 1 -1]/(8*sqr2);
    end
end

%
% Copyright (c) 1995. Jonathan Buckheit, Shaobing Chen and David Donoho
%
% Modified by Maureen Clerc and Jerome Kalifa, 1997
% [e-mail], [e-mail]
%


%
% Part of WaveLab Version 802
% Built Sunday, October 3, 1999 8:52:27 AM
% This is Copyrighted Material
% For Copying permissions see COPYING.m
% Comments? e-mail [e-mail]
%

function wcoef = FWT_SBS(x,L,qmf,dqmf)
% FWT_SBS -- Forward Wavelet Transform (symmetric extension, biorthogonal, symmetric)
%  Usage
%    wc = FWT_SBS(x,L,qmf,dqmf)
%  Inputs
%    x    1-d signal; arbitrary length
%    L    Coarsest Level of V_0;  L << J
%    qmf    quadrature mirror filter (symmetric)
%    dqmf   quadrature mirror filter (symmetric, dual of qmf)
%  Outputs
%    wc    1-d wavelet transform of x.
%
%  Description
%    1. qmf filter may be obtained from MakePBSFilter
%    2. usually, length(qmf) < 2^(L+1)
%    3. To reconstruct use IWT_SBS
%
%  See Also
%    IWT_SBS, MakePBSFilter
%
%  References
%   Based on the algorithm developed by Christopher Brislawn.
%   See "Classification of Symmetric Wavelet Transforms"
%

[n,J] = dyadlength(x);

wcoef = zeros(1,n);
beta = ShapeAsRow(x);  % take samples at finest scale as beta-coeffts

dp = dyadpartition(n);

for j=J-1:-1:L,
    [beta, alfa] = DownDyad_SBS(beta,qmf,dqmf);
    dyadj = (dp(j+1)+1):dp(j+2);
    wcoef(dyadj) = alfa;
end
wcoef(1:length(beta)) = beta;
wcoef = ShapeLike(wcoef,x);

%
% Copyright (c) 1996. Thomas P.Y. Yu
%



%
% Part of WaveLab Version 802
% Built Sunday, October 3, 1999 8:52:27 AM
% This is Copyrighted Material
% For Copying permissions see COPYING.m
% Comments? e-mail [e-mail]
%

function dp = dyadpartition(n)
% dyadpartition -- determine dyadic partition in wavelet transform of
%                  nondyadic signals

J = ceil(log2(n));

m = n;
for j=J-1:-1:0;
    if rem(m,2)==0,
        dps(j+1) = m/2;
        m = m/2;
    else
        dps(j+1) = (m-1)/2;
        m = (m+1)/2;
    end
end

dp = cumsum([1 dps]);

%
% Copyright (c) 1996. Thomas P.Y. Yu
%


%
% Part of WaveLab Version 802
% Built Sunday, October 3, 1999 8:52:27 AM
% This is Copyrighted Material
% For Copying permissions see COPYING.m
% Comments? e-mail [e-mail]
%

function [beta, alpha] = DownDyad_SBS(x,qmf,dqmf)
% DownDyad_SBS -- Symmetric Downsampling operator
%  Usage
%    [beta,alpha] = DownDyad_SBS(x,qmf,dqmf)
%  Inputs
%    x     1-d signal at fine scale
%    qmf   quadrature mirror filter
%    dqmf  dual quadrature mirror filter
%  Outputs
%    beta  coarse coefficients
%    alpha fine coefficients
%  See Also
%    FWT_SBS
%

% oddf = (rem(length(qmf),2)==1);
oddf = ~(qmf(1)==0 & qmf(length(qmf))~=0);
oddx = (rem(length(x),2)==1);

% symmetric extension of x
if oddf,
    ex = extend(x,1,1);
else
    ex = extend(x,2,2);
end

% convolution
ebeta = DownDyadLo_PBS(ex,qmf);
ealpha = DownDyadHi_PBS(ex,dqmf);

% project
if oddx,
    beta = ebeta(1:(length(x)+1)/2);
    alpha = ealpha(1:(length(x)-1)/2);
else
    beta = ebeta(1:length(x)/2);
    alpha = ealpha(1:length(x)/2);
end

%
% Copyright (c) 1996. Thomas P.Y. Yu
%


%
% Part of WaveLab Version 802
% Built Sunday, October 3, 1999 8:52:27 AM
% This is Copyrighted Material
% For Copying permissions see COPYING.m
% Comments? e-mail [e-mail]
%

function y = extend(x, par1, par2)
% extend -- perform various kinds of symmetric extension
%
if par1==1 & par2==1,
    y = [x x((length(x)-1):-1:2)];
elseif par1==1 & par2==2,
    y = [x x((length(x)-1):-1:1)];
elseif par1==2 & par2==1,
    y = [x x(length(x):-1:2)];
elseif par1==2 & par2==2,
    y = [x reverse(x)];
end

%
% Copyright (c) 1996. Thomas P.Y. Yu
%


%
% Part of WaveLab Version 802
% Built Sunday, October 3, 1999 8:52:27 AM
% This is Copyrighted Material
% For Copying permissions see COPYING.m
% Comments? e-mail [e-mail]
%

function d = DownDyadLo_PBS(x,qmf)
% DownDyadLo_PBS -- Lo-Pass Downsampling operator (periodized,symmetric)
%  Usage
%    d = DownDyadLo_PBS(x,sf)
%  Inputs
%    x    1-d signal at fine scale
%    sf   symmetric filter
%  Outputs
%    y    1-d signal at coarse scale
%
%  See Also
%    DownDyadHi_PBS, UpDyadHi_PBS, UpDyadLo_PBS, FWT_PBSi, symm_aconv
%
d = symm_aconv(qmf,x);
n = length(d);
d = d(1:2:(n-1));

%
% Copyright (c) 1995. David L. Donoho
%




%
% Part of WaveLab Version 802
% Built Sunday, October 3, 1999 8:52:27 AM
% This is Copyrighted Material
% For Copying permissions see COPYING.m
% Comments? e-mail [e-mail]
%

function y = symm_aconv(sf,x)
% symm_aconv -- Symmetric Convolution Tool for Two-Scale Transform
%  Usage
%    y = symm_aconv(sf,x)
%  Inputs
%    sf   symmetric filter
%    x    1-d signal
%  Outputs
%    y    filtered result
%
%  Description
%    Filtering by periodic convolution of x with the
%    time-reverse of sf.
%
%  See Also
%    symm_iconv, UpDyadHi_PBS, UpDyadLo_PBS, DownDyadHi_PBS, DownDyadLo_PBS
%

n = length(x);
p = length(sf);
if p < n,
    xpadded = [x x(1:p)];
else
    z = zeros(1,p);
    for i=1:p,
        imod = 1 + rem(i-1,n);
        z(i) = x(imod);
    end
    xpadded = [x z];
end

fflip = reverse(sf);
ypadded = filter(fflip,1,xpadded);

if p < n,
    y = [ypadded((n+1):(n+p)) ypadded((p+1):(n))];
else
    for i=1:n,
        imod = 1 + rem(p+i-1,n);
        y(imod) = ypadded(p+i);
    end
end

shift = (p-1)/ 2 ;
shift = 1 + rem(shift-1, n);
y = [y((1+shift):n) y(1:(shift))] ;

%
% Copyright (c) 1995. Shaobing Chen and David L. Donoho
%




%
% Part of WaveLab Version 802
% Built Sunday, October 3, 1999 8:52:27 AM
% This is Copyrighted Material
% For Copying permissions see COPYING.m
% Comments? e-mail [e-mail]
%

function d = DownDyadHi_PBS(x,qmf)
% DownDyadHi_PBS -- Hi-Pass Downsampling operator (periodized,symmetric)
%  Usage
%    d = DownDyadHi_PBS(x,sqmf)
%  Inputs