###########Copyright and Notice============== #Copyright since 2001, Bin Han, University of Alberta, #All Rights Reserved #The use of all the following software is permitted for #non-commercial, educational, and research use only. The software #and/or related materials are provided "as-is" without warranty of any #kind including any warranties of performance or merchantability or #fitness for a particular use or purpose or for any purpose #whatsoever, for the licensed product, however used. In no event shall #University of Alberta and/or Bin Han be liable for any damages and/or #costs, including but not limited to incidental or consequential #damages of any kind, including economic damage or injury to property #and lost profits, regardless of whether University of Alberta shall #be advised, have reason to know, or in fact shall know of the #possibility. User bears all risk relating to quality and performance #of the software and/or related materials. #Any use other than non-commercial, educational, or research, or any #redistribution in original or modified form requires prior written #authorization from the copyright holder. #Report errors, mistakes, and bugs to Bin Han at bhan@ualberta.ca. #Send comments/suggestions to Bin Han at bhan@ualberta.ca ##=================================================================== ## MAPLE 6.0 routines for CBC algorithms on constructing 1-dimensional ## biorthogonal multiwavelets with a general dilation factor. #This file consists of a set of useful MAPLE subroutines to handle 1D #matrix masks with a general dilation factor. ### Representation format of matrix mask: # poly(x)=sum(Mask[i]*x^i; i \in Z): #So poly is a matrix of polynomials with the variables (x) reserved #Usage: # the program uses Matrix and Vector from LinearAlgebra. # Do not use matrix and vector from linalg!!! ##D1PrintMask:=proc(poly,choice) ## Print a 1D matrix mask using various formats. ## Input: the matrix of polynomials of the matrix mask. ## Output: Display the matrix mask on screen in various formats. ## and output the matrix mask in array format ## choice<-1, Do not display, output matrix mask in numerical array ## format (6 decimal places) ## choice=-1, Do not display, output matrix mask in exact value array ## choice=0, Display the mask and ouput the exact value form, ## otherwise, Display and output numerical form (6 Digits) ##D1PrintElement:=proc(poly,xp,choice) ## Print the element of a 1D matrix mask at poistion [xp]. ## Input: the matrix of polynomials of the matrix mask, xp,yp ## Output: Display the element of the matrix mask on screen. ## choice>=0, exact format, otherwise, numerical to 6 digits ##D1BiorthEqs:=proc(primal_poly,dual_poly,dilation,orth,choice) ## Find the equations of the (bi)orthogonal relation: ## Primalpoly*\bar{Dualpoly} is an interpolatory mask w.r.t ## Dilation*Z. ## Input: primal_poly, dual_poly, dilation, orth, choice ## orth=0 means dual filter: the primal-dual relation ## orth=1 means orthogonal filter: lowpass-highpass relation ## choice=0, equations in set for exact dual relation ## otherwise, equations in set for dual relation without ## zero-position ##D1SumRuleEqs:=proc(poly,dilation,sr,knownvalue,choice) ## Find the equations for sum rules of order sr w.r.t. ## Dilation*Z. ## Input: poly,dilation,sr,choice,knownvalue ## Output: result -> the equations for sum rules of order sr ## choice:=0, then use standard definition of sum rules for ## scalar mask. That is, sum(Mask[i+Mj](i+Mj)^\mu, j\in Z^s) ## is independent of the integer i. ## choice:=1, then each coset equals the knownvalue. That is, ## sum(Mask[i+Mj](M^(-1)i+j)^\mu, j)=knownvalue[\mu]. ## choice=2, the y vector values are known for matrix mask. ## that is, yvect[\mu]=knownvalue[\mu] ## otherwise, we output the equations of sum rules with ## parameterized y vectors (in yvect) ## Do not assume that the sum of the mask is = |det(dilation)| ## The definition of sum rules is defined in Eqs. (1.6) in ## Bin Han, Approximation properties and construction of ## Hermite Interpolants and Biorthogonal Multiwavelets, ## Journal of Approximation Theory, Vol. 110, (2001) 18-53 ##D1CBCAlgorithm:=proc(poly,dilation,sr,choice) ## Given a matrix mask, find the \tilde y vector values by ## the CBC algorithm ## Input: poly,dilation,sr,choice ## Output: result -> the \tilde y vectors in the sum rules ## in Theorem 3.1 in the following paper ## Bin Han, Approximation properties and construction of ## Hermite Interpolants and Biorthogonal Multiwavelets, ## Journal of Approximation Theory, Vol. 110, (2001) 18-53 ## Choice=0, output the moment for scalar mask: ## h^a(\mu):=sum( \tilde a(i+Mj)(M^(-1)i+j)^\mu/\mu!, j) ## Choice=otherwise, output the \tilde y vectors in Theorem 3.1 ##D1MaskOfSym:=proc(sz,dilation,hermiteorder,interpolating,cc) ## Parametrize the mask using prescribed symmetry group {-1}. ## Input: sz,dilation,hermiteorder,interpolating ## Output: a matrix of polynomial for the parameterized matrix mask ## When interpolating=0, it is not an interpolanting mask ## w.r.t. dilation*Z, otherwise it is an interpolating matrix mask. ##cc is the name of the parameters: cc1, cc2, cc3,... ##D1YVectForHermiteMask:=proc(hmord,sr) ## Given hermiteorder and sum rule order, generate the y vector for ## refinable Hermite mask ##There are also many other useful maple routines in this file. ################################################################# ##Some directly related references: #For symmetry of a refinable function and mask with a dilation matrix # see the papers # Bin Han, Symmetry property and construction of wavelets with a # general dilation matrix, Linear Algebra and Its Applications, # (2001), to appear. # Bin Han, Computing the smoothness exponent of a symmetric # multivariate refinable function, (2001), preprint # Bin Han, Thomas P.-Y. Yu, Bruce Piper, ## Multivariate Refinable Hermite Interpolants, (2002), preprint. #For CBC algorithms on constructing biorthogonal multiwavelets. # CBC stands for Coset By Coset # These program is mainly based on the following paper # Bin Han, Hermite interpolants and biorthogonal multiwavelets # with arbitrary order of vanishing moments (1999) SPIE # Proc. Vol. 3813, pp. 147--161. # The CBC algorithm for the multivariate multiwavelets with # a general dilation matrix was established. # Also, see the following papers for the CBC algorithm # Bin Han, Analysis and Construction of Optimal Multivariate # Biorthogonal Wavelets With Compact Support, SIAM Journal on # Mathematical Analysis, Vol. 31, No.2 (1999/2000), 274--304. # CBC algorithm for interpolatory masks with dilation # matrix 2I_s was first introduced in this paper. # Di-Rong Chen, Bin Han and Sherman D. Riemenschneider, # Construction of Multivariate Biorthogonal Wavelets With # Arbitrary Vanishing Moments, Advances in Computational # Mathematics, Vol. 13 No. 2 (2000), 131-165. # CBC alogrithm for the scalar case with a general dilation # matrix was established. # Bin Han and Rong-Qing Jia, Quincunx Fundamental Refinable # Functions and Quincunx Biorthogonal Wavelets, Mathematics of # Computation, Vol. 71, No. 237, (2002), 165--196. # CBC algorithm for quincunx wavelets was discussed # Bin Han, Construction of multivariate biorthogonal wavelets by # CBC algorithm. Wavelet analysis and multiresolution methods # (Urbana-Champaign, IL, 1999), 105--143, Lecture Notes in Pure # and Appl. Math., 212, Dekker, New York, 2000. # Survey paper on CBC algorithm. # Bin Han, Hermite interpolants and biorthogonal multiwavelets # with arbitrary order of vanishing moments (1999) SPIE # Proc. Vol. 3813, pp. 147--161. # Survey paper on CBC algorithm for 1-dimensional multiwavelets. # Papers can be downloaded at ## http://www.ualberta.ca/~bhan/publ.htm # Program was developed by Bin Han at University of Alberta # Version 1 on January 6, 2002. # Initial Tests of the program have been done on January 6, 2002. # Report bugs, mistakes, errors, comments, suggestions etc. to # Bin Han, bhan@ualberta.ca, htpp://www.ualbert.ca/~bhan #################################################################### ################################################################ ##Main Maple 6.0 routines. ################################################################# read "basicall": default_D1dilation:=2: ##Find the support of a 1D matrix mask. ##Input: the matrix of polynomials of the matrix mask. ##Output: result -> support of mask=[supp[1],supp[2]] D1Support:=proc(poly) local i,j,multi,POLY,result: result:=Vector(2): multi:=Multiplicity(poly): for i from 1 to multi do for j from 1 to multi do POLY:=expand(poly[i,j]): result[1]:=max(result[1],degree(subs({x=1/x},POLY),x)): result[2]:=max(result[2],degree(POLY,x)): od:od: result[1]:=-result[1]: result:=result: end proc: ##Print a 1D matrix mask using various formats. ##Input: the matrix of polynomials of the matrix mask. ##Output: Display the matrix mask on screen in various formats. ## and output the matrix mask in array format ##choice<-1, Do not display, output matrix mask in numerical array ## format (6 decimal places) ##choice=-1, Do not display, output matrix mask in exact value array ##choice=0, Display the mask and ouput the exact value form, ##otherwise, Display and output numerical form (6 Digits) D1PrintMask:=proc(poly,choice) local i,j,ii,supp,multi,POLY,val,printmask,result: multi:=Multiplicity(poly): supp:=D1Support(poly): result:=Array(supp[1]..supp[2],1..multi,1..multi): if (choice>=0) then printmask:=Matrix(multi,(supp[2]-supp[1]+1)*multi): fi: for i from 1 to multi do for j from 1 to multi do POLY:=expand(poly[i,j]): for ii from supp[1] to supp[2] do val:=coeff(POLY,x,ii): if ((choice<(-1))or(choice>0)) then val:=evalf(val,6): fi: result[ii,i,j]:=val: if (choice>=0) then printmask[i,(ii-supp[1])*multi+j]:=val: fi: od: od:od: if (choice>=0) then print(printmask,supp): fi: result:=result: end proc: ##Print the element of a 1D matrix mask at poistion [xp]. ##Input: the matrix of polynomials of the matrix mask, xp,yp ##Output: Display the element of the matrix mask on screen. ##choice>=0, exact format, otherwise, numerical to 6 digits D1PrintElement:=proc(poly,xp,choice) local i,j,multi,POLY,val,result: multi:=Multiplicity(poly): result:=Matrix(multi,multi): for i from 1 to multi do for j from 1 to multi do POLY:=expand(poly[i,j]): val:=coeff(POLY,x,xp): if (choice<0) then val:=evalf(val,6): fi: result[i,j]:=val: od:od: print(result): result:=result: end proc: ##Find the equations of the (bi)orthogonal relation: ## Primalpoly*\bar{Dualpoly} is an interpolatory mask w.r.t ## Dilation*Z. ##Input: primal_poly, dual_poly, dilation, orth, choice ##orth=0 means dual filter: the primal-dual relation ##orth=1 means orthogonal filter: lowpass-highpass relation ##choice=0, equations in set for exact dual relation ##otherwise, equations in set for dual relation without zero-position D1BiorthEqs:=proc(primal_poly,dual_poly,dilation,orth,choice) local i,j,ii,m,supp,multi,POLY,result; multi:=Multiplicity(dual_poly): m:=abs(dilation): #absolute value of the dilation factor POLY:=Matrix(multi,multi): for i from 1 to multi do for j from 1 to multi do POLY[j,i]:=subs({x=1/x},dual_poly[i,j]): od:od: POLY:=.(primal_poly,POLY): supp:=D1Support(POLY): result:={}: if (choice=0) then for i from 1 to multi do POLY[i,i]:=POLY[i,i]-m*orth: od: else for i from 1 to multi do for j from 1 to multi do POLY[i,j]:=expand(POLY[i,j]-coeff(POLY[i,j],x,0)): od:od: fi: for i from 1 to multi do for j from 1 to multi do for ii from supp[1] to supp[2] do #if ii \in dilation*Z, then its coefficient must be 0 if ((abs(ii) mod m)=0) then result:=result union { coeff(POLY[i,j],x,ii) }: fi: od:od:od: result:=result: end proc: ##Find the moment array of a matrix mask. ##Input: poly,dilation,sr ##Output: result -> the array of the moments of a matrix mask: ##When choice=0, output the whole moments ## Ja[\mu]:=sum(Mask[j]*(M^(-1)j)^\mu, j) ## choice=otherwise, output the coset moments ## Ja[gep, \mu]:=sum(Mask[gep+Mj](M^(-1)gep+j)^\mu, j) D1MomentOfMask:=proc(poly,dilation,sr,choice) local i0,i,j,ii,iii,m,multi,POLY,supp,val,mask,tmp,result: m:=abs(dilation): multi:=Multiplicity(poly): supp:=D1Support(poly): mask:=Array(supp[1]..supp[2],1..multi,1..multi): for i from 1 to multi do for j from 1 to multi do POLY:=expand(poly[i,j]): for ii from supp[1] to supp[2] do mask[ii,i,j]:=coeff(POLY,x,ii): od: od:od: #compute the moment array of a matrix mask if (choice=0) then result:=Array(0..(sr-1),1..multi,1..multi): else result:=Array(0..(m-1),0..(sr-1),1..multi,1..multi): fi: for iii from 0 to (sr-1) do #\mu=iii for ii from supp[1] to supp[2] do if(iii=0) then val:=1: else val:=(ii/dilation)^iii: fi: for i from 1 to multi do for j from 1 to multi do tmp:=val*mask[ii,i,j]/(iii!): if (choice=0) then result[iii,i,j]:=result[iii,i,j]+tmp: else i0:=(ii+abs(ii)*m) mod m: result[i0,iii,i,j]:=result[i0,iii,i,j]+tmp: fi: od:od: od:od: result:=result: end proc: ##Find the equations for sum rules of order sr w.r.t. ## Dilation*Z. ##Input: poly,dilation,sr,choice,knownvalue ##Output: result -> the equations for sum rules of order sr ##choice:=0, then use standard definition of sum rules for ## scalar mask. That is, sum(Mask[i+Mj](i+Mj)^\mu, j\in Z^s) ## is independent of the integer i. ##choice:=1, then each coset equals the knownvalue. That is, ## sum(Mask[i+Mj](M^(-1)i+j)^\mu, j)=knownvalue[\mu]. ##choice=2, the y vector values are known for matrix mask. ## that is, yvect[\mu]=knownvalue[\mu] ##otherwise, we output the equations of sum rules with ## parameterized y vectors (in yvect) ##Do not assume that the sum of the mask is = |det(dilation)| ##The definition of sum rules is defined in Eqs. (1.6) in ## Bin Han, Approximation properties and construction of ## Hermite Interpolants and Biorthogonal Multiwavelets, ## Journal of Approximation Theory, Vol. 110, (2001) 18-53 D1SumRuleEqs:=proc(poly,dilation,sr,knownvalue,choice) local i,j,k,ii,jj,m,gep,multi,POLY,supp,val,Ja,tmp,yvect,result; m:=abs(dilation): multi:=Multiplicity(poly): supp:=D1Support(poly): Ja:=D1MomentOfMask(poly,dilation,sr,1): result:=Array(0..sr): if ((multi=1)and((choice=0)or(choice=1))) then #multi=1 for k from 0 to (sr-1) do #\mu=k result[k+1]:={}: for gep from 0 to (m-1) do if (choice=0) then result[k+1]:=result[k+1] union {Ja[gep,k,1,1]-Ja[0,k,1,1]}: else result[k+1]:=result[k+1] union {Ja[gep,k,1,1]-knownvalue[k]}: fi: od: od: else if (choice=2) then yvect:=knownvalue: else yvect:=Array(0..(sr-1),1..multi): jj:=1: for i from 0 to (sr-1) do #\mu=i for ii from 1 to multi do yvect[i,ii]:=evaln(yvect||jj): jj:=jj+1: od:od: fi: for k from 0 to (sr-1) do #\mu=k result[k+1]:={}: for gep from 0 to (m-1) do for i from 1 to multi do tmp:=-((1/dilation)^k)*yvect[k,i]: for ii from 0 to k do #\nu=ii <= \mu for j from 1 to multi do tmp:=tmp+(-1)^ii*Ja[gep,ii,j,i]*yvect[k-ii,j]: od:od: result[k+1]:=result[k+1] union {tmp}: od: od:od: fi: result[0]:={}: for i from 1 to (sr) do result[0]:=result[0] union convert(result[i],set): od: result:=result: end proc: ##Given a matrix mask, find the \tilde y vector values by ##the CBC algorithm ##Input: poly,dilation,sr,choice ##Output: result -> the \tilde y vectors in the sum rules ## in Theorem 3.1 in the following paper ## Bin Han, Approximation properties and construction of ## Hermite Interpolants and Biorthogonal Multiwavelets, ## Journal of Approximation Theory, Vol. 110, (2001) 18-53 ##Choice=0, output the moment for scalar mask: ## h^a(\mu):=sum( \tilde a(i+Mj)(M^(-1)i+j)^\mu/\mu!, j) ##Choice=otherwise, output the \tilde y vectors in Theorem 3.1 D1CBCAlgorithm:=proc(poly,dilation,sr,choice) local i,ii,iii,j,k,m,multi,Ja,Ja0MinusmI,yvect,SE,tmp,tmp1,result: multi:=Multiplicity(poly): m:=abs(dilation): Ja:=D1MomentOfMask(poly,dilation,sr,0): if ((multi=1)and(choice=0)) then if (Ja[0,1,1]=0) then printf(The sum of the primal mask must be nonzero!\n): printf(Program is terminated at subroutine D1CBCAlgorithm!\n): quit: fi: result:=Array(0..(sr-1)): result[0]:=1: for i from 1 to (sr-1) do #\mu=i tmp:=0: for ii from 0 to (i-1) do #\nu=ii < \mu tmp:=tmp+(-1)^(i-ii)*result[ii]*Ja[i-ii,1,1]: od: result[i]:=-tmp/Ja[0,1,1]: od: else result:=Array(0..(sr-1),1..multi): Ja0MinusmI:=matrix(multi,multi): for i from 1 to multi do for j from 1 to multi do Ja0MinusmI[i,j]:=-Ja[0,i,j]: od: Ja0MinusmI[i,i]:=Ja0MinusmI[i,i]+m: od: ##generate y0 vector here tmp1:=kernel(Ja0MinusmI, 'nulldimension'); if (nulldimension <> 1) then printf(There are something wrong with the mask!\n): printf(Program is terminated at subroutine D1CBCAlgorithm!\n): quit; fi: tmp:=tmp1[1]: for i from 1 to multi do result[0,i]:=tmp[i]: od: #generate all other \tilde y vectors (they must be unique) for iii from 1 to (sr-1) do #\mu=iii for i from 1 to multi do Ja0MinusmI[i,i]:=dilation^(1+iii)-Ja[0,i,i]: od: SE:=inverse(Ja0MinusmI): for i from 1 to multi do tmp[i]:=0: for ii from 0 to (iii-1) do #\nu=ii < \mu for k from 1 to multi do tmp[i]:=tmp[i]+dilation^(iii-ii)*Ja[iii-ii,i,k]*result[ii,k]: od:od: od: for i from 1 to multi do result[iii,i]:=0: for j from 1 to multi do result[iii,i]:=result[iii,i]+SE[i,j]*tmp[j]: od: od: od: fi: result:=result: end proc: ##Parametrize the mask using prescribed symmetry group {-1}. ##Input: sz,dilation,hermiteorder,interpolating ##Output: a matrix of polynomial for the parameterized matrix mask ##When interpolating=0, it is not an interpolanting mask ##w.r.t. dilation*Z, otherwise it is an interpolating matrix mask. ##cc is the name of the parameters: cc1, cc2, cc3,... D1MaskOfSym:=proc(sz,dilation,hermiteorder,interpolating,cc) local i,j,ii,jj,iii,jjj,k,ind,invdilation,multi,m,mask,VAR,result: multi:=hermiteorder+1: mask:=Array((-sz)..sz,1..multi,1..multi,0): m:=abs(dilation): result:=Matrix(multi): if (intepolating<>0) then for i from 1 to multi do mask[0,i,i]:=(1/dilation)^(i-1): od: fi: ind:=1: VAR:={}: for ii from 0 to sz do if((interpolating=0)or( ((ii+abs(ii)*m) mod m)<>0)) then for i from 1 to multi do for j from 1 to multi do mask[ii,i,j]:=evaln(cc||ind): mask[-ii,i,j]:=(-1)^(i+j)*mask[ii,i,j]: VAR:=VAR union { evaln(cc||ind) }: ind:=ind+1: od:od: fi: od: for i from 1 to multi do for j from 1 to multi do result[i,j]:=0: for ii from (-sz) to sz do result[i,j]:=result[i,j]+mask[ii,i,j]*(x^ii): od: od:od: save VAR, FreeParameter: result:=result: end proc: ##Given hermiteorder and sum rule order, generate the y vector for ##refinable Hermite mask D1YVectForHermiteMask:=proc(hmord,sr) local i,j,result: if(hmord=0) then result:=Array(0..(sr-1),0): result[0]:=1: else result:=Array(0..(sr-1),1..(hmord+1),0): for i from 0 to hmord do result[i,i+1]:=1: od: fi: result:=result: end proc: