perform_wavelet_transform.m

是一本非常不错的书

%  Usage
%    wc = FWT_PO(x,L,qmf)
%  Inputs
%    x    1-d signal; length(x) = 2^J
%    L    Coarsest Level of V_0;  L << J
%    qmf  quadrature mirror filter (orthonormal)
%  Outputs
%    wc    1-d wavelet transform of x.
%
%  Description
%    1. qmf filter may be obtained from MakeONFilter
%    2. usually, length(qmf) < 2^(L+1)
%    3. To reconstruct use IWT_PO
%
%  See Also
%    IWT_PO, MakeONFilter
%
[n,J] = dyadlength(x) ;
wcoef = zeros(1,n) ;
beta = ShapeAsRow(x);  %take samples at finest scale as beta-coeffts
for j=J-1:-1:L
    alfa = DownDyadHi(beta,qmf);
    wcoef(dyad(j)) = alfa;
    beta = DownDyadLo(beta,qmf) ;
end
wcoef(1:(2^L)) = beta;
wcoef = ShapeLike(wcoef,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 [n,J] = dyadlength(x)
% dyadlength -- Find length and dyadic length of array
%  Usage
%    [n,J] = dyadlength(x)
%  Inputs
%    x    array of length n = 2^J (hopefully)
%  Outputs
%    n    length(x)
%    J    least power of two greater than n
%
%  Side Effects
%    A warning is issued if n is not a power of 2.
%
%  See Also
%    quadlength, dyad, dyad2ix
%
n = length(x) ;
J = ceil(log(n)/log(2));
if 2^J ~= n ,
    disp('Warning in dyadlength: n != 2^J')
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 row = ShapeAsRow(sig)
% ShapeAsRow -- Make signal a row vector
%  Usage
%    row = ShapeAsRow(sig)
%  Inputs
%    sig     a row or column vector
%  Outputs
%    row     a row vector
%
%  See Also
%    ShapeLike
%
row = sig(:)';


%
% 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(x,qmf)
% DownDyadHi -- Hi-Pass Downsampling operator (periodized)
%  Usage
%    d = DownDyadHi(x,f)
%  Inputs
%    x    1-d signal at fine scale
%    f    filter
%  Outputs
%    y    1-d signal at coarse scale
%
%  See Also
%    DownDyadLo, UpDyadHi, UpDyadLo, FWT_PO, iconv
%
d = iconv( MirrorFilt(qmf),lshift(x));
n = length(d);
d = d(1:2:(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 y = MirrorFilt(x)
% MirrorFilt -- Apply (-1)^t modulation
%  Usage
%    h = MirrorFilt(l)
%  Inputs
%    l   1-d signal
%  Outputs
%    h   1-d signal with DC frequency content shifted
%        to Nyquist frequency
%
%  Description
%    h(t) = (-1)^(t-1)  * x(t),  1 <= t <= length(x)
%
%  See Also
%    DyadDownHi
%

y = -( (-1).^(1:length(x)) ).*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 = lshift(x)
% lshift -- Circular left shift of 1-d signal
%  Usage
%    l = lshift(x)
%  Inputs
%    x   1-d signal
%  Outputs
%    l   1-d signal
%        l(i) = x(i+1) except l(n) = x(1)
%

y = [ x( 2:length(x) ) x(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 y = iconv(f,x)
% iconv -- Convolution Tool for Two-Scale Transform
%  Usage
%    y = iconv(f,x)
%  Inputs
%    f   filter
%    x   1-d signal
%  Outputs
%    y   filtered result
%
%  Description
%    Filtering by periodic convolution of x with f
%
%  See Also
%    aconv, UpDyadHi, UpDyadLo, DownDyadHi, DownDyadLo
%
n = length(x);
p = length(f);
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(f,1,xpadded);
y = ypadded((p+1):(n+p));

%
% 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 i = dyad(j)
% dyad -- Index entire j-th dyad of 1-d wavelet xform
%  Usage
%    ix = dyad(j);
%  Inputs
%    j     integer
%  Outputs
%    ix    list of all indices of wavelet coeffts at j-th level
%
i = (2^(j)+1):(2^(j+1)) ;

%
% 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 d = DownDyadLo(x,qmf)
% DownDyadLo -- Lo-Pass Downsampling operator (periodized)
%  Usage
%    d = DownDyadLo(x,f)
%  Inputs
%    x    1-d signal at fine scale
%    f    filter
%  Outputs
%    y    1-d signal at coarse scale
%
%  See Also
%    DownDyadHi, UpDyadHi, UpDyadLo, FWT_PO, aconv
%
d = aconv(qmf,x);
n = length(d);
d = d(1:2:(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 y = aconv(f,x)
% aconv -- Convolution Tool for Two-Scale Transform
%  Usage
%    y = aconv(f,x)
%  Inputs
%    f    filter
%    x    1-d signal
%  Outputs
%    y    filtered result
%
%  Description
%    Filtering by periodic convolution of x with the
%    time-reverse of f.
%
%  See Also
%    iconv, UpDyadHi, UpDyadLo, DownDyadHi, DownDyadLo
%

n = length(x);
p = length(f);
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(f);
ypadded = filter(fflip,1,xpadded);
y = ypadded(p:(n+p-1));

%
% 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 vec = ShapeLike(sig,proto)
% ShapeLike -- Make 1-d signal with given shape
%  Usage
%    vec = ShapeLike(sig,proto)
%  Inputs
%    sig      a row or column vector
%    proto    a prototype shape (row or column vector)
%  Outputs
%    vec      a vector with contents taken from sig
%             and same shape as proto
%
%  See Also
%    ShapeAsRow
%
sp = size(proto);
ss = size(sig);
if( sp(1)>1 & sp(2)>1 )
    disp('Weird proto argument to ShapeLike')
elseif ss(1)>1 & ss(2) > 1,
    disp('Weird sig argument to ShapeLike')
else
    if(sp(1) > 1),
        if ss(1) > 1,
            vec = sig;
        else
            vec = sig(:);
        end
    else
        if ss(2) > 1,
            vec = sig;
        else
            vec = sig(:)';
        end
    end
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 = IWT_PO(wc,L,qmf)
% IWT_PO -- Inverse Wavelet Transform (periodized, orthogonal)
%  Usage
%    x = IWT_PO(wc,L,qmf)
%  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
%  Outputs
%    x      1-d signal reconstructed from wc
%
%  Description
%    Suppose wc = FWT_PO(x,L,qmf) where qmf is an orthonormal quad. mirror
%    filter, e.g. one made by MakeONFilter. Then x can be reconstructed by
%      x = IWT_PO(wc,L,qmf)
%
%  See Also
%    FWT_PO, MakeONFilter
%
wcoef = ShapeAsRow(wc);
x = wcoef(1:2^L);
[n,J] = dyadlength(wcoef);
for j=L:J-1
    x = UpDyadLo(x,qmf) + UpDyadHi(wcoef(dyad(j)),qmf)  ;
end
x = ShapeLike(x,wc);

%
% 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 = UpDyadLo(x,qmf)
% UpDyadLo -- Lo-Pass Upsampling operator; periodized
%  Usage
%    u = UpDyadLo(d,f)
%  Inputs
%    d    1-d signal at coarser scale
%    f    filter
%  Outputs
%    u    1-d signal at finer scale
%
%  See Also
%    DownDyadLo, DownDyadHi, UpDyadHi, IWT_PO, iconv
%
y =  iconv(qmf, 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 = UpSample(x,s)
% UpSample -- Upsampling operator
%  Usage
%    u = UpSample(d[,s])
%  Inputs
%    d   1-d signal, of length n
%    s   upsampling scale, default = 2
%  Outputs
%    u   1-d signal, of length s*n with zeros
%        interpolating alternate samples
%        u(s*i-1) = d(i), i=1,...,n
%

if nargin == 1, s = 2; end
n = length(x)*s;
y = zeros(1,n);
y(1:s:(n-s+1) )=x;