perform_wavelet_transf.m.svn-base

signal procesing toolbox

function x = perform_wavelet_transf(x, Jmin, dir, options)

% perform_wavelet_transf - peform fast lifting transform
%
%   y = perform_wavelet_transf(x, Jmin, dir, options);
%
%   Implement 1D and 2D symmetric wavelets with symmetric boundary treatements, using
%   a lifting implementation.
%
%   h = options.filter gives the coefficients of the lifting filter.
%   You can use h='linear' or h='7-9' to select automatically biorthogonal
%   transform with 2 and 4 vanishing moments.
%
%   You can set options.ti=1 to compute a translation invariant wavelet
%   transform.
%
%   You can set options.separable=1 to compute a separable 2D wavelet
%   transform.
%
%   Copyright (c) 2008 Gabriel Peyre

options.null = 0;
h = getoptions(options, 'filter', '9-7');
if isstr(h)
	% P/U/P/U/etc the last coefficient is scaling
   switch lower(h)
       case {'linear' '5-3'}
           h = [1/2 1/4 sqrt(2)];
       case {'9-7' '7-9'}
           h = [1.586134342 -.05298011854 -.8829110762 .4435068522 1.149604398];
       otherwise 
           error('Unknown filter');
   end
end

% detect dimensionality 
% d = getoptions(options,'ndims',-1);
d = -1;
if isfield(options, 'ndims')
    ndims = options.ndims;
end
if d<0
    if size(x,1)==1 || size(x,2)==1
        d = 1;
    else
        d = 2;
    end
end

if d==1 && size(x,1)==1
    warning('For 1D transform, the vector should be a column vector.');
    x = x(:);
end

ti = getoptions(options, 'ti', 0);

separable = getoptions(options, 'separable', 0);
if d==2 && separable==1
    options.ti = 0;
    if ti==1
        warning('Separable does not works for translation invariant transform');
    end
    % perform a separable wavelet transform
    n = size(x,1);
    if dir==1
        for i=1:n
            x(:,i) = perform_wavelet_transf(x(:,i),Jmin,dir,options);
        end
        for i=1:n
            x(i,:) = perform_wavelet_transf(x(i,:)',Jmin,dir,options)';
        end
    else
        for i=1:n
            x(i,:) = perform_wavelet_transf(x(i,:)',Jmin,dir,options)';
        end
        for i=1:n
            x(:,i) = perform_wavelet_transf(x(:,i),Jmin,dir,options);
        end
    end
    return;
end


% number of lifting steps
n = size(x,1);
m = (length(h)-1)/2;
Jmax = log2(n)-1;
jlist = Jmax:-1:Jmin;
if dir==-1
    jlist = jlist(end:-1:1);
end

if ti==0
    % subsampled
    for j=jlist
        if d==1
            x(1:2^(j+1),:) = lifting_step( x(1:2^(j+1),:), h, dir);
        else
            x(1:2^(j+1),1:2^(j+1)) = lifting_step( x(1:2^(j+1),1:2^(j+1)), h, dir);
            x(1:2^(j+1),1:2^(j+1)) = lifting_step( x(1:2^(j+1),1:2^(j+1))', h, dir)';
        end
    end
else
    % TI
    nJ = Jmax-Jmin+1;
    if dir==1 && d==1
        x = repmat(x,[1 1 nJ+1]);
    elseif dir==1 && d==2
        x = repmat(x,[1 1 3*nJ+1]);
    end
    for j=jlist
        dist = 2^(Jmax-j);
        if d==1
            if dir==1
                x(:,:,[1 j-Jmin+2]) = lifting_step_ti( x(:,:,1), h, dir, dist);
            else
                x(:,:,1) = lifting_step_ti( x(:,:,[1 j-Jmin+2]), h, dir, dist);                
            end
        else
            dj = 3*(j-Jmin);
            if dir==1
                x(:,:,[1 dj+2]) = lifting_step_ti( x(:,:,1), h, dir, dist);
                
                x(:,:,[1 dj+3]) = lifting_step_ti( x(:,:,1)', h, dir, dist);
                x(:,:,1) = x(:,:,1)'; x(:,:,dj+3) = x(:,:,dj+3)';
                
                x(:,:,[dj+[2 4]]) = lifting_step_ti( x(:,:,dj+2)', h, dir, dist);
                x(:,:,dj+2) = x(:,:,dj+2)';
                x(:,:,dj+4) = x(:,:,dj+4)';                
            else                
                x(:,:,dj+2) = x(:,:,dj+2)';
                x(:,:,dj+4) = x(:,:,dj+4)';
                x(:,:,dj+2) = lifting_step_ti( x(:,:,[dj+[2 4]]), h, dir, dist)';

                x(:,:,1) = x(:,:,1)'; 
                x(:,:,dj+3) = x(:,:,dj+3)';
                x(:,:,1) = lifting_step_ti( x(:,:,[1 dj+3]), h, dir, dist)';
                
                x(:,:,1) = lifting_step_ti( x(:,:,[1 dj+2]), h, dir, dist);
           end            
        end
    end
    if dir==-1
        x = x(:,:,1);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function x = lifting_step(x, h, dir)

% number of lifting steps
m = (length(h)-1)/2;

if dir==1
    % split
    d = x(2:2:end,:);
    x = x(1:2:end,:);
    for i=1:m
        d = d - h(2*i-1) * ( x + x([2:end end],:) );
        x = x + h(2*i  ) * ( d + d([1 1:end-1],:) );
    end
    x = [x*h(end);d/h(end)];
else
    % retrieve detail coefs
    d = x(end/2+1:end,:)*h(end);
    x = x(1:end/2,:)/h(end);
    for i=m:-1:1
        x = x - h(2*i  ) * ( d + d([1 1:end-1],:) );
        d = d + h(2*i-1) * ( x + x([2:end end],:) );
    end
    % merge
    x1 = zeros(size(x,1)*2,size(x,2));
    x1(1:2:end,:) = x;
    x1(2:2:end,:) = d;
    x = x1;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function x = lifting_step_ti(x, h, dir, dist)

% number of lifting steps
m = (length(h)-1)/2;
n = size(x,1);

s1 = (1:n)+dist;
s2 = (1:n)-dist;
% boundary conditions
s1(s1>n) = 2*n-s1(s1>n); s1(s1<1) = 2-s1(s1<1);
s2(s2>n) = 2*n-s2(s2>n); s2(s2<1) = 2-s2(s2<1);

% s1 = [2 1:n-1]; s2 = [2:n n-1];

if dir==1
    % split
    d = x;
    for i=1:m
        d = d - h(2*i-1) * ( x(s1,:) + x(s2,:) );
        x = x + h(2*i  ) * ( d(s1,:) + d(s2,:) );
    end
    x = cat(3,x*h(end),d/h(end));
else
    % retrieve detail coefs
    d = x(:,:,2)*h(end);
    x = x(:,:,1)/h(end);
    for i=m:-1:1
        x = x - h(2*i  ) * ( d(s1,:) + d(s2,:) );
        d = d + h(2*i-1) * ( x(s1,:) + x(s2,:) );
    end
    % merge
    x = (x+d)/2;
end