enopv.m

matlab程序

function X = enopv(X,Level,Degree,Tol,Dim)
%ENOPV  Essentially Non-Oscillatory point-value decomposition.
%   Y = ENOPV(X,L,N) computes the L-stage Essentially Non-Oscillatory
%   (ENO) decomposition of signal X using point-value discretization
%   and Nth-degree interpolation.  The size of X must be such that
%   length(X) = 2^K + 1 for some integer K.  If not, the right end
%   is padded using constant extrapolation.  By default, the
%   decomposition is performed down the columns.
%
%   ENOPV(X,-L,N) is the inverse transform, reversing L levels.
%
%   For approximation, 
%   ENOPV(X,L,N,T) with T > 0 uses the modified decomposition procedure 
%   such that the maximum error between reconstruction and X is 
%   bounded by T.  If T = 0, standard decomposition is done.
%
%   ENOPV(X,L,N,T,DIM) applies the transform across the dimension DIM.
%
%   Examples:
%   Y = enopv(X,4,3);  % 4-stage decomposition with cubic interpolation
%   R = enopv(X,-4,3); % reconstruct
%
%   Y = enopv(X,4,3,0.001);  % approximately decompose with tolerance = 0.001
%   A = enopv(X,-4,3);       % synthesize the approximation
%   plot([X,A]);
%
%   Reference:
%   [1] F. Arandiga and R. Donat.  ``Nonlinear Multiscale
%       Decompositions: The Approach of A. Harten.''  Numerical
%       Algorithms 23 (2000) 175-216.
%
%   See also ENOCA, ENOPV2, ENOINT.

% Pascal Getreuer 2005

if nargin < 3, error('Not enough input arguments.'); end
if nargin < 4, Tol = 0; end
if nargin < 5, Dim = min(find(size(X) ~= 1)); end

XSize = size(X);
N = XSize(Dim);
Perm = [Dim:max(length(XSize),Dim) 1:Dim-1];
X = reshape(permute(X,Perm),N,prod(XSize)/N);
N = pow2(ceil(log2(N-1))) + 1 - N;

if N ~= 0
   warning('Signal padded to 2^K+1 length.');
   XSize(Dim) = XSize(Dim) + N;
   X = X([1:size(X,1),size(X,1)*ones(1,N)],:);
end

if Level > 0
   if Tol == 0
      for k = 1:Level
         N = (size(X,1) - 1)*pow2(1-k) + 1;
         A = X(1:2:N,:);
         R = X(1:N,:) - enoint(A,Degree);
         X(1:N,:) = [A;R(2:2:N-1,:)];
      end
   else    
      Z = X;
      Y = X(1:2^Level:size(X,1),:);
      X = Y;
      
      for k = Level:-1:1
         Y = enoint(Y,Degree);
         N = size(Y,1);         
         D = Z(1+2^(k-1):2^k:size(Z,1),:) - Y(2:2:N,:);
         D = D.*(abs(D) > Tol);
         Y(2:2:N,:) = Y(2:2:N,:) + D;
         X = [X;D];
      end
   end
else
   for k = Level:-1
      N = (size(X,1) - 1)*pow2(1+k) + 1;
      P = enoint(X(1:(N+1)/2,:),Degree);
      X(2:2:N,:) = P(2:2:N,:) + X((N+1)/2+1:N,:);
      X(1:2:N,:) = P(1:2:N,:);
   end
end

X = ipermute(reshape(X,XSize(Perm)),Perm);