perform_wavelet_transform.m

是一本非常不错的书

%    x    1-d signal at fine scale
%    sqmf symmetric filter
%  Outputs
%    y    1-d signal at coarse scale
%
%  See Also
%    DownDyadLo_PBS, UpDyadHi_PBS, UpDyadLo_PBS, FWT_PBS, symm_iconv
%
d = symm_iconv( MirrorSymmFilt(qmf),lshift(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 = MirrorSymmFilt(x)
% MirrorSymmFilt -- apply (-1)^t modulation to symmetric filter
%  Usage
%    h = MirrorSymmFilt(l)
%  Inputs
%    l   symmetric filter
%  Outputs
%    h   symmetric filter with DC frequency content shifted
%        to Nyquist frequency
%
%  Description
%    h(t) = (-1)^t  * x(t),  -k <= t <= k ; length(x)=2k+1
%
%  See Also
%    DownDyadHi_PBS
%
k = (length(x)-1)/2;
y = ( (-1).^((-k):k) ) .*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 = symm_iconv(sf,x)
% symm_iconv -- Symmetric Convolution Tool for Two-Scale Transform
%  Usage
%    y = iconv(sf,x)
%  Inputs
%    sf  symmetric filter
%    x   1-d signal
%  Output
%    y    filtered result
%
%  Description
%    Filtering by periodic convolution of x with sf
%
%  See Also
%    symm_aconv, UpDyadHi_PBS, UpDyadLo_PBS, DownDyadHi_PBS, DownDyadLo_PBS
%
n = length(x);
p = length(sf);
if p <= n,
    xpadded = [x((n+1-p):n) x];
else
    z = zeros(1,p);
    for i=1:p,
        imod = 1 + rem(p*n -p + i-1,n);
        z(i) = x(imod);
    end
    xpadded = [z x];
end
ypadded = filter(sf,1,xpadded);
y = ypadded((p+1):(n+p));

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



%
% 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 x = IWT_SBS(wc,L,qmf,dqmf)
% iwt_po -- Inverse Wavelet Transform (symmetric extension, biorthogonal, symmetric)
%  Usage
%    x = IWT_SBS(wc,L,qmf,dqmf)
%  Inputs
%    wc     1-d wavelet transform: length(wc)= 2^J.
%    L      Coarsest scale (2^(-L) = scale of V_0); L << J;
%    qmf     quadrature mirror filter
%    dqmf    dual quadrature mirror filter (symmetric, dual of qmf)
%  Outputs
%    x      1-d signal reconstructed from wc
%  Description
%    Suppose wc = FWT_SBS(x,L,qmf,dqmf) where qmf and dqmf are orthonormal
%    quad. mirror filters made by MakeBioFilter.  Then x can be reconstructed
%    by
%      x = IWT_SBS(wc,L,qmf,dqmf)
%  See Also:
%    FWT_SBS, MakeBioFilter
%

wcoef = ShapeAsRow(wc);
[n,J] = dyadlength(wcoef);

dp = dyadpartition(n);

x = wcoef(1:dp(L+1));

for j=L:J-1,
    dyadj = (dp(j+1)+1):dp(j+2);
    x = UpDyad_SBS(x, wcoef(dyadj), qmf, dqmf);
end
x = ShapeLike(x,wc);

%
% 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 x = UpDyad_SBS(beta,alpha,qmf,dqmf)
% UpDyad_SBS -- Symmetric Upsampling operator
%  Usage
%    x = UpDyad_SBS(beta,alpha,qmf,dqmf)
%  Inputs
%    beta  coarse coefficients
%    alpha fine coefficients
%    qmf   quadrature mirror filter
%    dqmf  dual quadrature mirror filter
%  Outputs
%    x     1-d signal at fine scale
%  See Also
%    DownDyad_SBS, IWT_SBS
%

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

L = length(beta)+length(alpha);

if oddf,
    if oddx,
        ebeta = extend(beta,1,1);
        ealpha = extend(alpha,2,2);
    else
        ebeta = extend(beta,2,1);
        ealpha = extend(alpha,1,2);
    end
else
    if oddx,
        ebeta = extend(beta,1,2);
        ealpha = [alpha 0 -reverse(alpha)];
    else
        ebeta = extend(beta,2,2);
        ealpha = [alpha -reverse(alpha)];
    end
end

coarse = UpDyadLo_PBS(ebeta,dqmf);
fine = UpDyadHi_PBS(ealpha,qmf);
x = coarse + fine;
x = x(1:L);

%
% 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 = UpDyadLo_PBS(x,qmf)
% UpDyadLo_PBS -- Lo-Pass Upsampling operator; periodized
%  Usage
%    u = UpDyadLo_PBS(d,sf)
%  Inputs
%    d    1-d signal at coarser scale
%    sf   symmetric filter
%  Outputs
%    u    1-d signal at finer scale
%
%  See Also
%    DownDyadLo_PBS , DownDyadHi_PBS , UpDyadHi_PBS, IWT_PBS, symm_iconv
%
y =  symm_iconv(qmf, UpSample(x,2) );

%
% 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 = UpDyadHi_PBS(x,qmf)
% UpDyadHi_PBS -- Hi-Pass Upsampling operator; periodized
%  Usage
%    u = UpDyadHi_PBS(d,f)
%  Inputs
%    d    1-d signal at coarser scale
%    sf   symmetric filter
%  Outputs
%    u    1-d signal at finer scale
%
%  See Also
%    DownDyadLo_PBS, DownDyadHi_PBS, UpDyadLo_PBS, IWT_PBS, symm_aconv
%

y = symm_aconv( MirrorSymmFilt(qmf), rshift( UpSample(x,2) ) );

%
% 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 wc = FWT2_SBS(x,L,qmf,dqmf)
% FWT2_SBS -- 2-dimensional wavelet transform
%              (symmetric extension, bi-orthogonal)
%  Usage
%    wc = FWT2_SBS(x,L,qmf,dqmf)
%  Inputs
%    x     2-d image (n by n array, n arbitrary)
%    L     coarsest level
%    qmf   low-pass quadrature mirror filter
%    dqmf  high-pass dual quadrature mirror filter
%  Output
%    wc    2-d wavelet transform
%  Description
%    A two-dimensional Wavelet Transform is computed for the
%    matrix x. To reconstruct, use IWT2_SBS.
%  See Also
%    IWT2_SBS
%

[m,J] = dyadlength(x(:,1));
[n,K] = dyadlength(x(1,:));
wc = x;
mc = m;
nc = n;

J = min([J,K]);

for jscal=J-1:-1:L,

    if rem(mc,2)==0,
        top = (mc/2+1):mc;
        bot = 1:(mc/2);
    else
        top = ((mc+1)/2+1):mc;
        bot = 1:((mc+1)/2);
    end
    if rem(nc,2)==0,
        right = (nc/2+1):nc;
        left = 1:(nc/2);
    else
        right = ((nc+1)/2+1):nc;
        left = 1:((nc+1)/2);
    end

    for ix=1:mc,
        row = wc(ix,1:nc);
        [beta,alpha] = DownDyad_SBS(row,qmf,dqmf);
        wc(ix,left) = beta;
        wc(ix,right) = alpha;
    end
    for iy=1:nc,
        column = wc(1:mc,iy)';
        [beta,alpha] = DownDyad_SBS(column,qmf,dqmf);
        wc(bot,iy) = beta';
        wc(top,iy) = alpha';
    end
    mc = bot(length(bot));
    nc = left(length(left));
end







%
% 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 x = IWT2_SBS(wc,L,qmf,dqmf)
% IWT2_SBS -- Inverse 2d Wavelet Transform
%            (symmetric extention, bi-orthogonal)
%  Usage
%    x = IWT2_SBS(wc,L,qmf,dqmf)
%  Inputs
%      wc    2-d wavelet transform [n by n array, n arbitrary]
%      L     coarse level
%      qmf   low-pass quadrature mirror filter
%      dqmf  high-pas dual quadrature mirror filter
%  Outputs
%      x     2-d signal reconstructed from wc
%  Description
%      If wc is the result of a forward 2d wavelet transform, with
%           wc = FWT2_SBS(x,L,qmf,dqmf)
%      then x = IWT2_SBS(wc,L,qmf,dqmf) reconstructs x exactly if qmf is a nice
%      quadrature mirror filter, e.g. one made by MakeBioFilter
%  See Also:
%    FWT2_SBS, MakeBioFilter
%

[m,J] = dyadlength(wc(:,1));
[n,K] = dyadlength(wc(1,:));
% assume m==n, J==K

x = wc;

dpm = dyadpartition(m);

for jscal=L:J-1,
    bot = 1:dpm(jscal+1);
    top = (dpm(jscal+1)+1):dpm(jscal+2);
    all = [bot top];

    nc = length(all);

    for iy=1:nc,
        x(all,iy) =  UpDyad_SBS(x(bot,iy)', x(top,iy)', qmf, dqmf)';
    end
    for ix=1:nc,
        x(ix,all) = UpDyad_SBS(x(ix,bot), x(ix,top), qmf, dqmf);
    end
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 wc = FWT2_PO(x,L,qmf)
% FWT2_PO -- 2-d MRA wavelet transform (periodized, orthogonal)
%  Usage
%    wc = FWT2_PO(x,L,qmf)
%  Inputs
%    x     2-d image (n by n array, n dyadic)
%    L     coarse level
%    qmf   quadrature mirror filter
%  Outputs
%    wc    2-d wavelet transform
%
%  Description
%    A two-dimensional Wavelet Transform is computed for the
%    array x.  To reconstruct, use IWT2_PO.
%
%  See Also
%    IWT2_PO, MakeONFilter
%
[n,J] = quadlength(x);
wc = x;
nc = n;
for jscal=J-1:-1:L,
    top = (nc/2+1):nc; bot = 1:(nc/2);
    for ix=1:nc,
        row = wc(ix,1:nc);
        wc(ix,bot) = DownDyadLo(row,qmf);
        wc(ix,top) = DownDyadHi(row,qmf);
    end
    for iy=1:nc,
        row = wc(1:nc,iy)';
        wc(top,iy) = DownDyadHi(row,qmf)';
        wc(bot,iy) = DownDyadLo(row,qmf)';
    end
    nc = nc/2;
end

%
% Copyright (c) 1993. 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 x = IWT2_PO(wc,L,qmf)
% IWT2_PO -- Inverse 2-d MRA wavelet transform (periodized, orthogonal)
%  Usage
%    x = IWT2_PO(wc,L,qmf)
%  Inputs
%    wc    2-d wavelet transform [n by n array, n dyadic]
%    L     coarse level
%    qmf   quadrature mirror filter
%  Outputs
%    x     2-d signal reconstructed from wc
%