function [a,d,xout,err]=wtlift(x,tipe,N,M,intel) % [a,d,xout,err]=WTLIFT(x,tipe,N,M,intel) % % Fast discrete wavelet transform by (integer) lifting. % SINGLE iteration of the filter bank! % % INPUT: % % x Input signal % tipe Type of wavelet used ('CDF') % N Number of vanishing moments (primal step) % M Number of vanishing moments (dual step) % intel 1 for INTEGER lifting % 0 for REGULAR lifting % % OUTPUT: % % a Cell array with scaling coefficients % d Cell array with wavelet coefficients % xout Reconstructed signal after forward/inverse transform pair % err RMSE of the reconstructed signal % % EXAMPLE: % % x=real(dopnoise(1024,200,60,10,70,128)); % [a,d,xout,err]=wtlift(x+rand(size(x))+1); % % See also LIFTCO, FBANK, ABANK, SBANK. % % % Last modified by fjsimons-at-alum.mit.edu, 05/22/2012 t0=clock; defval('intel',1) defval('tipe','CDF') defval('N',2) defval('M',4) defval('x',[0 0 0 0 1 -1 2 5 1 7 0 0 0 0 0 0]') x=x(:); % Find wavelet coefficients in the catalog [h0,f0,Pa,Ua,Kp,Ku]=wc(tipe,[N M]); % In the case of CDF [2,4] we get: % Pa=[1/2 1/2]; % Ua=[-3 19 19 -3]/64; % Kp=sqrt(2) % Ku=1/sqrt(2); % and we need no more from WC. % If using integer lifting, must have scale factors as lifting steps if intel==1 & ~strcmp(tipe,'CDFI') disp('No precomputed lifting steps stored; attempt to reconstruct') % Here we now need to absorb the scaling as four more lifting steps % If indeed Ku equals 1/Kp, see Daubechies and Sweldens, 1998, Delft % book p 145 % disp('Extra lifting steps, may use CDFI instead') Pac=Pa; Uac=Ua; clear Pa Ua Pa{1}=Pac; Ua{1}=Uac; Pa{2}=-1; Pa{3}=Ku; Ua{2}=Kp-1; Ua{3}=Kp-Kp^2; % Verify that this whole thing works here: difer([Kp 0 ; 0 Ku]-[1 Ua{3} ; 0 1]*[1 0 ; -Pa{3} 1]*... [1 Ua{2} ; 0 1]*[1 0 ; -Pa{2} 1]); Ku=1; Kp=1; end % Only for interior points % SPLIT ----------------------------------------- d=x(even(x)); a=x(~even(x)); % Number of lifting steps------------------------ M=1; if iscell(Pa); M=length(Pa); end % FORWARD TRANSFORM % Loop over lifting steps------------------------ for index=1:M if iscell(Pa) P=Pa{index}; U=Ua{index}; else P=Pa; U=Ua; end Pl=length(P); Ul=length(U); % PREDICT --------------------------------------- for l=ceil(Pl/2):ceil(length(x)/2)-floor(Pl/2)-(Pl==1)*mod(length(x),2) Lp=l+[(1-ceil(Pl/2)):1:floor(Pl/2)]; if intel==1; d(l)=d(l)-floor(P(:)'*a(Lp)+1/2); end if intel==0; d(l)=d(l)-P(:)'*a(Lp); end end % UPDATE ---------------------------------------- for l=1+floor(Ul/2):floor(length(x)/2)-ceil(Ul/2)+1 Lu=l-[floor(Ul/2):-1:(1-ceil(Ul/2))]; if intel==1; a(l)=a(l)+floor(U(:)'*d(Lu)+1/2); end if intel==0; a(l)=a(l)+U(:)'*d(Lu); end end end % SCALE ----------------------------------------- d=d*Ku; a=a*Kp; % Try the inverse transform by running it backwards % UNDO SCALE ------------------------------------ y0=a/Kp; y1=d/Ku; if nargout>=2 % INVERSE TRANSFORM % Reverse loop over lifting steps------------------- for index=M:-1:1 if iscell(Pa) P=Pa{index}; U=Ua{index}; else P=Pa; U=Ua; end Pl=length(P); Ul=length(U); % UNDO UPDATE ----------------------------------- % floor(length(x)/2) for l=1+floor(Ul/2):floor(length(x)/2)-ceil(Ul/2)+1 Lu=l-[floor(Ul/2):-1:(1-ceil(Ul/2))]; if intel==1; y0(l)=y0(l)-floor(U(:)'*y1(Lu)+1/2); end if intel==0; y0(l)=y0(l)-U(:)'*y1(Lu); end end % UNDO PREDICT ---------------------------------- % ceil(length(x)/2) % ceil(Pl/2):ceil(length(x)/2)-floor(Pl/2)-(Pl==1)*mod(length(x),2) for l=ceil(Pl/2):ceil(length(x)/2)-floor(Pl/2)-(Pl==1)*mod(length(x),2) Lp=l+[(1-ceil(Pl/2)):floor(Pl/2)]; if intel==1; y1(l)=y1(l)+floor(P(:)'*y0(Lp)+1/2); end if intel==0; y1(l)=y1(l)+P(:)'*y0(Lp); end end end % MERGE ----------------------------------------- xout(~even(x),:)=y0; xout(even(x),:)=y1; err=sqrt(mean((xout(:)-x(:)).^2)); if err>1e-10 warning(sprintf('Reconstruction failed with %8.2e',err)) else % disp(sprintf('WTLIFT Reconstruction error %8.4e',err)) end % disp(sprintf('WTLIFT (Analysis & Synthesis) took %8.4f s',... % etime(clock,t0))) end