% dwt (Discrete wavelet transform) by Dan Piehl (2006)
%
% This is a discrete wavelet tranform based on the Haar wavelet.
% It can do one or two dimensional transforms, and requires the
% specification of a level parameter.  If a negative level
% is supplied, the wavelet inverse is performed.
% (For instance, use +8 and -8 for a 256 point signal)
%
% Note: To filter color images (here x is a 3-plane double array input)
%
% First convert to wavelets   wx=dwt(x,8);
% Then apply soft threshold   wz=thresh(wx,0.4,1);
% Then convert back           z=dwt(wz,-8);
% Optional low-pass step      z=lowpass(z);
% Clean up invalid pixels     z=imclip(z);

function b = dwt(a,level)

  % Examine parameters and return error messages, if needed
  numDimsA = ndims(a);          % Get the number of input dimensions
  S1=size(a,1);                 % Get the number of rows
  S2=size(a,2);                 % Get the number of columns
  S3=size(a,3);                 % Get the number of color planes
  if (nargin < 1 || level==0)
      error('MATLAB:wavelet:NoInputs',['No input arguments specified.'])
  end
  if (numDimsA<1 || numDimsA>3) % Only 1-3 dimensional arrays are allowed
      error('MATLAB:wavelet:NoInputs',['Invalid array dimensions.']) 
  end
  if (level==0 || level>32 || level-floor(level)>0)
      error('MATLAB:wavelet:NoInputs',['Invalid level number.'])
  end
  if not ( (S1==1 && S2>1) || (S1>1 && S2==1) || (S1>1 && S2>1) )
      error('MATLAB:wavelet:NoInputs',['Invalid array dimensions.'])
  end
  if not (S1==1 || (S1>1 && S1>=2^abs(level)))
      error('MATLAB:wavelet:NoInputs',['Invalid number of rows.'])
  end
  if not (S2==1 || (S2>1 && S2>=2^abs(level)))
      error('MATLAB:wavelet:NoInputs',['Invalid number of columns.'])
  end
  if (log(S1)/log(2)-floor(log(S1)/log(2))>0)
      error('MATLAB:wavelet:NoInputs',['Array dimensions must be powers of 2.'])
  end
  if (log(S2)/log(2)-floor(log(S2)/log(2))>0)
      error('MATLAB:wavelet:NoInputs',['Array dimensions must be powers of 2.'])
  end

  % Perform the forward or reverse wavelet transform
  b=a;                        % Reset return value
  iteration=0;                % Reset level counter
  if (level>0)                % For "forward" wavelet transform...            
    while ((S1>1 || S2>1) && iteration<level)  % Loop through all levels
      iteration=iteration+1;  % Increment iteration counter
      H1=(S1/2);              % Compute half of remaining array
      H2=(S2/2);
      if (size(a,1)==1)       % If it is a single row...
        t=b;                  % Process 1-D Haar
        b( 1,      1:H2,:) = (t(1,1:2:S2,:) + t(1,2:2:S2,:)) / sqrt(2);
        b( 1, (H2+1):S2,:) = (t(1,1:2:S2,:) - t(1,2:2:S2,:)) / sqrt(2);
      elseif (size(a,2)==1)   % If it is a single column...
        t=b;                  % Process 1-D Haar
        b(      1:H1, 1,:) = (t(1:2:S1,1,:) + t(2:2:S1,1,:)) / sqrt(2);
        b( (H1+1):S1, 1,:) = (t(1:2:S1,1,:) - t(2:2:S1,1,:)) / sqrt(2);
      else                    % For 2-D transformation...
        A11 = b(1:2:S1,1:2:S2,:)/2;  % Store partial arrays
        A12 = b(1:2:S1,2:2:S2,:)/2;
        A21 = b(2:2:S1,1:2:S2,:)/2;
        A22 = b(2:2:S1,2:2:S2,:)/2;
        b(      1:H1,      1:H2,:) = A11 + A12 + A21 + A22;  % Compute result
        b(      1:H1, (H2+1):S2,:) = A11 - A12 + A21 - A22;
        b( (H1+1):S1,      1:H2,:) = A11 + A12 - A21 - A22;
        b( (H1+1):S1, (H2+1):S2,:) = A11 - A12 - A21 + A22;
      end
      S1=S1/2;                % Drop down to the next level
      S2=S2/2;
    end
  elseif (level<0)            % For an inverse transform...
    for iteration=1:(-level)  % Determine the smallest scale
      S1=S1/2;
      S2=S2/2;
    end
    while (S1<size(a,1) || S2<size(a,2)) % Loop until array is rebuilt
      S1=S1*2;                % Double the scale factors
      S2=S2*2; 
      H1=(S1/2);              % Compute the half scale
      H2=(S2/2);
      if (size(a,1)==1)
        t=b;
        b(1,1:2:S2,:) = (t(1,1:H2,:) + t(1,(H2+1):S2,:)) / sqrt(2);
        b(1,2:2:S2,:) = (t(1,1:H2,:) - t(1,(H2+1):S2,:)) / sqrt(2);
      elseif (size(a,2)==1)
        t=b;
        b(1:2:S1,1,:) = (t(1:H1,1,:) + t((H2+1):S1,1,:)) / sqrt(2);
        b(2:2:S1,1,:) = (t(1:H1,1,:) - t((H2+1):S1,1,:)) / sqrt(2);
      else
        A11=b(      1:H1,      1:H2,:)/2;  % Compute parial arrays
        A12=b(      1:H1, (H2+1):S2,:)/2;
        A21=b( (H1+1):S1,      1:H2,:)/2;
        A22=b( (H1+1):S1, (H2+1):S2,:)/2;
        b(1:2:S1,1:2:S2,:) = A11 + A12 + A21 + A22; % Store results
        b(1:2:S1,2:2:S2,:) = A11 - A12 + A21 - A22;
        b(2:2:S1,1:2:S2,:) = A11 + A12 - A21 - A22;
        b(2:2:S1,2:2:S2,:) = A11 - A12 - A21 + A22;
      end
    end
  end