wavelet.m

EZW加入简单的熵编码

   Family = 'Reverse spline';
case 'lc5/11-a'
   Seq = {0,0;-[1,1]/2,1;[1,1]/4,0;[1,-1,-1,1]/32,2};   
   ScaleS = sqrt(2);
   ScaleD = -1/sqrt(2);
   Family = 'Low complexity';
case 'lc6/14'
   Seq = {0,0;-1,0;[-1,8,1]/16,1;[1,-6,0,6,-1]/16,2};   
   ScaleS = sqrt(2);
   ScaleD = -1/sqrt(2);
   Family = 'Low complexity';
case 'lc13/7-c'
   Seq = {0,0;[1,-9,-9,1]/16,2;[-1,5,5,-1]/16,1};   
   ScaleS = sqrt(2);
   ScaleD = -1/sqrt(2);
   Family = 'Low complexity';
otherwise
   Seq = {};
   return;
end

if ~isempty(findstr(lower(WaveletName),'rspline'))
   [Seq,ScaleS,ScaleD] = seqdual(Seq,ScaleS,ScaleD);
   Family = 'Reverse spline';
end

return;


function [Seq,ScaleS,ScaleD] = seqdual(Seq,ScaleS,ScaleD)
% Dual of a lifting sequence

L = size(Seq,1);

for k = 1:L
   % f'(z) = -f(z^-1)
   Seq{k,2} = -(Seq{k,2} - length(Seq{k,1}) + 1);
   Seq{k,1} = -fliplr(Seq{k,1});
end

if all(Seq{1,1} == 0)
   Seq = reshape({Seq{2:end,:}},L-1,2);
else
   [Seq{1:L+1,:}] = deal(0,Seq{1:L,1},0,Seq{1:L,2});
end

ScaleS = 1/ScaleS;
ScaleD = 1/ScaleD;
return;


function [h,g] = seq2hg(Seq,ScaleS,ScaleD,Dual)
% Find wavelet filters from lifting sequence
if Dual, [Seq,ScaleS,ScaleD] = seqdual(Seq,ScaleS,ScaleD); end
if rem(size(Seq,1),2), [Seq{size(Seq,1)+1,:}] = deal(0,0); end

h = {1,0};
g = {1,1};

for k = 1:2:size(Seq,1)
   h = lp_lift(h,g,{Seq{k,:}});
   g = lp_lift(g,h,{Seq{k+1,:}});
end

h = {ScaleS*h{1},h{2}};
g = {ScaleD*g{1},g{2}};

if Dual
   h{2} = -(h{2} - length(h{1}) + 1);
   h{1} = fliplr(h{1});

   g{2} = -(g{2} - length(g{1}) + 1);
   g{1} = fliplr(g{1});
end

return;


function a = lp_lift(a,b,c)
% a(z) = a(z) + b(z) c(z^2)

d = zeros(1,length(c{1})*2-1);
d(1:2:end) = c{1};
d = conv(b{1},d);
z = b{2}+c{2}*2;
zmax = max(a{2},z);
f = [zeros(1,zmax-a{2}),a{1},zeros(1,a{2} - length(a{1}) - z + length(d))];
i = zmax-z + (1:length(d));
f(i) = f(i) + d;

if all(abs(f) < 1e-12)
   a = {0,0};
else
   i = find(abs(f)/max(abs(f)) > 1e-10);
   i1 = min(i);
   a = {f(i1:max(i)),zmax-i1+1};
end
return;


function X = xfir(B,Z,X,Ext)
%XFIR  Noncausal FIR filtering with boundary handling.
%   Y = XFIR(B,Z,X,EXT) filters X with FIR filter B with leading
%   delay -Z along the columns of X.  EXT specifies the  boundary
%   handling.  Special handling  is done for one and two-tap filters.

% Pascal Getreuer 2005-2006

N = size(X);

% Special handling for short filters
if length(B) == 1 & Z == 0
   if B == 0
      X = zeros(size(X));
   elseif B ~= 1
      X = B*X;
   end
   return;
end

% Compute the number of samples to add to each end of the signal
pl = max(length(B)-1-Z,0);       % Padding on the left end
pr = max(Z,0);                   % Padding on the right end

switch lower(Ext)
case {'sym','wsws'}   % Symmetric extension, WSWS
   if all([pl,pr] < N(1))
         X = filter(B,1,X([pl+1:-1:2,1:N(1),N(1)-1:-1:N(1)-pr],:,:),[],1);
         X = X(Z+pl+1:Z+pl+N(1),:,:);
      return;
   else
      i = [1:N(1),N(1)-1:-1:2];
      Ns = 2*N(1) - 2 + (N(1) == 1);
      i = i([rem(pl*(Ns-1):pl*Ns-1,Ns)+1,1:N(1),rem(N(1):N(1)+pr-1,Ns)+1]);
   end
case {'symh','hshs'}  % Symmetric extension, HSHS
   if all([pl,pr] < N(1))
      i = [pl:-1:1,1:N(1),N(1):-1:N(1)-pr+1];
   else
      i = [1:N(1),N(1):-1:1];
      Ns = 2*N(1);
      i = i([rem(pl*(Ns-1):pl*Ns-1,Ns)+1,1:N(1),rem(N(1):N(1)+pr-1,Ns)+1]);
   end
case 'wshs'           % Symmetric extension, WSHS
   if all([pl,pr] < N(1))
      i = [pl+1:-1:2,1:N(1),N(1):-1:N(1)-pr+1];
   else
      i = [1:N(1),N(1):-1:2];
      Ns = 2*N(1) - 1;
      i = i([rem(pl*(Ns-1):pl*Ns-1,Ns)+1,1:N(1),rem(N(1):N(1)+pr-1,Ns)+1]);
   end
case 'hsws'           % Symmetric extension, HSWS
   if all([pl,pr] < N(1))
      i = [pl:-1:1,1:N(1),N(1)-1:-1:N(1)-pr];
   else
      i = [1:N(1),N(1)-1:-1:1];
      Ns = 2*N(1) - 1;
      i = i([rem(pl*(Ns-1):pl*Ns-1,Ns)+1,1:N(1),rem(N(1):N(1)+pr-1,Ns)+1]);
   end
case 'zpd'
   Ml = N; Ml(1) = pl;
   Mr = N; Mr(1) = pr;
   
   X = filter(B,1,[zeros(Ml);X;zeros(Mr)],[],1);
   X = X(Z+pl+1:Z+pl+N(1),:,:);
   return;
case 'per'            % Periodic
   i = [rem(pl*(N(1)-1):pl*N(1)-1,N(1))+1,1:N(1),rem(0:pr-1,N(1))+1];
case 'sp0'            % Constant extrapolation
   i = [ones(1,pl),1:N(1),N(1)+zeros(1,pr)];
case 'asym'           % Asymmetric extension
   i1 = [ones(1,pl),1:N(1),N(1)+zeros(1,pr)];

   if all([pl,pr] < N(1))
      i2 = [pl+1:-1:2,1:N(1),N(1)-1:-1:N(1)-pr];
   else
      i2 = [1:N(1),N(1)-1:-1:2];
      Ns = 2*N(1) - 2 + (N(1) == 1);
      i2 = i2([rem(pl*(Ns-1):pl*Ns-1,Ns)+1,1:N(1),rem(N(1):N(1)+pr-1,Ns)+1]);
   end
   
   X = filter(B,1,2*X(i1,:,:) - X(i2,:,:),[],1);
   X = X(Z+pl+1:Z+pl+N(1),:,:);
   return;
otherwise
   error(['Unknown boundary handling, ''',Ext,'''.']);
end

X = filter(B,1,X(i,:,:),[],1);
X = X(Z+pl+1:Z+pl+N(1),:,:);
return;


function [x1,phi,x2,psi] = cascade(h,g,dx)
% Wavelet cascade algorithm

c = h{1}*2/sum(h{1});
x = 0:dx:length(c) - 1;
x1 = x - h{2};
phi0 = 1 - abs(linspace(-1,1,length(x))).';

ii = []; jj = []; s = [];

for k = 1:length(c)
   xk = 2*x - (k-1);
   i = find(xk >= 0 & xk <= length(c) - 1);
   ii = [ii,i];
   jj = [jj,floor(xk(i)/dx)+1];
   s = [s,c(k)+zeros(size(i))];
end

% Construct a sparse linear operator that iterates the dilation equation
Dilation = sparse(ii,jj,s,length(x),length(x));

for N = 1:30
   phi = Dilation*phi0;
   if norm(phi - phi0,inf) < 1e-5, break; end
   phi0 = phi;
end

if norm(phi) == 0
   phi = ones(size(phi))*sqrt(2);   % Special case for Haar scaling function
else
   phi = phi/(norm(phi)*sqrt(dx));  % Rescale result
end

if nargout > 2
   phi2 = phi(1:2:end);  % phi2 is approximately phi(2x)

   if length(c) == 2
      L = length(phi2);
   else
      L = ceil(0.5/dx);
   end

   % Construct psi from translates of phi2
   c = g{1};
   psi = zeros(length(phi2)+L*(length(c)-1),1);
   x2 = (0:length(psi)-1)*dx - g{2} - 0*h{2}/2;

   for k = 1:length(c)
      i = (1:length(phi2)) + L*(k-1);
      psi(i) = psi(i) + c(k)*phi2;
   end
end
return;


function s = filterstr(a,K)
% Convert a filter to a string

[n,d] = rat(K/sqrt(2));

if d < 50
   a{1} = a{1}/sqrt(2);   % Scale filter by sqrt(2)
   s = '( ';
else
   s = '';
end

Scale = [pow2(1:15),10,20,160,280,inf];

for i = 1:length(Scale)
   if norm(round(a{1}*Scale(i))/Scale(i) - a{1},inf) < 1e-9
      a{1} = a{1}*Scale(i);  % Scale filter by a power of 2 or 160
      s = '( ';
      break;
   end
end

z = a{2};
LineOff = 0;

for k = 1:length(a{1})
   v = a{1}(k);

   if v ~= 0  % Only display nonzero coefficients
      if k > 1
         s2 = [' ',char(44-sign(v)),' '];
         v = abs(v);
      else
         s2 = '';
      end

      s2 = sprintf('%s%g',s2,v);

      if z == 1
         s2 = sprintf('%s z',s2);
      elseif z ~= 0
         s2 = sprintf('%s z^%d',s2,z);
      end

      if length(s) + length(s2) > 72 + LineOff  % Wrap long lines
         s2 = [char(10),'        ',s2];
         LineOff = length(s);
      end

      s = [s,s2];
   end

   z = z - 1;
end

if s(1) == '('
   s = [s,' )'];

   if d < 50, s = [s,' sqrt(2)']; end

   if i < length(Scale)
      s = sprintf('%s/%d',s,Scale(i));
   end
end

return;

function N = numvanish(g)
% Determine the number of vanishing moments from highpass filter g(z)

for N = 0:length(g)-1  % Count the number of roots at z = 1
   [g,r] = deconv(g,[1,-1]);
   if norm(r,inf) > 1e-7, break; end
end
return;