enoca.m

matlab程序

function X = enoca(X,Level,Degree,Dim)
%ENOCA  Essentially Non-Oscillatory cell-average decomposition.
%   Y = ENOCA(X,L,N) computes the L-stage Essentially Non-Oscillatory
%   (ENO) decomposition of signal X using cell-average discretization
%   and Nth-degree interpolation.  The signal length must be divisible
%   by 2^L.
%
%   ENOCA(X,-L,N) is the inverse transform, reversing L levels.
%
%   ENOCA(X,L,N,DIM) applies the transform across the dimension DIM.
%
%   Example:
%   % Transform with 3 stages and 5th-degree interpolation
%   Y = enoca(X,3,5);   % Decompose X
%   R = enoca(Y,-3,5);  % Recover X from Y
%
%   Reference:
%   [1] F. Arandiga and R. Donat.  ``Nonlinear Multiscale
%       Decompositions: The Approach of A. Harten.''  Numerical
%       Algorithms 23 (2000) 175-216.
%
%   See also ENOPV, ENOCA2, ENOINT.

% Pascal Getreuer 2005

if nargin < 3, error('Not enough input arguments.'); end
if nargin < 4, 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);

if rem(N,pow2(abs(Level))), error('Invalid input size.'); end

if Level > 0
   for k = 1:Level
      N = size(X,1)*pow2(1-k);
      A = (X(1:2:N,:) + X(2:2:N,:))/2;
      P = diff(enoint(filter(2,[1,-1],[zeros(1,size(X,2));A],[],1),Degree),1);
      X(1:N,:) = [A;X(2:2:N,:) - P(2:2:N,:)];
   end
else
   for k = Level:-1
      N = size(X,1)*pow2(1+k);
      P = diff(enoint(filter(2,[1,-1],[zeros(1,size(X,2));X(1:N/2,:)],[],1),Degree),1);
      X([1:2:N,2:2:N],:) = P([1:2:N,2:2:N],:) + [-X(N/2+1:N,:);X(N/2+1:N,:)];
   end
end

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