enoint.m

matlab程序

function Y = enoint(X, Degree, Dim)
%ENOINT  Essentially Non-Oscillatory interpolation.
%  Y = ENOINT(X, N) performs Nth-degree polynomial interpolation on X
%  with nonhierarchical ENO stencil selection to do reconstruction in
%  one-dimensional ENO.  If X is a matrix, interpolation is done along
%  the columns.
%
%  ENOINT(X, N, DIM) performs interpolation over dimension DIM.
%
%  See also ENOPV, ENOCA.

% Pascal Getreuer 2005

if nargin < 2, error('Not enough input arguments.'); end
if nargin < 3, 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 = size(X);
% Create polynomial singularity filter
PolyFilt = 1;

for k = 1:Degree
   PolyFilt = conv(PolyFilt,[0.5,-0.5]);   
end

% Create interpolation filters
A = zeros(Degree+1);
b = zeros(Degree+1,Degree);

for k = 0:Degree
   A(k+1,:) = ((0:Degree)/Degree).^k;
   b(k+1,:) = ((0.5:1:Degree)/Degree).^k;
end

IntFilt = A\b;
IntFilt = IntFilt - (A\(A*IntFilt - b));

% Determine stencil shifts (nonhierarchical selection)
pX = abs(filter(PolyFilt,1, ...
   [ones(Degree-1,N(2))*nan;X;ones(Degree-1,N(2))*nan]));
pX = pX(Degree+1:size(pX,1),:);
Shift = zeros(N(1)-1,N(2));
i = 0:Degree-1;

for k = 1:N(1)-1
   [tmp,Shift(k,:)] = min(pX(k + i,:),[],1);
end

% Interpolate
Y(1:2:N(1)*2-1,:) = X;
tmp = zeros(N(1)-1,N(2));

for k = 1:Degree
   i = find(Shift == k);
   iX = filter(IntFilt(:,k),1,X,[],1);
   iX = iX(2:N(1),:);
   tmp(i) = iX(i + k - 1);
end

Y(2:2:end,:) = tmp;

XSize(Dim) = N(1)*2-1;
Y = ipermute(reshape(Y,XSize(Perm)),Perm);