Skip to content

# csdms-contrib/slepian_echo

Fetching contributors…
Cannot retrieve contributors at this time
129 lines (121 sloc) 4.07 KB
 function f=angularD6WT(x,levels,precon,tipe,mr) % f=ANGULARD6WT(x,[n1 n2],precon,tipe,mr) % % Perform D6 (perhaps preconditioned) interval wavelet transform in the % first two indices of a three-dimensional array up to [n1 n2] scales. % % INPUT: % % x The three-dimensional array, dimensions must be powers of two % n1,n2 The number of levels in the first two directions % n1=0, n2=0 is the identity if precon=[0 0] % precon Array of length 2 identifying preconditioning % 1 Precondition (NOT an orthogonal transform) so a sequence of % ones and a linearly increasing sequence are mapped to a % sequence of zeroes in the wavelet bands % 0 Don't precondition (YES, this is an orthogonal transform) % so a sequence of ones and a linearly increasing sequence are % not mapped to all zeroes at the edges in the wavelet bands, % i.e. the first and last two coefficients are not zero % tipe 'forward'|'inverse'|'transpose'|'inversetranspose', always of % the same type in both indices % mr 0 The "original" improper, dimensionally sequential transform % 1 The "proper" multiresolution transform [default] % % OUTPUT: % % f The three-dimensional array with wavelets and scaling % functions. Put scaling coefficients in front, followed by % wavelet coefficients. % % SEE ALSO: PRECOND6, D6BOXCOF, D6BOXSTEP, D6BOXSTEPI % % EXAMPLE: % % angularD6WT('demo1') % % Inspired by Ignace Loris (igloris@vub.ac.be) on 26/06/2009 % Last modified by fjsimons-at-alum.mit.edu, 10/22/2010 if ~isstr(x) % Get the levels of the decomposition defval('levels',[3 3]) defval('tipe','forward') defval('mr',1) % Get and check the number of decompositions n1=levels(1); n2=levels(2); if mr==1 && n1~=n2 error('Number of scales must be equal in both dimensions') else n=n1; end % Get the coefficients cofs=d6boxcof; % Precondition, maybe switch tipe case {'forward','inversetranspose'} x=preconD6(x,precon,tipe,cofs); end % Prepare for output f=x; % Now do loop over levels switch tipe case {'forward','inversetranspose'} if mr==0 % In the first dimension for level=1:n1; f=d6boxstep(f,level,1,cofs); end % In the second dimension for level=1:n2; f=d6boxstep(f,level,2,cofs); end elseif mr==1 for level=1:n % Work on shorter and shorter matrices lo1=1:size(f,1)/2^(level-1); lo2=1:size(f,2)/2^(level-1); % Which makes this look like it's always at the first level f(lo1,lo2,:)=d6boxstep(f(lo1,lo2,:),1,2,cofs); f(lo1,lo2,:)=d6boxstep(f(lo1,lo2,:),1,1,cofs); end end case {'inverse','transpose'} if mr==0 % In the second dimension for level=n2:-1:1; f=d6boxstepi(f,level,2,cofs); end for level=n1:-1:1; f=d6boxstepi(f,level,1,cofs); end elseif mr==1 for level=n:-1:1 % Work on shorter and shorter matrices lo1=1:size(f,1)/2^(level-1); lo2=1:size(f,2)/2^(level-1); % Which makes this look like it's always at the first level f(lo1,lo2,:)=d6boxstepi(f(lo1,lo2,:),1,2,cofs); f(lo1,lo2,:)=d6boxstepi(f(lo1,lo2,:),1,1,cofs); end end end % Precondition, maybe switch tipe case {'inverse','transpose'} f=preconD6(f,precon,tipe,cofs); end elseif strcmp(x,'demo1') % Test Ignace's code for the non-multiresolution case whereitsat=fullfile(getenv('MFILES'),'ignaceloris'); whereiamat=fullfile(getenv('MFILES'),'wavelets'); x=rand(256,256,32); precon=[1 0]; ff=angularD6WT(x,[4 4],precon,'forward',1); fi=angularD6WT(x,[4 4],precon,'inverse',1); fit=angularD6WT(x,[4 4],precon,'inversetranspose',1); ft=angularD6WT(x,[4 4],precon,'transpose',1); rmpath(whereiamat) addpath(whereitsat) ffil=angularD6MRWT(x,[4 4],precon,'forward'); fiil=angularD6MRWT(x,[4 4],precon,'inverse'); fitil=angularD6MRWT(x,[4 4],precon,'inversetranspose'); ftil=angularD6MRWT(x,[4 4],precon,'transpose'); minmax(ff(:)-ffil(:)) minmax(fi(:)-fiil(:)) minmax(fit(:)-fitil(:)) minmax(ft(:)-ftil(:)) rmpath(whereitsat) addpath(whereiamat) end
You can’t perform that action at this time.