www.gusucode.com > matlab非局部均值工具箱 > matlab非局部均值工具箱/matlab非局部均值工具箱/toolbox_nlmeans/toolbox/perform_wavelet_transform.m
function y = perform_wavelet_transform(x, Jmin, dir, options) % perform_wavelet_transform - wrapper to wavelab Wavelet transform (1D/2D and orthogonal/biorthogonal). % % y = perform_wavelet_transform(x, Jmin, dir, options); % % 'x' is either a 1D or a 2D array. % 'Jmin' is the minimum scale (i.e. the coarse channel is of size 2^Jmin % in 1D). % 'dir' is +1 for fwd transform and -1 for bwd. % 'options.wavelet_vm' is the number of Vanishing moment (both for primal and dual). % 'options.wavelet_type' can be % 'daubechies', 'symmlet', 'battle', 'biorthogonal'. % % Typical use : % M = <load your image here>; % Jmin = 4; % options.wavelet_type = 'biorthogonal_swapped'; % options.wavelet_vm = 4; % MW = perform_wavelet_transform(M, Jmin, +1, options); % Mt = <perform some modification on MW> % M = perform_wavelet_transform(Mt, Jmin, -1, options); % % 'y' is an array of the same size as 'x'. This means that for the 2D % we are stuck to the wavelab coding style, i.e. the result % of each transform is an array organized using Mallat's ordering % (whereas Matlab official toolbox use a 1D ordering for the 2D transform). % % Here the transform automaticaly select symmetric boundary condition % if you use a symmetric filter. If your filter is not symmetric % (e.g. Dauechies filters) then as the output must have same length % as the input, the boundary condition are automatically set to periodic. % % You do not need Wavelab to use this function (the Wavelab .m file are % included in this script). However, for faster execution time, you % should install the mex file within the Wavelab distribution. % http://www-stat.stanford.edu/~wavelab/ % % Copyright (c) 2005 Gabriel Peyr? if nargin<3 dir = 1; end if nargin<2 Jmin = 3; end options.null = 0; if isfield(options, 'wavelet_type') wavelet_type = options.wavelet_type; else wavelet_type = 'daubechies'; end if isfield(options, 'wavelet_vm') VM = options.wavelet_vm; else VM = 4; end if isfield(options, 'ti') % for translation-invariant transform ti = options.ti; else ti = 0; end ndim = length(size(x)); if ndim==2 && ( size(x,2)==1 || size(x,1)==1 ) ndim=1; end % for color images if ndims(x)>2 y = x; for i=1:size(x,3) y(:,:,i) = perform_wavelet_transform(x(:,:,i), Jmin, dir, options); end return; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % WAVELAB %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % generate filters switch lower(wavelet_type) case 'daubechies' qmf = MakeONFilter('Daubechies',VM*2); % in Wavelab, 2nd argument is VM*2 for Daubechies... no comment ... case 'haar' qmf = MakeONFilter('Haar'); % in Wavelab, 2nd argument is VM*2 for Daubechies... no comment ... case 'symmlet' qmf = MakeONFilter('Symmlet',VM); case 'battle' qmf = MakeONFilter('Battle',VM-1); dqmf = qmf; % we need dual filter case 'biorthogonal' [qmf,dqmf] = MakeBSFilter( 'CDF', [VM,VM] ); case 'biorthogonal_swapped' [dqmf,qmf] = MakeBSFilter( 'CDF', [VM,VM] ); otherwise error('Unknown transform.'); end % translation invariant transform if ti if ndim==2 && size(x,2)<50 ndim = 1; end if ndim==1 if dir==1 y = TI2Stat( FWT_TI(x,Jmin,qmf) ); else y = IWT_TI( Stat2TI(x),qmf); end elseif ndim==2 if dir==1 y = FWT2_TI(x,Jmin,qmf); else y = IWT2_TI(x,Jmin,qmf); end end return; end % perform transform if ~exist('dqmf') %%% ORTHOGONAL %%% if ndim==1 if dir==1 y = FWT_PO(x,Jmin,qmf); else y = IWT_PO(x,Jmin,qmf); end elseif ndim==2 if dir==1 y = FWT2_PO(x,Jmin,qmf); else y = IWT2_PO(x,Jmin,qmf); end end else %%% BI-ORTHOGONAL %%% if ndim==1 if dir==1 y = FWT_SBS(x,Jmin,qmf,dqmf); else y = IWT_SBS(x,Jmin,qmf,dqmf); end elseif ndim==2 if dir==1 y = FWT2_SBS(x,Jmin,qmf,dqmf); else y = IWT2_SBS(x,Jmin,qmf,dqmf); end end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % WAVELAB DISTRIBUTION -- http://www-stat.stanford.edu/~wavelab/ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % WAVELAB DISTRIBUTION -- http://www-stat.stanford.edu/~wavelab/ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function f = MakeONFilter(Type,Par) % MakeONFilter -- Generate Orthonormal QMF Filter for Wavelet Transform % Usage % qmf = MakeONFilter(Type,Par) % Inputs % Type string, 'Haar', 'Beylkin', 'Coiflet', 'Daubechies', % 'Symmlet', 'Vaidyanathan','Battle' % Par integer, it is a parameter related to the support and vanishing % moments of the wavelets, explained below for each wavelet. % % Outputs % qmf quadrature mirror filter % % Description % The Haar filter (which could be considered a Daubechies-2) was the % first wavelet, though not called as such, and is discontinuous. % % The Beylkin filter places roots for the frequency response function % close to the Nyquist frequency on the real axis. % % The Coiflet filters are designed to give both the mother and father % wavelets 2*Par vanishing moments; here Par may be one of 1,2,3,4 or 5. % % The Daubechies filters are minimal phase filters that generate wavelets % which have a minimal support for a given number of vanishing moments. % They are indexed by their length, Par, which may be one of % 4,6,8,10,12,14,16,18 or 20. The number of vanishing moments is par/2. % % Symmlets are also wavelets within a minimum size support for a given % number of vanishing moments, but they are as symmetrical as possible, % as opposed to the Daubechies filters which are highly asymmetrical. % They are indexed by Par, which specifies the number of vanishing % moments and is equal to half the size of the support. It ranges % from 4 to 10. % % The Vaidyanathan filter gives an exact reconstruction, but does not % satisfy any moment condition. The filter has been optimized for % speech coding. % % The Battle-Lemarie filter generate spline orthogonal wavelet basis. % The parameter Par gives the degree of the spline. The number of % vanishing moments is Par+1. % % See Also % FWT_PO, IWT_PO, FWT2_PO, IWT2_PO, WPAnalysis % % References % The books by Daubechies and Wickerhauser. % if strcmp(Type,'Haar'), f = [1 1] ./ sqrt(2); end if strcmp(Type,'Beylkin'), f = [ .099305765374 .424215360813 .699825214057 ... .449718251149 -.110927598348 -.264497231446 ... .026900308804 .155538731877 -.017520746267 ... -.088543630623 .019679866044 .042916387274 ... -.017460408696 -.014365807969 .010040411845 ... .001484234782 -.002736031626 .000640485329 ]; end if strcmp(Type,'Coiflet'), if Par==1, f = [ .038580777748 -.126969125396 -.077161555496 ... .607491641386 .745687558934 .226584265197 ]; end if Par==2, f = [ .016387336463 -.041464936782 -.067372554722 ... .386110066823 .812723635450 .417005184424 ... -.076488599078 -.059434418646 .023680171947 ... .005611434819 -.001823208871 -.000720549445 ]; end if Par==3, f = [ -.003793512864 .007782596426 .023452696142 ... -.065771911281 -.061123390003 .405176902410 ... .793777222626 .428483476378 -.071799821619 ... -.082301927106 .034555027573 .015880544864 ... -.009007976137 -.002574517688 .001117518771 ... .000466216960 -.000070983303 -.000034599773 ]; end if Par==4, f = [ .000892313668 -.001629492013 -.007346166328 ... .016068943964 .026682300156 -.081266699680 ... -.056077313316 .415308407030 .782238930920 ... .434386056491 -.066627474263 -.096220442034 ... .039334427123 .025082261845 -.015211731527 ... -.005658286686 .003751436157 .001266561929 ... -.000589020757 -.000259974552 .000062339034 ... .000031229876 -.000003259680 -.000001784985 ]; end if Par==5, f = [ -.000212080863 .000358589677 .002178236305 ... -.004159358782 -.010131117538 .023408156762 ... .028168029062 -.091920010549 -.052043163216 ... .421566206729 .774289603740 .437991626228 ... -.062035963906 -.105574208706 .041289208741 ... .032683574283 -.019761779012 -.009164231153 ... .006764185419 .002433373209 -.001662863769 ... -.000638131296 .000302259520 .000140541149 ... -.000041340484 -.000021315014 .000003734597 ... .000002063806 -.000000167408 -.000000095158 ]; end end if strcmp(Type,'Daubechies'), if Par==4, f = [ .482962913145 .836516303738 ... .224143868042 -.129409522551 ]; end if Par==6, f = [ .332670552950 .806891509311 ... .459877502118 -.135011020010 ... -.085441273882 .035226291882 ]; end if Par==8, f = [ .230377813309 .714846570553 ... .630880767930 -.027983769417 ... -.187034811719 .030841381836 ... .032883011667 -.010597401785 ]; end if Par==10, f = [ .160102397974 .603829269797 .724308528438 ... .138428145901 -.242294887066 -.032244869585 ... .077571493840 -.006241490213 -.012580751999 ... .003335725285 ]; end if Par==12, f = [ .111540743350 .494623890398 .751133908021 ... .315250351709 -.226264693965 -.129766867567 ... .097501605587 .027522865530 -.031582039317 ... .000553842201 .004777257511 -.001077301085 ]; end if Par==14, f = [ .077852054085 .396539319482 .729132090846 ... .469782287405 -.143906003929 -.224036184994 ... .071309219267 .080612609151 -.038029936935 ... -.016574541631 .012550998556 .000429577973 ... -.001801640704 .000353713800 ]; end if Par==16, f = [ .054415842243 .312871590914 .675630736297 ... .585354683654 -.015829105256 -.284015542962 ... .000472484574 .128747426620 -.017369301002 ... -.044088253931 .013981027917 .008746094047 ... -.004870352993 -.000391740373 .000675449406 ... -.000117476784 ]; end if Par==18, f = [ .038077947364 .243834674613 .604823123690 ... .657288078051 .133197385825 -.293273783279 ... -.096840783223 .148540749338 .030725681479 ... -.067632829061 .000250947115 .022361662124 ... -.004723204758 -.004281503682 .001847646883 ... .000230385764 -.000251963189 .000039347320 ]; end if Par==20, f = [ .026670057901 .188176800078 .527201188932 ... .688459039454 .281172343661 -.249846424327 ... -.195946274377 .127369340336 .093057364604 ... -.071394147166 -.029457536822 .033212674059 ... .003606553567 -.010733175483 .001395351747 ... .001992405295 -.000685856695 -.000116466855 ... .000093588670 -.000013264203 ]; end end if strcmp(Type,'Symmlet'), if Par==4, f = [ -.107148901418 -.041910965125 .703739068656 ... 1.136658243408 .421234534204 -.140317624179 ... -.017824701442 .045570345896 ]; end if Par==5, f = [ .038654795955 .041746864422 -.055344186117 ... .281990696854 1.023052966894 .896581648380 ... .023478923136 -.247951362613 -.029842499869 ... .027632152958 ]; end if Par==6, f = [ .021784700327 .004936612372 -.166863215412 ... -.068323121587 .694457972958 1.113892783926 ... .477904371333 -.102724969862 -.029783751299 ... .063250562660 .002499922093 -.011031867509 ]; end if Par==7, f = [ .003792658534 -.001481225915 -.017870431651 ... .043155452582 .096014767936 -.070078291222 ... .024665659489 .758162601964 1.085782709814 ... .408183939725 -.198056706807 -.152463871896 ... .005671342686 .014521394762 ]; end if Par==8, f = [ .002672793393 -.000428394300 -.021145686528 ... .005386388754 .069490465911 -.038493521263 ... -.073462508761 .515398670374 1.099106630537 ... .680745347190 -.086653615406 -.202648655286 ... .010758611751 .044823623042 -.000766690896 ... -.004783458512 ]; end if Par==9, f = [ .001512487309 -.000669141509 -.014515578553 ... .012528896242 .087791251554 -.025786445930 ... -.270893783503 .049882830959 .873048407349 ... 1.015259790832 .337658923602 -.077172161097 ... .000825140929 .042744433602 -.016303351226 ... -.018769396836 .000876502539 .001981193736 ]; end if Par==10, f = [ .001089170447 .000135245020 -.012220642630 ... -.002072363923 .064950924579 .016418869426 ... -.225558972234 -.100240215031 .667071338154 ... 1.088251530500 .542813011213 -.050256540092 ... -.045240772218 .070703567550 .008152816799 ... -.028786231926 -.001137535314 .006495728375 ... .000080661204 -.000649589896 ]; end end if strcmp(Type,'Vaidyanathan'), f = [ -.000062906118 .000343631905 -.000453956620 ... -.000944897136 .002843834547 .000708137504 ... -.008839103409 .003153847056 .019687215010 ... -.014853448005 -.035470398607 .038742619293 ... .055892523691 -.077709750902 -.083928884366 ... .131971661417 .135084227129 -.194450471766 ... -.263494802488 .201612161775 .635601059872 ... .572797793211 .250184129505 .045799334111 ]; end if strcmp(Type,'Battle'), if Par == 1, g = [0.578163 0.280931 -0.0488618 -0.0367309 ... 0.012003 0.00706442 -0.00274588 -0.00155701 ... 0.000652922 0.000361781 -0.000158601 -0.0000867523 ]; end if Par == 3, g = [0.541736 0.30683 -0.035498 -0.0778079 ... 0.0226846 0.0297468 -0.0121455 -0.0127154 ... 0.00614143 0.00579932 -0.00307863 -0.00274529 ... 0.00154624 0.00133086 -0.000780468 -0.00065562 ... 0.000395946 0.000326749 -0.000201818 -0.000164264 ... 0.000103307 ]; end if Par == 5, g = [0.528374 0.312869 -0.0261771 -0.0914068 ... 0.0208414 0.0433544 -0.0148537 -0.0229951 ... 0.00990635 0.0128754 -0.00639886 -0.00746848 ... 0.00407882 0.00444002 -0.00258816 -0.00268646 ... 0.00164132 0.00164659 -0.00104207 -0.00101912 ... 0.000662836 0.000635563 -0.000422485 -0.000398759 ... 0.000269842 0.000251419 -0.000172685 -0.000159168 ... 0.000110709 0.000101113 ]; end l = length(g); f = zeros(1,2*l-1); f(l:2*l-1) = g; f(1:l-1) = reverse(g(2:l)); end f = f ./ norm(f); % % Copyright (c) 1993-5. Jonathan Buckheit and David Donoho % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function wcoef = FWT_PO(x,L,qmf) % FWT_PO -- Forward Wavelet Transform (periodized, orthogonal) % Usage % wc = FWT_PO(x,L,qmf) % Inputs % x 1-d signal; length(x) = 2^J % L Coarsest Level of V_0; L << J % qmf quadrature mirror filter (orthonormal) % Outputs % wc 1-d wavelet transform of x. % % Description % 1. qmf filter may be obtained from MakeONFilter % 2. usually, length(qmf) < 2^(L+1) % 3. To reconstruct use IWT_PO % % See Also % IWT_PO, MakeONFilter % [n,J] = dyadlength(x) ; wcoef = zeros(1,n) ; beta = ShapeAsRow(x); %take samples at finest scale as beta-coeffts for j=J-1:-1:L alfa = DownDyadHi(beta,qmf); wcoef(dyad(j)) = alfa; beta = DownDyadLo(beta,qmf) ; end wcoef(1:(2^L)) = beta; wcoef = ShapeLike(wcoef,x); % % Copyright (c) 1993. Iain M. Johnstone % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function [n,J] = dyadlength(x) % dyadlength -- Find length and dyadic length of array % Usage % [n,J] = dyadlength(x) % Inputs % x array of length n = 2^J (hopefully) % Outputs % n length(x) % J least power of two greater than n % % Side Effects % A warning is issued if n is not a power of 2. % % See Also % quadlength, dyad, dyad2ix % n = length(x) ; J = ceil(log(n)/log(2)); if 2^J ~= n , disp('Warning in dyadlength: n != 2^J') end % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function row = ShapeAsRow(sig) % ShapeAsRow -- Make signal a row vector % Usage % row = ShapeAsRow(sig) % Inputs % sig a row or column vector % Outputs % row a row vector % % See Also % ShapeLike % row = sig(:)'; % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function d = DownDyadHi(x,qmf) % DownDyadHi -- Hi-Pass Downsampling operator (periodized) % Usage % d = DownDyadHi(x,f) % Inputs % x 1-d signal at fine scale % f filter % Outputs % y 1-d signal at coarse scale % % See Also % DownDyadLo, UpDyadHi, UpDyadLo, FWT_PO, iconv % d = iconv( MirrorFilt(qmf),lshift(x)); n = length(d); d = d(1:2:(n-1)); % % Copyright (c) 1993. Iain M. Johnstone % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function y = MirrorFilt(x) % MirrorFilt -- Apply (-1)^t modulation % Usage % h = MirrorFilt(l) % Inputs % l 1-d signal % Outputs % h 1-d signal with DC frequency content shifted % to Nyquist frequency % % Description % h(t) = (-1)^(t-1) * x(t), 1 <= t <= length(x) % % See Also % DyadDownHi % y = -( (-1).^(1:length(x)) ).*x; % % Copyright (c) 1993. Iain M. Johnstone % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function y = lshift(x) % lshift -- Circular left shift of 1-d signal % Usage % l = lshift(x) % Inputs % x 1-d signal % Outputs % l 1-d signal % l(i) = x(i+1) except l(n) = x(1) % y = [ x( 2:length(x) ) x(1) ]; % % Copyright (c) 1993. Iain M. Johnstone % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function y = iconv(f,x) % iconv -- Convolution Tool for Two-Scale Transform % Usage % y = iconv(f,x) % Inputs % f filter % x 1-d signal % Outputs % y filtered result % % Description % Filtering by periodic convolution of x with f % % See Also % aconv, UpDyadHi, UpDyadLo, DownDyadHi, DownDyadLo % n = length(x); p = length(f); if p <= n, xpadded = [x((n+1-p):n) x]; else z = zeros(1,p); for i=1:p, imod = 1 + rem(p*n -p + i-1,n); z(i) = x(imod); end xpadded = [z x]; end ypadded = filter(f,1,xpadded); y = ypadded((p+1):(n+p)); % % Copyright (c) 1993. David L. Donoho % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function i = dyad(j) % dyad -- Index entire j-th dyad of 1-d wavelet xform % Usage % ix = dyad(j); % Inputs % j integer % Outputs % ix list of all indices of wavelet coeffts at j-th level % i = (2^(j)+1):(2^(j+1)) ; % % Copyright (c) 1993. David L. Donoho % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function d = DownDyadLo(x,qmf) % DownDyadLo -- Lo-Pass Downsampling operator (periodized) % Usage % d = DownDyadLo(x,f) % Inputs % x 1-d signal at fine scale % f filter % Outputs % y 1-d signal at coarse scale % % See Also % DownDyadHi, UpDyadHi, UpDyadLo, FWT_PO, aconv % d = aconv(qmf,x); n = length(d); d = d(1:2:(n-1)); % % Copyright (c) 1993. Iain M. Johnstone % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function y = aconv(f,x) % aconv -- Convolution Tool for Two-Scale Transform % Usage % y = aconv(f,x) % Inputs % f filter % x 1-d signal % Outputs % y filtered result % % Description % Filtering by periodic convolution of x with the % time-reverse of f. % % See Also % iconv, UpDyadHi, UpDyadLo, DownDyadHi, DownDyadLo % n = length(x); p = length(f); if p < n, xpadded = [x x(1:p)]; else z = zeros(1,p); for i=1:p, imod = 1 + rem(i-1,n); z(i) = x(imod); end xpadded = [x z]; end fflip = reverse(f); ypadded = filter(fflip,1,xpadded); y = ypadded(p:(n+p-1)); % % Copyright (c) 1993. David L. Donoho % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function vec = ShapeLike(sig,proto) % ShapeLike -- Make 1-d signal with given shape % Usage % vec = ShapeLike(sig,proto) % Inputs % sig a row or column vector % proto a prototype shape (row or column vector) % Outputs % vec a vector with contents taken from sig % and same shape as proto % % See Also % ShapeAsRow % sp = size(proto); ss = size(sig); if( sp(1)>1 & sp(2)>1 ) disp('Weird proto argument to ShapeLike') elseif ss(1)>1 & ss(2) > 1, disp('Weird sig argument to ShapeLike') else if(sp(1) > 1), if ss(1) > 1, vec = sig; else vec = sig(:); end else if ss(2) > 1, vec = sig; else vec = sig(:)'; end end end % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function x = IWT_PO(wc,L,qmf) % IWT_PO -- Inverse Wavelet Transform (periodized, orthogonal) % Usage % x = IWT_PO(wc,L,qmf) % Inputs % wc 1-d wavelet transform: length(wc) = 2^J. % L Coarsest scale (2^(-L) = scale of V_0); L << J; % qmf quadrature mirror filter % Outputs % x 1-d signal reconstructed from wc % % Description % Suppose wc = FWT_PO(x,L,qmf) where qmf is an orthonormal quad. mirror % filter, e.g. one made by MakeONFilter. Then x can be reconstructed by % x = IWT_PO(wc,L,qmf) % % See Also % FWT_PO, MakeONFilter % wcoef = ShapeAsRow(wc); x = wcoef(1:2^L); [n,J] = dyadlength(wcoef); for j=L:J-1 x = UpDyadLo(x,qmf) + UpDyadHi(wcoef(dyad(j)),qmf) ; end x = ShapeLike(x,wc); % % Copyright (c) 1993. Iain M. Johnstone % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function y = UpDyadLo(x,qmf) % UpDyadLo -- Lo-Pass Upsampling operator; periodized % Usage % u = UpDyadLo(d,f) % Inputs % d 1-d signal at coarser scale % f filter % Outputs % u 1-d signal at finer scale % % See Also % DownDyadLo, DownDyadHi, UpDyadHi, IWT_PO, iconv % y = iconv(qmf, UpSample(x) ); % % Copyright (c) 1993. Iain M. Johnstone % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function y = UpSample(x,s) % UpSample -- Upsampling operator % Usage % u = UpSample(d[,s]) % Inputs % d 1-d signal, of length n % s upsampling scale, default = 2 % Outputs % u 1-d signal, of length s*n with zeros % interpolating alternate samples % u(s*i-1) = d(i), i=1,...,n % if nargin == 1, s = 2; end n = length(x)*s; y = zeros(1,n); y(1:s:(n-s+1) )=x; % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function y = UpDyadHi(x,qmf) % UpDyadHi -- Hi-Pass Upsampling operator; periodized % Usage % u = UpDyadHi(d,f) % Inputs % d 1-d signal at coarser scale % f filter % Outputs % u 1-d signal at finer scale % % See Also % DownDyadLo, DownDyadHi, UpDyadLo, IWT_PO, aconv % y = aconv( MirrorFilt(qmf), rshift( UpSample(x) ) ); % % Copyright (c) 1993. Iain M. Johnstone % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function y = rshift(x) % rshift -- Circular right shift of 1-d signal % Usage % r = rshift(x) % Inputs % x 1-d signal % Outputs % r 1-d signal % r(i) = x(i-1) except r(1) = x(n) % n = length(x); y = [ x(n) x( 1: (n-1) )]; % % Copyright (c) 1993. Iain M. Johnstone % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function [qmf,dqmf] = MakeBSFilter(Type,Par) % MakeBSFilter -- Generate Biorthonormal QMF Filter Pairs % Usage % [qmf,dqmf] = MakeBSFilter(Type,Par) % Inputs % Type string, one of: % 'Triangle' % 'Interpolating' 'Deslauriers' (the two are same) % 'Average-Interpolating' % 'CDF' (spline biorthogonal filters in Daubechies's book) % 'Villasenor' (Villasenor's 5 best biorthogonal filters) % Par integer list, e.g. if Type ='Deslauriers', Par=3 specifies % Deslauriers-Dubuc filter, polynomial degree 3 % Outputs % qmf quadrature mirror filter (odd length, symmetric) % dqmf dual quadrature mirror filter (odd length, symmetric) % % See Also % FWT_PBS, IWT_PBS, FWT2_PBS, IWT2_PBS % % References % I. Daubechies, "Ten Lectures on Wavelets." % % G. Deslauriers and S. Dubuc, "Symmetric Iterative Interpolating Processes." % % D. Donoho, "Smooth Wavelet Decompositions with Blocky Coefficient Kernels." % % J. Villasenor, B. Belzer and J. Liao, "Wavelet Filter Evaluation for % Image Compression." % if nargin < 2, Par = 0; end sqr2 = sqrt(2); if strcmp(Type,'Triangle'), qmf = [0 1 0]; dqmf = [.5 1 .5]; elseif strcmp(Type,'Interpolating') | strcmp(Type,'Deslauriers'), qmf = [0 1 0]; dqmf = MakeDDFilter(Par)'; dqmf = dqmf(1:(length(dqmf)-1)); elseif strcmp(Type,'Average-Interpolating'), qmf = [0 .5 .5] ; dqmf = [0 ; MakeAIFilter(Par)]'; elseif strcmp(Type,'CDF'), if Par(1)==1, dqmf = [0 .5 .5] .* sqr2; if Par(2) == 1, qmf = [0 .5 .5] .* sqr2; elseif Par(2) == 3, qmf = [0 -1 1 8 8 1 -1] .* sqr2 / 16; elseif Par(2) == 5, qmf = [0 3 -3 -22 22 128 128 22 -22 -3 3].*sqr2/256; end elseif Par(1)==2, dqmf = [.25 .5 .25] .* sqr2; if Par(2)==2, qmf = [-.125 .25 .75 .25 -.125] .* sqr2; elseif Par(2)==4, qmf = [3 -6 -16 38 90 38 -16 -6 3] .* (sqr2/128); elseif Par(2)==6, qmf = [-5 10 34 -78 -123 324 700 324 -123 -78 34 10 -5 ] .* (sqr2/1024); elseif Par(2)==8, qmf = [35 -70 -300 670 1228 -3126 -3796 10718 22050 ... 10718 -3796 -3126 1228 670 -300 -70 35 ] .* (sqr2/32768); end elseif Par(1)==3, dqmf = [0 .125 .375 .375 .125] .* sqr2; if Par(2) == 1, qmf = [0 -.25 .75 .75 -.25] .* sqr2; elseif Par(2) == 3, qmf = [0 3 -9 -7 45 45 -7 -9 3] .* sqr2/64; elseif Par(2) == 5, qmf = [0 -5 15 19 -97 -26 350 350 -26 -97 19 15 -5] .* sqr2/512; elseif Par(2) == 7, qmf = [0 35 -105 -195 865 363 -3489 -307 11025 11025 -307 -3489 363 865 -195 -105 35] .* sqr2/16384; elseif Par(2) == 9, qmf = [0 -63 189 469 -1911 -1308 9188 1140 -29676 190 87318 87318 190 -29676 ... 1140 9188 -1308 -1911 469 189 -63] .* sqr2/131072; end elseif Par(1)==4, dqmf = [.026748757411 -.016864118443 -.078223266529 .266864118443 .602949018236 ... .266864118443 -.078223266529 -.016864118443 .026748757411] .*sqr2; if Par(2) == 4, qmf = [0 -.045635881557 -.028771763114 .295635881557 .557543526229 ... .295635881557 -.028771763114 -.045635881557 0] .*sqr2; end end elseif strcmp(Type,'Villasenor'), if Par == 1, % The "7-9 filters" qmf = [.037828455506995 -.023849465019380 -.11062440441842 .37740285561265]; qmf = [qmf .85269867900940 reverse(qmf)]; dqmf = [-.064538882628938 -.040689417609558 .41809227322221]; dqmf = [dqmf .78848561640566 reverse(dqmf)]; elseif Par == 2, qmf = [-.008473 .003759 .047282 -.033475 -.068867 .383269 .767245 .383269 -.068867... -.033475 .047282 .003759 -.008473]; dqmf = [0.014182 0.006292 -0.108737 -0.069163 0.448109 .832848 .448109 -.069163 -.108737 .006292 .014182]; elseif Par == 3, qmf = [0 -.129078 .047699 .788486 .788486 .047699 -.129078]; dqmf = [0 .018914 .006989 -.067237 .133389 .615051 .615051 .133389 -.067237 .006989 .018914]; elseif Par == 4, qmf = [-1 2 6 2 -1] / (4*sqr2); dqmf = [1 2 1] / (2*sqr2); elseif Par == 5, qmf = [0 1 1]/sqr2; dqmf = [0 -1 1 8 8 1 -1]/(8*sqr2); end end % % Copyright (c) 1995. Jonathan Buckheit, Shaobing Chen and David Donoho % % Modified by Maureen Clerc and Jerome Kalifa, 1997 % clerc@cmapx.polytechnique.fr, kalifa@cmapx.polytechnique.fr % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function wcoef = FWT_SBS(x,L,qmf,dqmf) % FWT_SBS -- Forward Wavelet Transform (symmetric extension, biorthogonal, symmetric) % Usage % wc = FWT_SBS(x,L,qmf,dqmf) % Inputs % x 1-d signal; arbitrary length % L Coarsest Level of V_0; L << J % qmf quadrature mirror filter (symmetric) % dqmf quadrature mirror filter (symmetric, dual of qmf) % Outputs % wc 1-d wavelet transform of x. % % Description % 1. qmf filter may be obtained from MakePBSFilter % 2. usually, length(qmf) < 2^(L+1) % 3. To reconstruct use IWT_SBS % % See Also % IWT_SBS, MakePBSFilter % % References % Based on the algorithm developed by Christopher Brislawn. % See "Classification of Symmetric Wavelet Transforms" % [n,J] = dyadlength(x); wcoef = zeros(1,n); beta = ShapeAsRow(x); % take samples at finest scale as beta-coeffts dp = dyadpartition(n); for j=J-1:-1:L, [beta, alfa] = DownDyad_SBS(beta,qmf,dqmf); dyadj = (dp(j+1)+1):dp(j+2); wcoef(dyadj) = alfa; end wcoef(1:length(beta)) = beta; wcoef = ShapeLike(wcoef,x); % % Copyright (c) 1996. Thomas P.Y. Yu % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function dp = dyadpartition(n) % dyadpartition -- determine dyadic partition in wavelet transform of % nondyadic signals J = ceil(log2(n)); m = n; for j=J-1:-1:0; if rem(m,2)==0, dps(j+1) = m/2; m = m/2; else dps(j+1) = (m-1)/2; m = (m+1)/2; end end dp = cumsum([1 dps]); % % Copyright (c) 1996. Thomas P.Y. Yu % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function [beta, alpha] = DownDyad_SBS(x,qmf,dqmf) % DownDyad_SBS -- Symmetric Downsampling operator % Usage % [beta,alpha] = DownDyad_SBS(x,qmf,dqmf) % Inputs % x 1-d signal at fine scale % qmf quadrature mirror filter % dqmf dual quadrature mirror filter % Outputs % beta coarse coefficients % alpha fine coefficients % See Also % FWT_SBS % % oddf = (rem(length(qmf),2)==1); oddf = ~(qmf(1)==0 & qmf(length(qmf))~=0); oddx = (rem(length(x),2)==1); % symmetric extension of x if oddf, ex = extend(x,1,1); else ex = extend(x,2,2); end % convolution ebeta = DownDyadLo_PBS(ex,qmf); ealpha = DownDyadHi_PBS(ex,dqmf); % project if oddx, beta = ebeta(1:(length(x)+1)/2); alpha = ealpha(1:(length(x)-1)/2); else beta = ebeta(1:length(x)/2); alpha = ealpha(1:length(x)/2); end % % Copyright (c) 1996. Thomas P.Y. Yu % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function y = extend(x, par1, par2) % extend -- perform various kinds of symmetric extension % if par1==1 & par2==1, y = [x x((length(x)-1):-1:2)]; elseif par1==1 & par2==2, y = [x x((length(x)-1):-1:1)]; elseif par1==2 & par2==1, y = [x x(length(x):-1:2)]; elseif par1==2 & par2==2, y = [x reverse(x)]; end % % Copyright (c) 1996. Thomas P.Y. Yu % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function d = DownDyadLo_PBS(x,qmf) % DownDyadLo_PBS -- Lo-Pass Downsampling operator (periodized,symmetric) % Usage % d = DownDyadLo_PBS(x,sf) % Inputs % x 1-d signal at fine scale % sf symmetric filter % Outputs % y 1-d signal at coarse scale % % See Also % DownDyadHi_PBS, UpDyadHi_PBS, UpDyadLo_PBS, FWT_PBSi, symm_aconv % d = symm_aconv(qmf,x); n = length(d); d = d(1:2:(n-1)); % % Copyright (c) 1995. David L. Donoho % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function y = symm_aconv(sf,x) % symm_aconv -- Symmetric Convolution Tool for Two-Scale Transform % Usage % y = symm_aconv(sf,x) % Inputs % sf symmetric filter % x 1-d signal % Outputs % y filtered result % % Description % Filtering by periodic convolution of x with the % time-reverse of sf. % % See Also % symm_iconv, UpDyadHi_PBS, UpDyadLo_PBS, DownDyadHi_PBS, DownDyadLo_PBS % n = length(x); p = length(sf); if p < n, xpadded = [x x(1:p)]; else z = zeros(1,p); for i=1:p, imod = 1 + rem(i-1,n); z(i) = x(imod); end xpadded = [x z]; end fflip = reverse(sf); ypadded = filter(fflip,1,xpadded); if p < n, y = [ypadded((n+1):(n+p)) ypadded((p+1):(n))]; else for i=1:n, imod = 1 + rem(p+i-1,n); y(imod) = ypadded(p+i); end end shift = (p-1)/ 2 ; shift = 1 + rem(shift-1, n); y = [y((1+shift):n) y(1:(shift))] ; % % Copyright (c) 1995. Shaobing Chen and David L. Donoho % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function d = DownDyadHi_PBS(x,qmf) % DownDyadHi_PBS -- Hi-Pass Downsampling operator (periodized,symmetric) % Usage % d = DownDyadHi_PBS(x,sqmf) % Inputs % x 1-d signal at fine scale % sqmf symmetric filter % Outputs % y 1-d signal at coarse scale % % See Also % DownDyadLo_PBS, UpDyadHi_PBS, UpDyadLo_PBS, FWT_PBS, symm_iconv % d = symm_iconv( MirrorSymmFilt(qmf),lshift(x)); n = length(d); d = d(1:2:(n-1)); % % Copyright (c) 1995. David L. Donoho % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function y = MirrorSymmFilt(x) % MirrorSymmFilt -- apply (-1)^t modulation to symmetric filter % Usage % h = MirrorSymmFilt(l) % Inputs % l symmetric filter % Outputs % h symmetric filter with DC frequency content shifted % to Nyquist frequency % % Description % h(t) = (-1)^t * x(t), -k <= t <= k ; length(x)=2k+1 % % See Also % DownDyadHi_PBS % k = (length(x)-1)/2; y = ( (-1).^((-k):k) ) .*x; % % Copyright (c) 1993. Iain M. Johnstone % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function y = symm_iconv(sf,x) % symm_iconv -- Symmetric Convolution Tool for Two-Scale Transform % Usage % y = iconv(sf,x) % Inputs % sf symmetric filter % x 1-d signal % Output % y filtered result % % Description % Filtering by periodic convolution of x with sf % % See Also % symm_aconv, UpDyadHi_PBS, UpDyadLo_PBS, DownDyadHi_PBS, DownDyadLo_PBS % n = length(x); p = length(sf); if p <= n, xpadded = [x((n+1-p):n) x]; else z = zeros(1,p); for i=1:p, imod = 1 + rem(p*n -p + i-1,n); z(i) = x(imod); end xpadded = [z x]; end ypadded = filter(sf,1,xpadded); y = ypadded((p+1):(n+p)); shift = (p+1)/2; shift = 1 + rem(shift-1, n); y = [y(shift:n) y(1:(shift-1))]; % % Copyright (c) 1995. Shaobing Chen and David L. Donoho % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function x = IWT_SBS(wc,L,qmf,dqmf) % iwt_po -- Inverse Wavelet Transform (symmetric extension, biorthogonal, symmetric) % Usage % x = IWT_SBS(wc,L,qmf,dqmf) % Inputs % wc 1-d wavelet transform: length(wc)= 2^J. % L Coarsest scale (2^(-L) = scale of V_0); L << J; % qmf quadrature mirror filter % dqmf dual quadrature mirror filter (symmetric, dual of qmf) % Outputs % x 1-d signal reconstructed from wc % Description % Suppose wc = FWT_SBS(x,L,qmf,dqmf) where qmf and dqmf are orthonormal % quad. mirror filters made by MakeBioFilter. Then x can be reconstructed % by % x = IWT_SBS(wc,L,qmf,dqmf) % See Also: % FWT_SBS, MakeBioFilter % wcoef = ShapeAsRow(wc); [n,J] = dyadlength(wcoef); dp = dyadpartition(n); x = wcoef(1:dp(L+1)); for j=L:J-1, dyadj = (dp(j+1)+1):dp(j+2); x = UpDyad_SBS(x, wcoef(dyadj), qmf, dqmf); end x = ShapeLike(x,wc); % % Copyright (c) 1996. Thomas P.Y. Yu % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function x = UpDyad_SBS(beta,alpha,qmf,dqmf) % UpDyad_SBS -- Symmetric Upsampling operator % Usage % x = UpDyad_SBS(beta,alpha,qmf,dqmf) % Inputs % beta coarse coefficients % alpha fine coefficients % qmf quadrature mirror filter % dqmf dual quadrature mirror filter % Outputs % x 1-d signal at fine scale % See Also % DownDyad_SBS, IWT_SBS % % oddf = (rem(length(qmf),2)==1); oddf = ~(qmf(1)==0 & qmf(length(qmf))~=0); oddx = (length(beta) ~= length(alpha)); L = length(beta)+length(alpha); if oddf, if oddx, ebeta = extend(beta,1,1); ealpha = extend(alpha,2,2); else ebeta = extend(beta,2,1); ealpha = extend(alpha,1,2); end else if oddx, ebeta = extend(beta,1,2); ealpha = [alpha 0 -reverse(alpha)]; else ebeta = extend(beta,2,2); ealpha = [alpha -reverse(alpha)]; end end coarse = UpDyadLo_PBS(ebeta,dqmf); fine = UpDyadHi_PBS(ealpha,qmf); x = coarse + fine; x = x(1:L); % % Copyright (c) 1996. Thomas P.Y. Yu % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function y = UpDyadLo_PBS(x,qmf) % UpDyadLo_PBS -- Lo-Pass Upsampling operator; periodized % Usage % u = UpDyadLo_PBS(d,sf) % Inputs % d 1-d signal at coarser scale % sf symmetric filter % Outputs % u 1-d signal at finer scale % % See Also % DownDyadLo_PBS , DownDyadHi_PBS , UpDyadHi_PBS, IWT_PBS, symm_iconv % y = symm_iconv(qmf, UpSample(x,2) ); % % Copyright (c) 1995. David L. Donoho % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function y = UpDyadHi_PBS(x,qmf) % UpDyadHi_PBS -- Hi-Pass Upsampling operator; periodized % Usage % u = UpDyadHi_PBS(d,f) % Inputs % d 1-d signal at coarser scale % sf symmetric filter % Outputs % u 1-d signal at finer scale % % See Also % DownDyadLo_PBS, DownDyadHi_PBS, UpDyadLo_PBS, IWT_PBS, symm_aconv % y = symm_aconv( MirrorSymmFilt(qmf), rshift( UpSample(x,2) ) ); % % Copyright (c) 1995. David L. Donoho % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function wc = FWT2_SBS(x,L,qmf,dqmf) % FWT2_SBS -- 2-dimensional wavelet transform % (symmetric extension, bi-orthogonal) % Usage % wc = FWT2_SBS(x,L,qmf,dqmf) % Inputs % x 2-d image (n by n array, n arbitrary) % L coarsest level % qmf low-pass quadrature mirror filter % dqmf high-pass dual quadrature mirror filter % Output % wc 2-d wavelet transform % Description % A two-dimensional Wavelet Transform is computed for the % matrix x. To reconstruct, use IWT2_SBS. % See Also % IWT2_SBS % [m,J] = dyadlength(x(:,1)); [n,K] = dyadlength(x(1,:)); wc = x; mc = m; nc = n; J = min([J,K]); for jscal=J-1:-1:L, if rem(mc,2)==0, top = (mc/2+1):mc; bot = 1:(mc/2); else top = ((mc+1)/2+1):mc; bot = 1:((mc+1)/2); end if rem(nc,2)==0, right = (nc/2+1):nc; left = 1:(nc/2); else right = ((nc+1)/2+1):nc; left = 1:((nc+1)/2); end for ix=1:mc, row = wc(ix,1:nc); [beta,alpha] = DownDyad_SBS(row,qmf,dqmf); wc(ix,left) = beta; wc(ix,right) = alpha; end for iy=1:nc, column = wc(1:mc,iy)'; [beta,alpha] = DownDyad_SBS(column,qmf,dqmf); wc(bot,iy) = beta'; wc(top,iy) = alpha'; end mc = bot(length(bot)); nc = left(length(left)); end % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function x = IWT2_SBS(wc,L,qmf,dqmf) % IWT2_SBS -- Inverse 2d Wavelet Transform % (symmetric extention, bi-orthogonal) % Usage % x = IWT2_SBS(wc,L,qmf,dqmf) % Inputs % wc 2-d wavelet transform [n by n array, n arbitrary] % L coarse level % qmf low-pass quadrature mirror filter % dqmf high-pas dual quadrature mirror filter % Outputs % x 2-d signal reconstructed from wc % Description % If wc is the result of a forward 2d wavelet transform, with % wc = FWT2_SBS(x,L,qmf,dqmf) % then x = IWT2_SBS(wc,L,qmf,dqmf) reconstructs x exactly if qmf is a nice % quadrature mirror filter, e.g. one made by MakeBioFilter % See Also: % FWT2_SBS, MakeBioFilter % [m,J] = dyadlength(wc(:,1)); [n,K] = dyadlength(wc(1,:)); % assume m==n, J==K x = wc; dpm = dyadpartition(m); for jscal=L:J-1, bot = 1:dpm(jscal+1); top = (dpm(jscal+1)+1):dpm(jscal+2); all = [bot top]; nc = length(all); for iy=1:nc, x(all,iy) = UpDyad_SBS(x(bot,iy)', x(top,iy)', qmf, dqmf)'; end for ix=1:nc, x(ix,all) = UpDyad_SBS(x(ix,bot), x(ix,top), qmf, dqmf); end end % % Copyright (c) 1996. Thomas P.Y. Yu % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function wc = FWT2_PO(x,L,qmf) % FWT2_PO -- 2-d MRA wavelet transform (periodized, orthogonal) % Usage % wc = FWT2_PO(x,L,qmf) % Inputs % x 2-d image (n by n array, n dyadic) % L coarse level % qmf quadrature mirror filter % Outputs % wc 2-d wavelet transform % % Description % A two-dimensional Wavelet Transform is computed for the % array x. To reconstruct, use IWT2_PO. % % See Also % IWT2_PO, MakeONFilter % [n,J] = quadlength(x); wc = x; nc = n; for jscal=J-1:-1:L, top = (nc/2+1):nc; bot = 1:(nc/2); for ix=1:nc, row = wc(ix,1:nc); wc(ix,bot) = DownDyadLo(row,qmf); wc(ix,top) = DownDyadHi(row,qmf); end for iy=1:nc, row = wc(1:nc,iy)'; wc(top,iy) = DownDyadHi(row,qmf)'; wc(bot,iy) = DownDyadLo(row,qmf)'; end nc = nc/2; end % % Copyright (c) 1993. David L. Donoho % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function x = IWT2_PO(wc,L,qmf) % IWT2_PO -- Inverse 2-d MRA wavelet transform (periodized, orthogonal) % Usage % x = IWT2_PO(wc,L,qmf) % Inputs % wc 2-d wavelet transform [n by n array, n dyadic] % L coarse level % qmf quadrature mirror filter % Outputs % x 2-d signal reconstructed from wc % % Description % If wc is the result of a forward 2d wavelet transform, with % wc = FWT2_PO(x,L,qmf), then x = IWT2_PO(wc,L,qmf) reconstructs x % exactly if qmf is a nice qmf, e.g. one made by MakeONFilter. % % See Also % FWT2_PO, MakeONFilter % [n,J] = quadlength(wc); x = wc; nc = 2^(L+1); for jscal=L:J-1, top = (nc/2+1):nc; bot = 1:(nc/2); all = 1:nc; for iy=1:nc, x(all,iy) = UpDyadLo(x(bot,iy)',qmf)' ... + UpDyadHi(x(top,iy)',qmf)'; end for ix=1:nc, x(ix,all) = UpDyadLo(x(ix,bot),qmf) ... + UpDyadHi(x(ix,top),qmf); end nc = 2*nc; end % % Copyright (c) 1993. David L. Donoho % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function [n,J] = quadlength(x) % quadlength -- Find length and dyadic length of square matrix % Usage % [n,J] = quadlength(x) % Inputs % x 2-d image; size(n,n), n = 2^J (hopefully) % Outputs % n length(x) % J least power of two greater than n % % Side Effects % A warning message is issue if n is not a power of 2, % or if x is not a square matrix. % s = size(x); n = s(1); if s(2) ~= s(1), disp('Warning in quadlength: nr != nc') end k = 1 ; J = 0; while k < n , k=2*k; J = 1+J ; end ; if k ~= n , disp('Warning in quadlength: n != 2^J') end % % Copyright (c) 1993. David L. Donoho % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function wp = FWT_TI(x,L,qmf) % FWT_TI -- translation invariant forward wavelet transform % Usage % TIWT = FWT_TI(x,L,qmf) % Inputs % x array of dyadic length n=2^J % L degree of coarsest scale % qmf orthonormal quadrature mirror filter % Outputs % TIWT stationary wavelet transform table % formally same data structure as packet table % % See Also % IWT_TI % [n,J] = dyadlength(x); D = J-L; wp = zeros(n,D+1); x = ShapeAsRow(x); % wp(:,1) = x'; for d=0:(D-1), for b=0:(2^d-1), s = wp(packet(d,b,n),1)'; hsr = DownDyadHi(s,qmf); hsl = DownDyadHi(rshift(s),qmf); lsr = DownDyadLo(s,qmf); lsl = DownDyadLo(rshift(s),qmf); wp(packet(d+1,2*b ,n),d+2) = hsr'; wp(packet(d+1,2*b+1,n),d+2) = hsl'; wp(packet(d+1,2*b ,n),1 ) = lsr'; wp(packet(d+1,2*b+1,n),1 ) = lsl'; end end % % Copyright (c) 1994. David L. Donoho % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function tiwt = FWT2_TI(x,L,qmf) % FWT_TI -- 2-D translation invariant forward wavelet transform % Usage % TIWT = FWT2_TI(x,L,qmf) % Inputs % x 2-d image (n by n real array, n dyadic) % L degree of coarsest scale % qmf orthonormal quadrature mirror filter % Outputs % TIWT translation-invariant wavelet transform table, (3*(J-L)+1)*n by n % % See Also % IWT2_TI, IWT2_TIMedian % [n,J] = quadlength(x); D = J-L; tiwt = zeros((3*D+1)*n, n); lastx = (3*D*n+1):(3*D*n+n); lasty = 1:n; tiwt(lastx,lasty) = x; % for d=0:(D-1), l = J-d-1; ns = 2^(J-d); for b1=0:(2^d-1), for b2=0:(2^d-1), s = tiwt(3*D*n+packet(d,b1,n), packet(d,b2,n)); wc00 = FWT2_PO(s,l,qmf); wc01 = FWT2_PO(CircularShift(s,0,1),l,qmf); wc10 = FWT2_PO(CircularShift(s,1,0),l,qmf); wc11 = FWT2_PO(CircularShift(s,1,1),l,qmf); index10 = packet(d+1,2*b1,n); index20 = packet(d+1,2*b2,n); index11 = packet(d+1,2*b1+1,n); index21 = packet(d+1,2*b2+1,n); % horizontal stuff tiwt(3*d*n + index10 , index20) = wc00(1:(ns/2),(ns/2+1):ns); tiwt(3*d*n + index11, index20) = wc01(1:(ns/2),(ns/2+1):ns); tiwt(3*d*n + index10 , index21) = wc10(1:(ns/2),(ns/2+1):ns); tiwt(3*d*n + index11 , index21) = wc11(1:(ns/2),(ns/2+1):ns); % vertical stuff tiwt((3*d+1)*n + index10 , index20) = wc00((ns/2+1):ns,1:(ns/2)); tiwt((3*d+1)*n + index11, index20) = wc01((ns/2+1):ns,1:(ns/2)); tiwt((3*d+1)*n + index10 , index21) = wc10((ns/2+1):ns,1:(ns/2)); tiwt((3*d+1)*n + index11 , index21) = wc11((ns/2+1):ns,1:(ns/2)); % diagonal stuff tiwt((3*d+2)*n + index10 , index20) = wc00((ns/2+1):ns,(ns/2+1):ns); tiwt((3*d+2)*n + index11, index20) = wc01((ns/2+1):ns,(ns/2+1):ns); tiwt((3*d+2)*n + index10 , index21) = wc10((ns/2+1):ns,(ns/2+1):ns); tiwt((3*d+2)*n + index11 , index21) = wc11((ns/2+1):ns,(ns/2+1):ns); % low freq stuff tiwt(3*D*n + index10 , index20) = wc00(1:(ns/2),1:(ns/2)); tiwt(3*D*n + index11, index20) = wc01(1:(ns/2),1:(ns/2)); tiwt(3*D*n + index10 , index21) = wc10(1:(ns/2),1:(ns/2)); tiwt(3*D*n + index11 , index21) = wc11(1:(ns/2),1:(ns/2)); end, end end % % Copyright (c) 1995. David L. Donoho and Thomas P.Y. Yu % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function p=packet(d,b,n) % packet -- Packet table indexing % Usage % p = packet(d,b,n) % Inputs % d depth of splitting in packet decomposition % b block index among 2^d possibilities at depth d % n length of signal % Outputs % p linear indices of all coeff's in that block % npack = 2^d; p = ( (b * (n/npack) + 1) : ((b+1)*n/npack ) ) ; % % Copyright (c) 1993. David L. Donoho % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function x = IWT_TI(pkt,qmf) % IWT_TI -- Invert translation invariant wavelet transform % Usage % x = IWT_TI(TIWT,qmf) % Inputs % TIWT translation-invariant wavelet transform table % qmf quadrature mirror filter % Outputs % x 1-d signal reconstructed from translation-invariant % transform TIWT % % See Also % FWT_TI % [n,D1] = size(pkt); D = D1-1; J = log2(n); L = J-D; % wp = pkt; % sig = wp(:,1)'; for d= D-1:-1:0, for b=0:(2^d-1) hsr = wp(packet(d+1,2*b ,n),d+2)'; hsl = wp(packet(d+1,2*b+1,n),d+2)'; lsr = sig(packet(d+1,2*b ,n) ); lsl = sig(packet(d+1,2*b+1,n) ); loterm = (UpDyadLo(lsr,qmf) + lshift(UpDyadLo(lsl,qmf)))/2; hiterm = (UpDyadHi(hsr,qmf) + lshift(UpDyadHi(hsl,qmf)))/2; sig(packet(d,b,n)) = loterm+hiterm; end end x = sig; % % Copyright (c) 1994. David L. Donoho % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function x = IWT2_TI(tiwt,L,qmf) % IWT2_TI -- Invert 2-d translation invariant wavelet transform % Usage % x = IWT2_TI(TIWT,qmf) % Inputs % TIWT translation-invariant wavelet transform table, (3*(J-L)+1)*n by n % L degree of coarsest scale % qmf quadrature mirror filter % Outputs % x 2-d image reconstructed from translation-invariant transform TIWT % % See Also % FWT2_TI, IWT2_TIMedian % [D1,n] = size(tiwt); J = log2(n); D = J-L; % lastx = (3*D*n+1):(3*D*n+n); lasty = 1:n; x = tiwt(lastx,lasty); for d=(D-1):-1:0, l = J-d-1; ns = 2^(J-d); for b1=0:(2^d-1), for b2=0:(2^d-1), index10 = packet(d+1,2*b1,n); index20 = packet(d+1,2*b2,n); index11 = packet(d+1,2*b1+1,n); index21 = packet(d+1,2*b2+1,n); wc00 = [x(index10,index20) tiwt(3*d*n+index10,index20) ; ... tiwt((3*d+1)*n+index10,index20) tiwt((3*d+2)*n+index10,index20)]; wc01 = [x(index11,index20) tiwt(3*d*n+index11,index20) ; ... tiwt((3*d+1)*n+index11,index20) tiwt((3*d+2)*n+index11,index20)]; wc10 = [x(index10,index21) tiwt(3*d*n+index10,index21) ; ... tiwt((3*d+1)*n+index10,index21) tiwt((3*d+2)*n+index10,index21)]; wc11 = [x(index11,index21) tiwt(3*d*n+index11,index21) ; ... tiwt((3*d+1)*n+index11,index21) tiwt((3*d+2)*n+index11,index21)]; x(packet(d,b1,n), packet(d,b2,n)) = ( IWT2_PO(wc00,l,qmf) + .... CircularShift(IWT2_PO(wc01,l,qmf),0,-1) + ... CircularShift(IWT2_PO(wc10,l,qmf),-1,0) + ... CircularShift(IWT2_PO(wc11,l,qmf),-1,-1) ) / 4; end, end end % % Copyright (c) 1995. David L. Donoho and Thomas P.Y. Yu % % % Part of WaveLab Version 802 % Built Sunday, October 3, 1999 8:52:27 AM % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu % function result = CircularShift(matrix, colshift, rowshift) % CIRCULARSHIFT: Circular shifting of a matrix/image, i.e., pixels that get % shifted off one side of the image are put back on the other side. % % result = circularShift(matrix, colshift, rowshift) % % EPS, DJH '96 lastrow = size(matrix, 1); lastcol = size(matrix, 2); result = matrix; % Shift the cols if (colshift>0) result = [result(:,[lastcol-colshift+1:lastcol]) ... result(:,[1:lastcol-colshift])]; else colshift = -colshift; result = [result(:,[colshift+1:lastcol]) ... result(:,[1:colshift])]; end % Shift the rows if (rowshift>0) result = [result([lastrow-rowshift+1:lastrow],:) ; ... result([1:lastrow-rowshift],:)]; else rowshift = -rowshift; result = [result([rowshift+1:lastrow],:) ; ... result([1:rowshift],:)]; end function x = reverse(x) x = x(end:-1:1); function StatWT = TI2Stat(TIWT) % TI2Stat -- Convert Translation-Invariant Transform to Stationary Wavelet Transform % Usage % StatWT = TI2Stat(TIWT) % Inputs % TIWT translation invariant table from FWT_TI % Outputs % StatWT stationary wavelet transform table table as FWT_Stat % % See Also % Stat2TI, FWT_TI, FWT_Stat % StatWT = TIWT; [n,D1] = size(StatWT); D = D1-1; J = log2(n); L = J-D; % index = 1; for d=1:D, nb = 2^d; nk = n/nb; index = [ (index+nb/2); index]; index = index(:)'; for b= 0:(nb-1), StatWT(d*n + (index(b+1):nb:n)) = TIWT(d*n + packet(d,b,n)); end end for b= 0:(nb-1), StatWT((index(b+1):nb:n)) = TIWT(packet(d,b,n)); end % % Copyright (c) 1994. Shaobing Chen % function TIWT = Stat2TI(StatWT) % Stat2TI -- Convert Stationary Wavelet Transform to Translation-Invariant Transform % Usage % TIWT = Stat2TI(StatWT) % Inputs % StatWT stationary wavelet transform table as FWT_Stat % Outputs % TIWT translation-invariant transform table as FWT_TI % % See Also % Stat2TI, FWT_TI, FWT_Stat % TIWT = StatWT; [n,D1] = size(StatWT); D = D1-1; index = 1; for d=1:D, nb = 2^d; nk = n/nb; index = [ (index+nb/2); index]; index = index(:)'; for b= 0:(nb-1), TIWT(d*n + packet(d,b,n)) = StatWT(d*n + (index(b+1):nb:n)); end end for b= 0:(nb-1), TIWT(packet(d,b,n)) = StatWT((index(b+1):nb:n)); end % % Copyright (c) 1994. Shaobing Chen % % % Part of Wavelab Version 850 % Built Tue Jan 3 13:20:40 EST 2006 % This is Copyrighted Material % For Copying permissions see COPYING.m % Comments? e-mail wavelab@stat.stanford.edu