Discrete Linear Canonical Transform Based on Hyperdifferential Operators

Discrete Linear Canonical T ransform Based on Hyp erdifferen tial Op erators Aykut Ko¸ c, Burak Bartan, and Haldun M. Ozaktas, ∗ † ‡ F ebruary 14, 2019 Abstract Linear canonical transforms (LCTs) are of imp or- tance in man y areas of science and engineering with man y applications. Therefore a satisfactory discrete implemen tation is of considerable interest. Although there are metho ds that link the samples of the input signal to the samples of the linear canonical trans- formed output signal, no widely-accepted definition of the discrete LCT has been established. W e in- tro duce a new approach to defining the discrete lin- ear canonical transform (DLCT) b y employing oper- ator theory . Op erators are abstract entities that can ha ve b oth contin uous and discrete concrete manifes- tations. Generating the con tinuous and discrete man- ifestations of LCTs from the same abstract op erator framew ork allo ws us to define the con tinuous and dis- crete transforms in a structurally analogous manner. By utilizing h yp erdifferen tial op erators, w e obtain a DLCT matrix which is totally compatible with the theory of the discrete F ourier transform (DFT) and its dual and circulan t structure, which mak es further analytical manipulations and progress p ossible. The prop osed DLCT is to the contin uous LCT, what the DFT is to the contin uous F ourier transform (FT). The DLCT of the signal is obtained simply by mul- tiplying the v ector holding the samples of the input signal by the DLCT matrix. ∗ Aykut Ko¸ c (Corresponding Author) is with ASELSAN Research Cen ter, Ank ara, T urkey , e-mail: aykutk o c@aselsan.com.tr. † Burak Bartan is with Electrical Engineering Department, Stanford Univ ersity , Stanford, California ‡ Haldun M. Ozaktas is with Electrical Engineering Depart- ment, Bilkent University , TR-06800 Bilk ent, Ank ara, T urkey 1 In tro duction Linear canonical transforms (LCTs) are a family of linear in tegral transforms with three parameters, [25, 48, 54, 75]. The family of LCTs is a generaliza- tion of many imp ortant transforms such as the frac- tional F ourier transform (FR T), chirp m ultiplication (CM), chirp con volution (CC), and scaling op era- tions. F or certain v alues of the three parameters, the LCT reduces to these transforms or their com- binations. LCTs ha ve sev eral applications in signal pro cessing [25] and computational and applied math- ematics [19, 37], including fast and efficient optimal filtering [7], radar signal pro cessing [15, 16], sp eec h pro cessing [61], image represen tation [1], and image encryption and watermarking [41, 60, 67], to mention a small sample of published w orks. LCTs hav e also b een extensively studied for their applications in op- tics [4, 8–10, 48, 65, 66], electromagnetics, and classical and quantum mec hanics [25, 34, 42, 75]. In optical con texts, LCTs are commonly referred to as quadratic-phase in tegrals or quadratic-phase sys- tems [9, 47]. The so-called AB C D systems widely used in optics [29] are also represented by linear canonical transforms. They hav e also b een referred to b y other names: generalized Huygens integrals [65], generalized F resnel transforms [33, 52], special affine F ourier transforms [2, 3], extended fractional F ourier transforms [32], and Moshinsky-Quesne transforms [75]. Tw o-dimensional (2D) LCTs and complex- parametered LCTs (CLCTs) ha ve also been discussed in the literature, [21, 39, 40, 62]. Bilat- eral Laplace transforms, Bargmann transforms, 1 Gauss-W eierstrass transforms, [73 – 75], fractional Laplace transforms, [63, 69], and complex-ordered FR Ts [11, 12, 64, 71] are all sp ecial cases of CLCTs. The establishment of a discrete framew ork is es- sen tial to the deploymen t of LCTs in applications. There is considerable work on discrete or finite forms of fractional F ourier transforms, and, to a lesser de- gree, discrete or finite linear canonical transforms. Being one of the most imp ortan t sp ecial cases of LCTs, discretization and discrete versions of frac- tional F ourier transforms ha ve b een well studied and established [5, 6, 14, 20, 57 – 59, 70, 76, 78 – 80, 82]. As for the discretization or digital computation of LCTs, there are many approaches present in the lit- erature, [13, 26 – 28, 30, 31, 38, 46, 47, 53, 55, 56, 68, 72, 83 – 85]. Some of these [13, 27, 28, 38, 47, 55, 68, 83, 84] n umerically compute the con tinuous in tegral and es- tablish a direct mapping b et ween the samples of the con tinuous input function and the samples of the LCT-transformed contin uous output function. The metho ds in [27, 53, 68, 83, 84] directly con v ert the LCT in tegral to a summation and [13, 28, 38, 46, 47, 55] mak e use of decompositions in to elementary build- ing blo c ks. Moreov er, some approac hes focus on defining a discrete LCT (DLCT), which can then b e used to numerically approximate contin uous LCTs, in the same wa y that the discrete F ourier transform (DFT) is used to appro ximate contin uous F ourier transforms [26, 30, 31, 46, 53, 56, 68, 72, 85]. Algorithms in [46, 53, 68] also numerically appro ximate the con tin- uous LCTs in the same wa y the DFT approximates the contin uous FT. Based on the DLCT definition prop osed in [53], Refs. [27] and [28] propose efficient n umerical computation algorithms. Ref. [53] also in- cludes a comparison of the prop erties satisfied b y def- initions of DLCTs prop osed up to that date. Despite these works, no single definition has b een widely established as the definition of the DLCT. In this pap er, we presen t a different approach based on h yp erdifferen tial op erator theory [45, 48, 49, 75, 81], to obtain a definition of the DLCT. Wh y do w e pro- p ose to use op erator theory? Most approaches to discretization are naturally based on sampling of the con tinuous entities. How ever, sampling often do es not lead to a clean, discrete transform definition that satisfies op erational form ulas and exhibits desirable analytical prop erties such as unitarity and preserv a- tion of the group structure. So if our purp ose is not to merely numerically compute a contin uous trans- form, but to obtain a self-consistent discrete trans- form definition, it often turns out to be insufficient. A purely n umerical metho d can compute the con tin- uous transform accurately , but it do es not provide us with a definition on which further manipulation can b e done, and theoretical progress can build up on. W e w ant a discrete definition that is as analogous to the con tinuous definition as p ossible. (This is satisfied b y the discrete F ourier transform (DFT) and that is wh y the DFT is so established.) Ho w do es op erator theory help? Op erators are ab- stract en tities that can ha ve b oth contin uous and dis- crete concrete manifestations. Thus if we b egin from a contin uous entit y and can appropriately deduce the abstract operator underlying that entit y , then, that can form a basis for defining its discrete v er- sion. Since b oth the contin uous and discrete versions are based on the same abstract op erator, they can b e exp ected to exhibit similar structural c haracteris- tics and op erational prop erties to the exten t pos sible. The structure of relationships b et ween different en- tities can also b e preserved and can b e exp ected to mirror the relationships b et ween the abstract op er- ators. Thus we can obtain discrete entities that are not merely numerical appro ximations, but which ex- hibit desirable analytical and op erational properties. This is the rationale of the presen t paper. Our definition of the discrete LCT will b e pre- sen ted in the form of a matrix of size N × N which, up on multiplication, pro duces the DLCT of a dis- crete and finite signal of length N , expressed as a column v ector. The main difference from earlier ap- proac hes is that the definition is based on hyperdif- feren tial forms of the discrete coordinate multiplica- tion and differentiation op erators, which w e carefully define so that they are strictly F ourier duals related through the DFT matrix. Our definition provides a self-consistent, pure, and elegan t definition of the DLCT which is fully compatible with the theory of the discrete F ourier transform and its dual and cir- culan t structure. By self-consisten t we mean that the relations b et ween discrete entities should mirror those betw een contin uous en tities as muc h as p ossi- 2 ble, e.g. if the co ordinate mu ltiplication and differ- en tiation operators are dual in the con tinuous case, they should also be so in the discrete case. The dis- crete LCT should b e built up on these t wo op erators in the same wa y that the contin uous LCT is, and so forth. By dualit y w e mean that a kind of sym- metry b et w een the tw o domains is exactly satisfied (e.g. co ordinate m ultiplication in one domain is dif- feren tiation in the other, translation in one domain is phase multiplication in the other, etc.). All the dual prop erties of the F ourier transform (suc h as those in paren thesis ab o ve) can be derived from the dualit y of U and D [48], so first and foremost, this duality m ust b e maintained. One of the most imp ortan t features of our approach is that our definition maintains this structure by treating b oth domains totally symmet- rically . The pap er is organized as follo ws: Section 2 re- views the preliminaries and the definition and im- p ortan t prop erties of LCTs. Section 3 describ es the theory and deriv ations for the prop osed DLCT. The- oretical discussions on defining a discrete LCT and the properties of such a definition that need to ex- ist are given in Section 4. In Section 5, n umerical examples and comparisons are pro vided. Lastly , we conclude in Section 6. There is also an App endix in whic h w e hav e provided some pro ofs, necessary fun- damen tal information, justifications and implemen- tation details that are needed for the deriv ations in Section 3. 2 Preliminaries 2.1 Linear Canonical T ransform LCTs are unitary transforms sp ecified by a 2 × 2 pa- rameter matrix L . Because the determinant of L is required to be equal to 1, an LCT can also be uniquely sp ecified by three indep enden t parameters, often denoted b y α, β , γ . The elements A, B , C, D of the 2 × 2 matrix and α, β , γ are related b y: L =  A B C D  =  γ β 1 β − β + αγ β α β  =  α β − 1 β β − αγ β γ β  − 1 . (1) W e can define an LCT through either the parameter set ( A, B , C, D ) with the condition that AD − B C = 1 or the parameter set ( α, β , γ ). In this pap er, we restrict ourselves to the case where the parameters in both sets are all real. The definition of the LCT as a linear in tegral transform, using the second set of parameters, can b e written as: C L f ( u ) = p β e − iπ / 4 Z ∞ −∞ exp  iπ ( αu 2 − 2 β uu 0 + γ u 0 2 )  f ( u 0 ) du 0 . (2) Ev ery triplet ( α, β , γ ) corresp onds to a differen t LCT. W e denote the LCT op erator using C L where the sub- script L denotes the 2 × 2 parameter matrix. 2.2 Imp ortan t Prop erties The utility of the parameter set ( A, B , C, D ) is b est appreciated up on observing the concatenation prop- ert y: If any tw o LCTs are concatenated (applied one after the other), the resulting op eration is also an LCT whose 2 × 2 matrix is the pro duct of the 2 × 2 matrices of the t wo original LCTs. This can be stated as: C L f ( u ) = C L 1 C L 2 f ( u ) , (3) where L = L 1 L 2 . An imp ortant sp ecial case of this prop ert y is the rev ersibility prop ert y . It basically states that the 2 × 2 matrix for the inv erse of an LCT is again an LCT whose 2 × 2 matrix is the matrix inv erse of the original LCT: C L 2 C L 1 f ( u ) = f ( u ) , (4) if L 2 = L − 1 1 . 2.3 Sp ecial Linear Canonical T rans- forms W e now giv e some sp ecial transforms and operations, whic h are all sp ecial cases of LCTs. 3 2.3.1 Scaling The parameter matrix for the scaling operation is as follo ws L M =  M 0 0 1 M  =  1 M 0 0 M  − 1 . (5) F unctionally it can b e defined in the following w ay: C L M f ( u ) = M M f ( u ) = r 1 M f  u M  . (6) 2.3.2 F ractional F ourier T ransform The F ractional F ourier transform (FR T) is the gen- eralized version of the F ourier transform (FT). It has the following parameter matrix: L F a lc =  cos θ sin θ − sin θ cos θ  =  cos θ − sin θ sin θ cos θ  − 1 , (7) where θ = π a/ 2 and a is the fractional order. When a = 1, the FR T reduces to the FT. (It should b e noted that there is a sligh t difference b etw een the FR T thus defined ( F a lc ) and the more commonly used definition of the FR T ( F a ), [48].) The a th order fractional F ourier transform F a of the function f ( u ) ma y be defined as [48]: F a f ( u ) = Z ∞ −∞ K a ( u, u 0 ) f ( u 0 ) du 0 , K a ( u, u 0 ) = A θ exp  iπ ( u 2 cot θ − 2 uu 0 csc θ + u 0 2 cot θ )  , A θ = exp( − iπ sgn(sin θ ) / 4 + iθ / 2) | sin θ | 1 / 2 (8) 2.3.3 Chirp Multiplication The parameter matrix for the chirp m ultiplication op eration is L Q q =  1 0 − q 1  =  1 0 q 1  − 1 . (9) The c hirp multiplication op eration can be expressed as C Q q f ( u ) = Q q f ( u ) = exp( − iπ qu 2 ) f ( u ) . (10) Corresp onding formulas for chirp conv olution may b e found in [48]. 3 Discrete Linear Canonical T ransforms W e now presen t our developmen t of the DLCT based on hyperdifferential op erator theory . Our approach is based on decomp osing the LCT into simpler parts, finding the discrete versions of these parts b y using op erator theory , and then m ultiplying those to obtain the final DLCT matrix. Although there are several wa ys to decomp ose the LCT [38], here we c ho ose the Iw asaw a decomp osition since it includes a greater num b er of sp ecial LCTs than other decomp ositions, providing the opp ortu- nit y to discuss their hyperdifferential forms. The metho d of using hyperdifferential op erators outlined here can also b e applied to other decompositions. 3.1 The Iw asaw a Decomp osition The linear canonical transform (LCT) operator C L can b e expressed as com binations of other simpler op erators in many wa ys. Using scaling M M , chirp m ultiplication Q q and fractional F ourier F a op era- tors, it is p ossible to construct any linear canoni- cal transform. The Iw asaw a decomp osition w e will emplo y , breaks do wn an arbitrary LCT into a frac- tional F ourier transform follow ed by scaling follow ed b y c hirp multiplication, and can b e written in op er- ator notation as follo ws [25]: C L = Q q M M F a lc , (11) When each op erator is c haracterized b y their 2 × 2 LCT parameter matrix, the decomp osition lo oks lik e L =  A B C D  =  γ β 1 β − β + αγ β α β  =  1 0 − q 1   M 0 0 1 / M   cos aπ / 2 sin aπ / 2 − sin aπ / 2 cos aπ / 2  (12) 4 where a , q , M must b e chosen as: M =  p 1 + γ 2 /β , γ ≥ 0 , − p 1 + γ 2 /β , γ < 0 , (13) q = γ β 2 1 + γ 2 − α, (14) a = 2 π arccot γ . (15) This decomp osition can break down any arbitrary linear canonical transform in to a cascade of elemen- tary op erations. Our approac h will b e to find the N × N discrete transform matrix for each of these three operations and m ultiply them to obtain the dis- crete LCT matrix. 3.2 The Hyp erdifferen tial F orms The term hyperdifferential refers to ha ving differen- tial op erators in an exp onen t. In the LCT con text, w e only hav e second order co ordinate multiplication and differentiation op erators in the exp onent. Op er- ators represen ting an arbitrary LCT or all of its sp e- cial cases can b e generated b y exp onentiating these second order operators and these constitute the hy- p erdifferen tial forms of these transforms. There is corresp ondence among the integral transforms, hy- p erdifferen tial operators and the 2x2 parameter ma- trices that are given in the preliminaries section. An LCT can be represented by any one of these mathe- matical ob jects. More details can b e found in [75]. It is well established that the chirp m ultiplication op erator Q q , the scaling op erator M M , and the frac- tional F ourier transform op erator F a lc can all b e writ- ten in hyperdifferential forms as follo ws: [48, 75]: Q q = exp  − i 2 π q U 2 2  , (16) M M = exp  − i 2 π ln ( M ) U D + DU 2  , (17) F a lc = exp  − iaπ 2 U 2 + D 2 2  , (18) where U and D are the co ordinate multiplication and differen tiation op erators, resp ectiv ely . W e see that all three of the op erators we are working with can b e expressed in terms of these tw o building blo cks, whose contin uous manifestations are: U f ( u ) = uf ( u ) (19) D f ( u ) = 1 i 2 π d f ( u ) du , (20) where the ( i 2 π ) − 1 is included so that U and D are precisely F ourier duals (the effect of either in one do- main is its dual in the F ourier domain). This duality can b e expressed as follo ws: U = F D F − 1 . (21) 3.3 The Discrete Linear Canonical T ransform Our approach is based on requiring that, to the exten t p ossible, all the discrete entities we define observ e the same structural relationships as they do in abstract op erator form. W e wan t a discrete definition that is as analogous to the con tinuous definition as p ossible. T o ensure this, we define the discrete LCT and its sp ecial cases as the discrete manifestations of Eq. 11, Eq. 16, Eq. 17 and Eq. 18, with the abstract op era- tors b eing replaced by matrix op erators. This can b e written as follows: C L = Q q M M F a lc . (22) Q q = exp  − i 2 π q U 2 2  . (23) M M = exp  − i 2 π ln ( M ) UD + DU 2  . (24) F a lc = exp  − iaπ 2 U 2 + D 2 2  . (25) Note that exp() in the ab o ve equations are matrix ex- p onen tials and ho w they are computed is discussed in App endix C. Thus the discrete LCT matrix is given b y C L = exp  − i 2 π q U 2 2  × exp  − i 2 π ln ( M ) UD + DU 2  exp  − iaπ 2 U 2 + D 2 2  . (26) 5 The discrete LCT matrix is defined as the pro duct of the FR T, scaling, and chirp m ultiplication matrices, all of whic h are defined in terms of the U and D matrices. T o get the DLCT of a function of a discrete v ariable, we just need to write it as a column vector and multiply it with the DLCT matrix C L . Th us it is seen that all rests on the differen tiation and co ordinate m ultiplication matrices D and U and computation of the matrix exp onen tials in Eq. 26. Th us, w e mov e on to how to obtain the U and D matrices. F or signals of discrete v ariables, the closest thing to differen tiation is finite differencing. Consider the follo wing definition: ˜ D h f ( u ) = 1 i 2 π f ( u + h/ 2) − f ( u − h/ 2) h . (27) If h → 0, then ˜ D h → D , since in this case the right- hand side approaches ( i 2 π ) − 1 d f ( u ) /du . Therefore, ˜ D h can b e interpreted as a finite difference operator. No w, using f ( u + h ) = exp( i 2 π h D ) f ( u ), whic h is another established result in op erator theory [48, 75], w e express Eq. 27 in hyperdifferential form: ˜ D h = 1 i 2 π e iπ h D − e − iπ h D h = 1 i 2 π 2 i sin( π h D ) h = sinc( h D ) D . (28) Note that if w e let h → 0 in the last equation and tak e the limit, we can v erify that ˜ D h → D from here as well. No w, we turn our attention to the task of defining ˜ U h . It is tempting to define the discrete version of the co ordinate m ultiplication matrix by simply form- ing a diagonal matrix with the diagonal entries b eing equal to the co ordinate v alues. How ever, up on closer insp ection we hav e decided that this could not be tak en for gran ted. In order to obtain the most self- consisten t form ulation p ossible, we must b e sure to main tain the structural symmetry b et w een U and D in all their manifestations. Therefore, we choose to define ˜ U h suc h that it is related to U , in exactly the same wa y as ˜ D h is related to D : ˜ U h = sinc( h U ) U , (29) from which we can observ e that as h → 0, we ha ve ˜ U h → U , as should b e. How ever, b ey ond that, it is also p ossible to sho w that, ˜ U h , when defined like this, satisfies the same duality expression Eq. 21 satisfied b y U and D : ˜ U h = F ˜ D h F − 1 . (30) T o see this, substitute ˜ D h in this equation: ˜ U h = F  1 i 2 π 2 i sin( π h D ) h  F − 1 = 1 i 2 π 2 i sin( π h U ) h = sinc( h U ) U . (31) When acting on a contin uous signal f ( u ), the op era- tor U b ecomes ˜ U h f ( u ) = 1 π sin( π hu ) h f ( u ) . (32) W e observe that the effect is not merely multiplying with the co ordinate v ariable. Had w e defined ˜ U h suc h that it corresp onds to m ultiplication with the co ordi- nate v ariable, we would ha ve destro yed the symmetry and duality b et ween U and D in passing to the dis- crete world. No w, by sampling Eq. 32, w e can obtain the matrix op erator to act on finite discrete signals. The sample p oin ts will be taken as u = nh to finally yield the U matrix defined as: U mn = ( √ N π sin  π N n  , for m = n 0 , for m 6 = n . (33) As alwa ys, the v alue of N should b e determined based on the time/space and frequency exten t of the signal, along with the required accuracy [23, 38, 50, 51]. F ur- ther detail is pro vided in Section 3.5. The matrix D , on the other hand, can b e calculated in terms of U by using the discrete version of the dualit y relation given in Eq. 21: D = F − 1 UF , (34) in whic h F is the matrix represen ting the unitary dis- crete F ourier transform (DFT) matrix. The elements F mn of the N -p oint unitary DFT matrix F can b e written in terms of W N = exp( − j 2 π / N ) as follows: F mn = 1 √ N W mn N . 6 When all is put together, the LCT of a signal x [ n ] of length N , represen ted b y the column vector x , is then computed by C L x , yielding an N × 1 output. F urther details of the developmen t of the U and D matrices and their applications may b e found in [36], whic h together with the presen t w ork, not only estab- lish a form ulation of these op erators that is fully con- sisten t with the theory of the DFT and its circulant structure, but also pa ve the w ay for the utilization of op erator theory in deriving other more sophisti- cated discrete operations. W e b eliev e these w orks are the first to apply op erator theory in defining discrete transforms. 3.4 Unitarity of the Discrete Linear Canonical T ransform One of the most essen tial properties of the kind of discrete transforms we are working with is unitarity . This leads to Parsev al type relationships and mani- fests itself as energy or pow er conserv ation in ph ysical applications. Here we pro ve that the prop osed DLCT definition is unitary by showing that the matrix C L giv en in Eq. 22 and more explicitly in Eq. 26 is unitary . Theorem 1. The discr ete LCT define d in Eq. 26 is unitary, with M , q , a chosen ac c or ding to Eqs. 13, 14, 15, and U and D define d ac c or ding to Eqs. 33 and 34. Before pro ceeding with the pro of, we first recall some fundamental definitions: A matrix A is said to b e Hermitian when A = A H holds, where A H denotes the conjugate transp ose of A , and is said to b e unitary when A − 1 = A H . Since C L is defined as the pro duct of three matrices, showing that each of them is unitary will suffice to show that C L is unitary . U and D are the fundamental matrices that giv e rise to those three components. W e will first show that these matrices are Hermitian. F rom that it will follow that the three m ultiplied matrices are all unitary . Theorem 2. The matric es U and D ar e Hermitian and the matric es define d in Eqs. 23, 24, 25 ar e uni- tary. Theorem 2 is prov ed in the App endix A from which Theorem 1 follows. 3.5 Discretization, Sampling and In- dexing W e introduce discretization by replacing the con tin- uous deriv ativ e with a finite difference, such that, as the finite interv al go es to zero, it approaches the con- tin uous deriv ative. Remem b ering that exp onen tia- tion etc. can b e expressed as p ow er series, the full LCT developmen t is then based on the following op- erations on this finite difference op eration: in version, fractional and ordinary F ourier transformation, re- p eated application, m ultiplication with a scalar and addition. No w, as the finite difference goes to a deriv ative, similar will hold for its rep eated applica- tions, as w ell as scalar m ultiplied and added v ersions. Lik ewise, w e know that the DFT appro ximates the con tinuous F ourier transform more and more closely as the sampling in terv al is reduced, so if this op era- tion is in succession with finite differencing, the re- sulting limit will b e the succession of F ourier trans- formation and con tinuous differen tiation. Similar ap- plies to fractional F ourier transformation, of whic h in version is a sp ecial case. In this paper we deal with finite-length signals of a discrete (integer) v ariable. (W e could equiv alently think of them as being defined on a circulan t domain, whic h w ould not make a difference in our argumen ts.) The length of our signal v ectors will b e denoted b y N . When N is even, they will b e defined on the interv al of in tegers [ − N 2 , N 2 − 1], and when N is odd, they will b e defined on the interv al of in tegers [ − N − 1 2 , N − 1 2 ]. W e will also consider an alternativ e, less-common ap- proac h based on the device of using “half integers.” In this approach, the domain is defined as the in ter- v al of unit-spaced half in tegers [ − N 2 + 0 . 5 , N 2 − 1 + 0 . 5] for ev en N and [ − N − 1 2 − 0 . 5 , N − 1 2 − 0 . 5] for o dd N . Although not very usual, there is nothing unnatu- ral ab out this w ay of indexing signals of a discrete v ariable; it is merely a particular wa y of b o okk eep- ing. Note that the indices are still spaced by unity , and there is merely a shift by 0 . 5 with the purpose of making the interv al symmetrical around the origin 7 when N is ev en (with the consequence that symme- try is lost when N is o dd). A few examples of w orks considering this w ay of indexing are [17, 22, 44, 70]. Consisten t with this literature, we will refer to the former approac h as the or dinary DFT and refer to the latter one, in which we use ”half integers”, as the c enter e d DFT. The DLCT deriv ation pro cedure w e presen ted has b een carefully written in a manner that it is consistent with b oth approac hes. Readers in terested in further details on this issue ma y refer to [36]. Ho w the n umber of samples N should b e chosen will b e determined by factors suc h as the temp o- ral or spatial exten t of the signal, the frequency ex- ten t of the signal and therefore the time- or space- bandwidth product. It will also dep end on the pre- cision with which the results need to be computed in that application. The choice of N is exogenous to our metho d. Nev ertheless, for completeness, let us elab orate on how the num b er of samples N is c ho- sen. If the temp oral or spatial exten t is ∆ x and the double-sided frequency extent is ∆ ν , then we should b e sampling with an in terv al of 1 / ∆ ν , which means ∆ x/ (1 / ∆ ν ) = ∆ x ∆ ν samples. W e call this n umber of samples N , the time- or space- bandwidth pro duct. If appropriate normalization as describ ed in [38] is ap- plied so that the time/space exten t and the frequency exten t are made equal in a dimensionless space, it follo ws that we should sample ov er an extent √ N with sampling in terv al h = 1 / √ N . Th us as w e in- crease N , we will b e making h smaller and smaller. Consequen tly , the finite difference op erator in Eq. 27 approac hes a contin uous deriv ative and the finite co- ordinate multiplication op erator will approac h the con tinuous co ordinate multiplication op erator. The matrix in Eq. 33 will approach U mn = n/ √ N , corre- sp onding to samples of contin uous co ordinate multi- plication. Since all our operators, including the LCT, are defined in terms of co ordinate multiplication and differen tiation through smo oth exp onen tial functions, they will all approach their con tinuous counterparts. 4 Discussions Con tinuous unitary LCTs represen ted b y the pa- rameter matrices L form the real symplectic group S p (2 , R ) with three indep enden t parameters [43]. The desirable prop erties of a discrete LCT mirror those of the contin uous LCT: unitarit y , preserv ation of group structure as expressed by the concatena- tion prop ert y (and its special case reversibilit y), re- duction to imp ortan t sp ecial cases and inv erses of sp ecial cases, and some satisfactory approximation of the contin uous transform. How ev er, a theorem from group theory [35, 77] precludes realization of this ideal: It is theoretically imp ossible to discretize all LCTs with a finite n umber of samples such that they are b oth unitary and they preserve the group struc- ture [35, 77]. More on the group-theoretical prop er- ties of LCTs can b e found in [48, 75, 77]. That said, no unitary DLCT definition can exhibit exact concatenation/reversibilit y prop erties. How- ev er, if the prop osed definition is to hav e practi- cal use, we can exp ect that these properties are at least approximately satisfied. In Section 3.4, we the- oretically prov ed that our prop osed DLCT is uni- tary , so that it cannot exactly satisfy the concate- nation/rev ersibility property . Therefore, in the next section, we will n umerically sho w that the concate- nation and reversibilit y prop erties are satisfied with a reasonable accuracy . W e will also show that, re- gardless of concatenation, the discrete transform pro- vides a reasonable appro ximation to the con tinuous LCT. Before mo ving on, it needs to be noted that our definition, by construction, reduces to the identit y , F ourier and fractional F ourier transforms, c hirp m ul- tiplication, and magnification (scaling). This result can b e trivially obtained by substituting the com bi- nation of v alues leading to the special cases for the parameters a , M , and q in Eq. 26. 5 Numerical Results and Com- parisons W e will numerically explore three different asp ects of the proposed DLCT definition: (i) approximation of the con tinuous LCT, (ii) concatenation of multiple 8 transforms, and (ii) reversibilit y . W e will carry out n umerical tests regarding these asp ects of the pro- p osed DLCT definition. As the example input functions, the dis- cretized versions of the c hirp ed pulse function exp( − π u 2 − iπ u 2 ), denoted F1, the trap ezoidal func- tion 1 . 5tri( u/ 3) − 0 . 5tri( u ), denoted F2 (tri( u ) = rect( u ) ∗ rect( u )), rectangular pulse function rect( u ), denoted F3, and the damp ed sine function exp( − 2 | u | ) sin(3 π u ), denoted F4, are used. The n um- b er of samples N are tak en as 256 and 1024 for tw o sets of numerical simulations. F our transforms, de- noted by T1, T2, T3, and T4, are considered, with parameters ( α, β , γ ) = ( − 3 , − 2 , − 1), ( − 0 . 8 , 3 , 1), ( − 1 . 8 , − 1 . 75 , − 1 . 3), and (0 . 3 , − 1 . 6 , − 0 . 9), resp ec- tiv ely . The LCTs T1, T2, T3 and T4 of the func- tions F1, F2, F3 and F4 hav e b een computed b oth b y the presen ted DLCT and b y a highly inefficient brute force numerical approach whic h is tak en as a reference. Throughout our numerical comparisons we use percentage mean squared error (MSE) as the p er- formance metric. It is defined as the energy of the difference normalized b y the energy of the reference, expressed as a p ercen tage. 5.1 Approximation of the Con tinuous LCT In this subsection, we fo cus on ho w well our metho d appro ximates the contin uous LCT. The “true” con- tin uous LCT of the original function is obtained b y highly inefficient brute force numerical in tegra- tion of the con tin uous LCT. The resulting percentage MSE scores, for both or dinary and c enter e d sampling sc hemes, turn out to b e giving v ery similar results, are tabulated in T able 1. Plots for some examples for the resulting DLCTs (T1 of F1, T2 of F2, T3 of F3 and T4 of F4) and the corresp onding references obtained by the brute force n umerical metho d hav e b een presen ted for both real and imaginary parts of the signals in Fig. 1. Although we use the same t wo v alues of N for all the signals we consider for fair comparison, normally the v alue of N should b e c hosen according to the exten t of the signals in b oth the time/space and fre- quency domains. The error is primarily determined b y how muc h of the signal falls outside of the extents implied by the chosen v alue of N . F or example, for F1, whic h has a v ery rapidly decaying Gaussian en- v elop e, very little falls outside so the errors are muc h smaller than for the others. In those cases where the results are not sufficiently accurate for the applica- tion at hand, it is p ossible to obtain higher accuracy b y increasing N. 5.2 Concatenation In order to test how well the concatenation prop ert y is satisfied, w e employ the follo wing pro cedure. Let us consider T1 and T2 as an example: First deriv e the DLCT matrices C L 1 and C L 2 for T1 and T2 sep- arately , following the pro cedure giv en in Section 3. Then, by using Eq. 1, we calculate the 2 × 2 LCT parameter matrices L 1 and L 2 for T1 and T2. Mul- tiplying these t wo matrices by using Eq. 3, w e ob- tain the 2 × 2 parameter matrix of the concatenated system L 12 = L 2 L 1 . Then, we obtain C L 12 from L 12 , again by using our proposed DLCT pro cedure. Finally , we compare the result of applying the con- catenated transform matrix C L 12 directly with the result of applying C L 1 and C L 2 consecutiv ely . More precisely , w e compare C L 12 x with C L 2 C L 1 x where a signal x [ n ] of length N is represented b y the col- umn v ector x . The resulting MSE differences are tabulated in T able 2 for sev eral such concatenations among T1, T2, T3, and T4. The ordinary sampling sc heme is used in these numerical calculations. 5.3 Reversibilit y T o test the rev ersibility prop ert y numerically , w e fol- lo w a similar pro cedure as in concatenation. This time the second LCTs in the cascade are the inv erses of the first ones. F or example, we compare x with C L − 1 1 C L 1 x . Again the ordinary sampling scheme is used in these calculations and the resulting MSE dif- ferences are tabulated in T able 2. 9 T able 1: Percen tage MSE Errors for Different F unctions and T ransforms (for b oth ordinary and centered sc hemes) Input N T1 (ord.) T2 (ord.) T3 (ord.) T4 (ord.) T1 (cen t.) T2 (cent.) T3 (cen t.) T4 (cent.) F1 256 9 . 82 × 10 − 4 4 . 72 × 10 − 3 6 . 78 × 10 − 4 3 . 93 × 10 − 2 9 . 82 × 10 − 4 4 . 71 × 10 − 3 6 . 78 × 10 − 4 3 . 93 × 10 − 2 1024 6 . 40 × 10 − 5 2 . 76 × 10 − 4 4 . 26 × 10 − 5 2 . 49 × 10 − 3 6 . 40 × 10 − 5 2 . 76 × 10 − 4 4 . 26 × 10 − 5 2 . 49 × 10 − 3 F2 256 4 . 31 10 . 6 1 . 95 6 . 65 4 . 31 10 . 6 1 . 96 6 . 65 1024 0 . 32 0 . 87 0 . 13 0 . 46 0 . 32 0 . 87 0 . 13 0 . 46 F3 256 2 . 49 1 . 55 2 . 84 2 . 85 2 . 02 1 . 45 2 . 37 2 . 66 1024 1 . 09 0 . 75 1 . 40 1 . 44 1 . 10 0 . 85 1 . 34 1 . 50 F4 256 1 . 34 0 . 64 2 . 29 6 . 77 1 . 35 0 . 63 2 . 30 6 . 79 1024 9 . 43 × 10 − 2 4 . 38 × 10 − 2 0 . 16 0 . 49 9 . 44 × 10 − 2 4 . 38 × 10 − 2 0 . 16 0 . 49 T able 2: Percen tage MSE Errors for Different Concatenations and Inv erses Input N T1-T2 T3-T4 T3-T1 T3-T2 T1-T1 − 1 T3-T3 − 1 F1 256 1 . 32 × 10 − 2 2 . 78 × 10 − 3 1 . 55 × 10 − 3 4 . 10 × 10 − 3 5 . 85 × 10 − 3 9 . 64 × 10 − 4 1024 6 . 82 × 10 − 4 1 . 71 × 10 − 4 9 . 58 × 10 − 5 2 . 79 × 10 − 4 3 . 85 × 10 − 4 6 . 29 × 10 − 5 F2 256 17 . 7 0 . 34 0 . 35 2 . 99 1 . 77 0 . 49 1024 1 . 64 2 . 47 × 10 − 2 2 . 43 × 10 − 2 0 . 23 0 . 11 3 . 48 × 10 − 2 F3 256 1 . 47 1 . 32 0 . 99 1 . 26 6 . 22 5 . 31 1024 1 . 14 1 . 05 1 . 01 1 . 26 5 . 67 4 . 16 F4 256 6 . 73 1 . 77 1 . 03 2 . 15 18 . 37 1 . 83 1024 0 . 28 0 . 14 8 . 16 × 10 − 2 0 . 17 2 . 12 0 . 23 6 Conclusion In this pap er, a definition of the discrete linear canon- ical transform (DLCT) based on hyperdifferential op- erator theory is prop osed. F or finite-length signals of a discrete v ariable, a unitary DLCT matrix is ob- tained so that the LCT-transformed version of the input signal can b e obtained by direct matrix m ul- tiplication. Given a vector holding the samples of a con tinuous-time signal, this DLCT matrix multiplies the v ector to obtain the approximate samples of the con tinuous-time LCT-transformed signal, similar to the DFT b eing used to approximate the contin uous- time F ourier transform. The adv an tage of a discrete transform is that it pro vides a basis for numerical computation. How- ev er, our exp ectations were more than that. The main goal of this work was to obtain a formulation of the discrete LCT based on self-consistent defini- tions of the discrete co ordinate m ultiplication and differen tiation op erators, that mirror the structure of their con tinuous counterparts. Care was taken to ensure that the discrete co ordinate multiplication and differen tiation op erators w ere strictly duals of eac h other, related through the DFT. The resulting DLCT matrix is totally compatible with the theory 10 of the discrete F ourier transform (DFT) and its dual and circulan t structure. Desirable prop erties of a dis- crete LCT definition suc h as unitarity , preserv ation of group structure, rev ersibility and approximation of the contin uous LCT were discussed b oth theoret- ically and numerically . One immediate p ossibilit y for future work is to explore the application of the metho d to alternative decomp ositions, such as those discussed in [31, 31, 38]. W e show ed in [38], that we could digitally com- pute the contin uous LCT to an accuracy limited by the uncertaint y relationship, with a fast algorithm. Ho wev er, this numerical computation method did not exhibit prop erties w e desire from a discrete definition. On the other hand, without a fast algorithm, applica- tion of the definition prop osed in the presen t pap er in volv es a matrix multiplication and thus has com- plexit y O ( N 2 ). The b est of both worlds would b e to find a fast algorithm for the definition prop osed in the present pap er. This would b e analogous to first defining the DFT and then deriving the FFT al- gorithm for its fast computation. How ever, suc h an algorithm is presently not av ailable and will require future work. In the mean time, fast computational metho ds as in [27, 30, 31, 38] can b e used in practical applications when sp eed is important. The compu- tational complexity of taking the DLCT of signals, whic h is a matrix multiplication with O ( N 2 ) com- plexit y , should not b e confused with the complexity of constructing the prop osed DLCT matrix, whic h has to b e done once for a particular LCT. The latter is discussed in App endix D. In the present pap er our emphasis was to define the DLCT in a manner that preserv es structural sim- ilarit y with the contin uous DLCT. The structure in question is ho w the LCT is defined in terms of co or- dinate m ultiplication and differentiation in terms of h yp erdifferen tial op erators, which we follow ed closely . Since everything rests on these tw o operators, their accuracy is what defines the accuracy of the metho d. W e chose the conceptually simplest first-order ap- pro ximations for these. Accuracy can b e increased either by increasing N , or by replacing these build- ing blocks with higher-order approximations. Thus, the h yp erdifferen tial form ulation provided here con- stitutes not only a theoretically pure approach to defining the DLCT, it serv es as a framew ork for high accuracy numerical computations. In conclusion, w e ha ve applied h yp erdifferential op- erator theory to the task of defining the discrete LCT in a manner that is fully consistent with the dual and circulan t structure of the DFT. Although sev eral def- initions for the DLCT hav e b een prop osed, a com- prehensiv e ev aluation of their relationships remains an imp ortan t sub ject for future work. W e believe our prop osed analytical approach can lead to further p ossible researc h directions in the theory of discrete transforms in general. App endix A Pro of of Unitarit y W e start with U giv en in Eq. 33. U is a real diagonal matrix, which implies it is Hermitian. The next step is to show D is also Hermitian. Starting from Eq. 34, w e can write D H = ( F − 1 UF ) H = F H U H ( F H ) H = F − 1 UF = D implying that D is also Hermitian. Now, we mo v e on to sho w that Q q , M M , and F a lc are unitary giv en U and D are Hermitian, by sho wing that their in verses and their Hermitians are equal. The inv erse of Q q is Q − 1 q = Q − q = exp  i 2 π q U 2 2  (35) while the Hermitian of Q q is Q H q = exp  i 2 π q ( U H ) 2 2  = exp  i 2 π q U 2 2  , (36) whic h are equal to eac h other. Similarly , one can follo w the same pro cedure for M M as follows: M − 1 M = M 1 / M = exp  − i 2 π ln (1 / M ) UD + DU 2  = exp  i 2 π ln ( M ) UD + DU 2  (37) and M H M = exp  i 2 π ln ( M ) ( UD + DU ) H 2  = exp  i 2 π ln ( M ) DU + UD 2  = M − 1 M . (38) 11 And, finally for F a lc w e can write: ( F a lc ) − 1 = F − a lc = exp  iaπ 2 U 2 + D 2 2  (39) and ( F a lc ) H = exp  iaπ 2 ( U 2 + D 2 ) H 2  = ( F a lc ) − 1 . (40) The first equalities in Eqs. 36, 38, and 40 can be sho wn b y considering p o wer expansion formula (Ap- p endix B). Thus we hav e prov en Theorem 2 and therefore Theorem 1. Justifications for the in terme- diate steps ab ov e will b e given in the App endix B. App endix B Some F undamen- tals of Op erator Theory Here w e provide further details regarding the deriv a- tions that app ear in Section 3 and App endix A. These deriv ations are mostly based on the follo wing elemen- tary definitions or results: (i) The integer p o wer of an operator is defined as its repeated application, e.g. A 3 = AAA . (ii) Therefore, an y p o wer of A commutes with itself, i.e. A n A = AA n . (iii) This leads to the fact that any p olynomial p ( A ) of A comm utes with A , i.e. p ( A ) A = A p ( A ). (iv) F unctions such as exp( A ) and sin( A ) can b e defined through p o w er series of exp( · ) and sin( · ), which are essentially lik e p olynomi- als, therefore these functions of A also comm ute with A . (v) Carrying this one step further, tw o different functions of A that can b e expressed as pow er series will also commute with each other, again as a conse- quence of (ii). (vi) The Hermitian of p ( A ), and thus also exp( A ) and sin( A ) can b e obtained by replacing A with its Hermitian inside the p o wer series. This follo ws from the fact that ( A n ) H = ( A H ) n . Eq. 31 follows directly from (iv) ab o ve. Eq. 32 fol- lo ws from the fact that the effect of U on a contin uous signal f ( u ) is to m ultiply it with u , and the fact that sin( U ) can b e written as a p o w er series of U . The steps in Eqs. 35 to 40 in the App endix A are most clearly established as follo ws. F or the first equalit y in Eq. 36, it follo ws from (vi) in the es- tablished facts ab o ve. With regards to Eq. 35, we observ e that Eqs. 9 and 10 show that the in verse of the chirp multiplication op erator is again a sim- ilar op erator but with negative parameter. Simi- lar observ ations can b e made for the other op er- ators by referring to their 2 × 2 matrices. Re- garding Eq. 35, this means that the inv erse of a c hirp m ultiplication op erator is of the same form but with negative parameter − q . So we need to show that exp( i 2 π q U 2 / 2) exp( − i 2 π q U 2 / 2) is equal to the iden tity . Here we can inv oke the Baker-Campbell- Hausdorff formula for matrices, [18, 24], which states that exp( A ) exp( B ) = exp( A + B + 1 / 2( AB − BA )) , (41) for t wo complex matrices A and B where both A and B commute with their commutator ( AB − BA ). In our case, A = − B , so that ( AB − BA ) = 0 . Therefore, the Baker-Campbell-Hausdorff formula’s condition is met since every matrix comm utes with the zero matrix. Finally , w e observe that the product on the left-hand side of the abov e identit y b ecomes equal to the exponential of the zero matrix and there- fore the identit y operator, pro ving the claim. Exactly the same argumen t applies for Eq. 37 and Eq. 39 since, although the exp onen ts are more complicated, in each case a minus sign is in tro duced to the exp o- nen t but otherwise the exponent remains the same. Therefore the exp onent of the original and the in v erse are merely negatives of eac h other and will commute, so that the pro duct of the original and inv erse matri- ces will b e the iden tity . The sinc( x ) = sin( π x ) / ( π x ) function has a p ow er series that is obtained by dividing the p o wer series of sin( π x ) b y ( π x ). F rom n umber (iv) of our elemen tary results, sinc( h D ) commutes with D , so b oth forms in Eq. 28 are the same. The same is true for Eq. 31. 12 App endix C Computation of the Matrix Exp o- nen tial Although it may b e view ed as an implementation de- tail, giv en that it lies at the heart of the prop osed metho d, it is worth clarifying how to compute the matrix exponential operation in Eq. 26. In practice, it is common to use MA TLAB ’s standard routines to compute matrix exponentials. Mathematically , the wa y in which matrix exp onen tials are obtained is through the w ell-known eigen decomp osition A = PDP − 1 (42) where D is a diagonal matrix that holds the eigen- v alues of A and P is the matrix holding the eigen- v ectors. Then, exp( A ) = P exp( D ) P − 1 where the exp() that operates on D is now simply an element- wise exp onen tiation op eration. When A has a full set of eigenv alues, this pro cedure works without any complication. Given Eqs. 33 and 34, and the unitar- it y of the DFT matrix F , the matrices U and D are ensured to hav e a full set of eigenv alues and eigen- v ectors, so there is no mathematical complication in using matrix exp onentials. App endix D Computational Cost of Con- structing the Prop osed DLCT Matrix Giv en a sp ecified precision (i.e., num b er of bits used in computations is fixed), to find the complexity of generating the matrix C L as a function of N , we first find the complexit y of computing the matrices U and D . The matrix U is generated using Eq. 33. This pro cess requires ev aluation of the sine function at N p oin ts and N multiplications b y the constant √ N /π . Since w e assume a fixed precision, w e can take the ev aluation of the sine function at a p oin t to b e of complexit y O (1). The complexit y of computing U is th us O ( N ). Secondly , to compute D using Eq. 34, w e need to compute the matrix F and F − 1 , b oth of which can b e written in terms of W N . In gener- ating F , we compute W N only once and compute its ( mn )’th pow er for the ( mn )’th en try . Computing the ( mn )’th en try for the matrices F and F − 1 requires t wo m ultiplications and one exp onen tiation, whic h are each taken to b e O (1). It follo ws that computing F and F − 1 eac h takes O ( N 2 ) computations. Finally , m ultiplying F − 1 with U is O ( N 2 ) since U is diagonal whereas multiplying F − 1 U with F is O ( N 2 log N ) (by using fast F ourier transform (FFT) algorithm and by noting that neither matrices are diagonal), resulting in an ov erall complexity of O ( N 2 log N ) for D . W e can no w mov e on to the complexities of com- puting the matrices Q q , M M , F a lc based on Eqs. 23, 24, and 25. Note that in Eqs. 23, 24, and 25, the scalar constants can be tak en outside the exp() func- tion, b e computed separately and then b e multiplied with the resulting matrix exp onen tials. This does not ha ve an effect on the computational complexity with resp ect to N . • Complexity of Q q : T aking the square of U is of complexit y O ( N ) since U is a diagonal matrix. W e can compute the matrix exponential of U 2 simply b y taking the exp onen tial of eac h diago- nal elemen t b ecause U 2 is also a diagonal ma- trix. This amounts to an ov erall computational complexit y of O ( N ). • Complexity of M M : One can compute b oth UD and DU in O ( N 2 ) time because U is a diag- onal matrix. Ho wev er, generating D increases the time to compute the argument of the exp() to O ( N 2 log N ). F urthermore, computing ma- trix exp onen tials as describ ed in App endix C is of complexit y O ( N 3 ). As a result, the o verall complexit y is O ( N 3 ). • Complexity of F a lc : This is the same as the com- plexit y of M M since it in volv es computing the matrix exp onen tial of a non-diagonal matrix. In conclusion, the o verall complexity for computing the matrix C L is O ( N 3 ). 13 References [1] Sparse representation of tw o- and three- dimensional images with fractional Fourier, Hartley , linear canonical, and Haar wa velet transforms. Exp ert Systems with Applic ations , 77:247 – 255, 2017. [2] S. Ab e and J. T. Sheridan. Generalization of the fractional Fourier transformation to an arbi- trary linear lossless transformation an op erator approac h. Journal of Physics A: Mathematic al and Gener al , 27(12):4179–4187, 1994. [3] S. Ab e and J. T. Sheridan. Optical operations on w av efunctions as the abelian subgroups of the sp ecial affine Fourier transformation. Opt. L ett. , 19:1801–1803, 1994. [4] T. Aliev a and M. J. Bastiaans. Prop erties of the canonical in tegral transformation. J. Opt. So c. A m. A , 24:3658–3665, 2007. [5] N. M. A takishiyev, L. E. Vicent, and K. B. W olf. Con tinuous vs. discrete fractional Fourier trans- forms. Journal of Computational and Applie d Mathematics , 107(1):73 – 95, 1999. [6] Laurence Barker, Cagatay Candan, T ugrul Hakioglu, M Alp er Kutay , and Haldun M Oza- ktas. The discrete harmonic oscillator, Harp er’s equation, and the discrete fractional Fourier transform. Journal of Physics A: Mathematic al and Gener al , 33(11):2209, 2000. [7] B. Barshan, M. A. Kutay , and H. M. Ozaktas. Optimal filtering with linear canonical transfor- mations. Optics Communic ations , 135(1-3):32 – 36, 1997. [8] M. J. Bastiaans. The Wigner distribution func- tion applied to optical signals and systems. Op- tics Communic ations , 25(1):26 – 30, 1978. [9] M. J. Bastiaans. Wigner distribution function and its application to first-order optics. J. Opt. So c. Am. , 69:1710–1716, 1979. [10] M. J. Bastiaans and T. Aliev a. Classification of lossless first-order optical systems and the linear canonical transformation. J. Opt. So c. Am. A , 24:1053–1062, 2007. [11] L. M. Bernardo. T alb ot self-imaging in fractional Fourier planes of real and complex orders. Optics Communic ations , 140:195–198, 1997. [12] L. M. Bernardo and O. D. D. Soares. Optical fractional Fourier transforms with complex or- ders. Appl. Opt. , 35(17):3163–3166, 1996. [13] R. G. Campos and J. Figueroa. A fast algorithm for the linear canonical transform. Signal Pr o- c essing , 91(6):1444–1447, 2011. [14] C. Candan, M. A. Kutay , and H. M. Ozaktas. The discrete fractional Fourier transform. Signal Pr o c essing, IEEE T r ansactions on , 48(5):1329 – 1337, 2000. [15] X. Chen, J. Guan, Y. Huang, N. Liu, and Y. He. Radon-linear canonical am bi- guit y function-based detection and estimation metho d for marine target with micromotion. IEEE T r ans. Ge oscienc e and R emote Sens. , 53(4):2225–2240, 2015. [16] X. Chen, J. Guan, N. Liu, W. Zhou, and Y. He. Detection of a low observ able sea-surface target with micromotion via the radon-linear canonical transform. IEEE Ge osci. R emote Sensing L ett. , 11(7):1225–1229, 2014. [17] Stuart Clary and Dale H. Mugler. Shifted Fourier matrices and their tridiagonal comm u- tors. SIAM Journal on Matrix Analysis and Ap- plic ations , 24(3):809–821, 2003. [18] C. Cohen-T annoudji, B. Diu, and F. Laloe. Quantum Me chanics . New Y ork: Wiley , 1979. [19] B. Da vies. Inte gr al T r ansforms and Their Ap- plic ations . Springer, New Y ork, 1978. [20] T. Erseghe, P . Kraniausk as, and G. Carioraro. Unified fractional Fourier transform and sam- pling theorem. Signal Pr o c essing, IEEE T r ans- actions on , 47(12):3419 –3423, dec. 1999. 14 [21] Qiang F eng and Bing-Zhao Li. Con volution and correlation theorems for the t wo-dimensional lin- ear canonical transform and its applications. IET Signal Pr o c essing , 10:125–132, 2016. [22] F.Alb erto Grun baum. The eigen vectors of the discrete Fourier transform: A version of the Her- mite functions. Journal of Mathematic al Analy- sis and Applic ations , 88(2):355 – 363, 1982. [23] T alha Cihad Gulcu and Haldun M. Ozaktas. Choice of quantization interv al for finite-energy fields. IEEE T r ansactions on Signal Pr o c essing , 66:2470–2479, 2018. [24] B. C. Hall. Lie Gr oups, Lie Algebr as, and R ep- r esentations (Chapter 3: The BakerCampb el l- Hausdorff F ormula) . New Y ork: Springer, 2003. [25] J. J. Healy , M. A. Kutay , H. M. Ozaktas, and J. T. Sheridan eds. Line ar Canonic al T r ans- forms: The ory and Applic ations . Springer New Y ork, New Y ork, NY, 2016. [26] J. J. Healy and J. T. Sheridan. Sampling and discretization of the linear canonical transform. Signal Pr o c ess. , 89(4):641–648, 2009. [27] J. J. Healy and J. T. Sheridan. F ast linear canon- ical transforms. J. Opt. So c. Am. A , 27(1):21– 30, 2010. [28] John J. Healy and John T. Sheridan. Reev alu- ation of the direct metho d of calculating F res- nel and other linear canonical transforms. Opt. L ett. , 35(7):947–949, Apr 2010. [29] E. Hec ht. Optics, 4th Ed. Addison W esley , 2001. [30] B. M. Hennelly and J. T. Sheridan. F ast n umeri- cal algorithm for the linear canonical transform. J. Opt. So c. A m. A , 22:928–937, 2005. [31] B. M. Hennelly and J. T. Sheridan. General- izing, optimizing, and inv enting n umerical algo- rithms for the fractional Fourier, F resnel, and linear canonical transforms. J. Opt. So c. Am. A , 22:917–927, 2005. [32] J. Hua, L. Liu, and G. Li. Extended frac- tional Fourier transforms. J. Opt. So c. Am. A , 14(12):3316–3322, 1997. [33] D. F. V. James and G. S. Agarwal. The gen- eralized F resnel transform and its application to optics. Optics Communic ations , 126(4-6):207 – 212, 1996. [34] C. Jung and H. Kruger. Represen tation of quan- tum mechanical wa vefunctions b y complex v al- ued extensions of classical canonical transfor- mation generators. J. Phys. A: Math. Gen. , 15:3509–3523, 1982. [35] A. W. Knapp. R epr esentation the ory of semisim- ple gr oups: An overview b ase d on examples . Princeton Universit y Press, 2001. [36] A. Ko¸ c, B. Bartan, and H. M. Ozaktas. Discrete scaling based on op erator theory . arXiv pr eprint arXiv:1805.03500 , 2018. [37] A. Ko¸ c, F. S. Oktem, H. M. Ozaktas, and M. A. Kuta y . Line ar Canonic al T r ansforms: The ory and Applic ations , chapter F ast Algorithms for Digital Computation of Linear Canonical T rans- forms, pages 293–327. Springer New Y ork, New Y ork, NY, 2016. [38] A. Ko¸ c, H. M. Ozaktas, C. Candan, and M. A. Kuta y . Digital computation of linear canoni- cal transforms. IEEE T r ans. Signal Pr o c ess. , 56(6):2383–2394, 2008. [39] A. Ko¸ c, H. M. Ozaktas, and L. Hesselink. F ast and accurate algorithm for the computation of complex linear canonical transforms. J. Opt. So c. Am. A , 27(9):1896–1908, 2010. [40] A. Ko¸ c, H. M. Ozaktas, and L. Hesselink. F ast and accurate computation of tw o-dimensional non-separable quadratic-phase integrals. J. Opt. So c. Am. A , 27(6):1288–1302, 2010. [41] B.Z. Li and Y.P . Shi. Image watermarking in the linear canonical transform domain. Mathe- matic al Pr oblems in Engine ering , 2014. 15 [42] M. Moshinsky . Canonical transformations and quan tum mechanics. SIAM J. Appl. Math. , 25(2):193–212, 1973. [43] M. Moshinsky and C. Quesne. Linear canon- ical transformations and their unitary repre- sen tations. Journal of Mathematic al Physics , 12(8):1772–1780, 1971. [44] D. H. Mugler. The cen tered discrete Fourier transform and a parallel implemen tation of the fft. In 2011 IEEE International Confer- enc e on A c oustics, Sp e e ch and Signal Pr o c essing (ICASSP) , pages 1725–1728, Ma y 2011. [45] M. Nazarathy and J. Shamir. F ourier optics de- scrib ed by op erator algebra. J. Opt. So c. Am. , 70:150–159, 1980. [46] F. S. Oktem and H. M. Ozaktas. Exact relation b et w een contin uous and discrete linear canoni- cal transforms. Signal Pr o c essing L etters, IEEE , 16(8):727 –730, 2009. [47] H. M. Ozaktas, A. Ko¸ c, I. Sari, and M. A. Ku- ta y . Efficien t computation of quadratic-phase in- tegrals in optics. Opt. L ett. , 31:35–37, 2006. [48] H. M. Ozaktas, Z. Zalevsky , and M. A. Kutay . The F r actional Fourier T r ansform with Appli- c ations in Optics and Signal Pr o c essing . New Y ork: Wiley , 2001. [49] Haldun M. Ozaktas, M. Alper Kutay , and Ca- gata y Candan. T r ansforms and Applic ations Handb o ok , chapter F ractional Fourier T rans- form, pages 14–1–14–28. CRC Press, Bo ca Ra- ton, New Y ork, NY, 2010. [50] A. Ozcelikk ale and H. M. Ozaktas. Represen ta- tion of optical fields using finite num b ers of bits. Optics letters , 37:2193–2195, 2012. [51] Aya Ozcelikk ale and Haldun M. Ozaktas. Op- timal representation of non-stationary random fields with finite num b ers of samples: A linear MMSE framework. Digital Signal Pr o c essing , 23(5):1602 – 1609, 2013. [52] C. Palma and V. Bagini. Extension of the F res- nel transform to abcd systems. J. Opt. So c. A m. A , 14(8):1774–1779, 1997. [53] S. C. P ei and J. J. Ding. Closed-form discrete fractional and affine Fourier transforms. IEEE T r ans. Signal Pr o c ess. , 48:1338–1353, 2000. [54] S. C. P ei and J. J. Ding. Eigenfunction of linear canonical transform. IEEE T r ans. Signal Pr o- c ess. , 50:11–26, 2002. [55] S. C. Pei and S. Huang. F ast discrete linear canonical transform based on cm-cc-cm decom- p osition and fft. IEEE T r ans. Signal Pr o c ess. , 64:855 – 866, 2016. [56] S. C. Pei and Y. C. Lai. Discrete linear canonical transforms based on dilated Hermite functions. J. Opt. So c. A m. A , 28(8):1695 –1708, 2011. [57] S. C. Pei and M. H. Y eh. Impro ved dis- crete fractional Fourier transform. Opt. L ett. , 22(14):1047–1049, 1997. [58] S. C. Pei, M. H. Y eh, and T. L. Luo. F rac- tional Fourier series expansion for finite signals and dual extension to discrete-time fractional Fourier transform. Signal Pr o c essing, IEEE T r ansactions on , 47(10):2883 –2888, 1999. [59] S. C. P ei, M. H. Y eh, and C. C. Tseng. Dis- crete fractional Fourier transform based on or- thogonal pro jections. Signal Pr o c essing, IEEE T r ansactions on , 47(5):1335 –1348, 1999. [60] Min Qi, Bing-Zhao Li, and Huafei Sun. Image w atermarking using p olar harmonic transform with parameters in SL(2,R). Signal Pr o c essing: Image Communic ation , 31:161 – 173, 2015. [61] W. Qiu, B.Z. Li, and X.W. Li. Sp eec h recov ery based on the linear canonical transform. Sp e e ch Communic ation , 55(1):40–50, 2013. [62] J. Ro drigo, T. Aliev a, and M. Luisa Calvo. Op- tical system design for orthosymplectic transfor- mations in phase space. J. Opt. So c. A m. A , 23:2494–2500, 2006. 16 [63] K. K. Sharma. F ractional Laplace transform. Signal, Image And Vide o Pr o c ess. , 4(3):377–379, 2009. [64] C. C. Shih. Optical interpretation of a complex-order Fourier transform. Opt. L ett. , 20(10):1178–1180, 1995. [65] A. E. Siegman. L asers . Mill V alley , California: Univ ersity Science Bo oks, 1986. [66] R. Simon and K. B. W olf. Structure of the set of paraxial optical systems. J. Opt. So c. Am. A , 17(2):342–355, 2000. [67] N. Singh and A. Sinha. Chaos based m ultiple image encryption using m ultiple canonical trans- forms. Optics and L aser T e chnolo gy , 42:724–731, 2010. [68] A. Stern. Wh y is the linear canonical transform so little known? IP Conf. Pr o c. , pages 225–234, 2006. [69] A. T orre. Linear and radial canonical transforms of fractional order. J. Compt. and Appl. Math. , 153:477–486, 2003. [70] J. G. V argas-Rubio and B. Santhanam. On the m ultiangle centered discrete fractional Fourier transform. IEEE Signal Pr o c ess. L ett. , 12(4):273–276, 2005. [71] C. W ang and B. Lu. Implemen tation of complex- order Fourier transforms in complex ab cd op- tical systems. Optics Communic ations , 203(1- 2):61 – 66, 2002. [72] Deyun W ei, Ruikui W ang, and Y uan-Min Li. Random discrete linear canonical transform. J. Opt. So c. Am. A , 33(12):2470–2476, Dec 2016. [73] K. B. W olf. Canonical transformations i. complex linear transforms. J. Math. Phys. , 15(8):1295–1301, 1974. [74] K. B. W olf. On self-recipro cal functions under a class of integral transforms. J. Math. Phys. , 18(5):1046–1051, 1977. [75] K. B. W olf. Inte gr al T r ansforms in Scienc e and Engine ering (Chapter 9: Construction and pr op erties of c anonic al tr ansforms) . New Y ork: Plen um Press, 1979. [76] K. B. W olf. Finite systems, fractional Fourier transforms and their finite phase spaces. Cze choslovak Journal of Physics , 55:1527–1534, 2005. [77] K. B. W olf. Line ar Canonic al T r ansforms: The ory and Applic ations , chapter Dev elopment of Linear Canonical T ransforms: A Historical Sk etch, pages 3–28. Springer New Y ork, New Y ork, NY, 2016. [78] K. B. W olf and G. Kr¨ otzsc h. Geometry and dy- namics in the fractional discrete Fourier trans- form. J. Opt. So c. Am. A , 24(3):651–658, 2007. [79] M. H. Y eh. Angular decomp ositions for the dis- crete fractional signal transforms. Signal Pr o- c essing , 85(3):537 – 547, 2005. [80] I. S. Y etik, M. A. Kuta y , H. Ozaktas, and H. M. Ozaktas. Contin uous and discrete fractional Fourier domain decomp osition. v olume 1, pages 93 –96 vol.1, 2000. [81] K. Y osida. Op er ational Calculus: A The ory of Hyp erfunctions . Springer, New Y ork, USA, 1984. [82] A. I. Zay ed and A. G. Garca. New sampling for- m ulae for the fractional Fourier transform. Sig- nal Pr o c essing , 77(1):111 – 114, 1999. [83] F eng Zhang, Ran T ao, and Y ue W ang. Discrete linear canonical transform computation by adap- tiv e metho d. Opt. Expr ess , 21(15):18138–18151, Jul 2013. [84] Juan Zhao, Ran T ao, and Y ue W ang. Sampling rate conv ersion for linear canonical transform. Signal Pr o c essing , 88(11):2825 – 2832, 2008. [85] Liang Zhao, John J. Healy , and John T. Sheri- dan. Unitary discrete linear canonical transform: analysis and application. Appl. Opt. , 52(7):C30– C36, Mar 2013. 17 -4 -3 -2 -1 0 1 2 3 4 -0.2 0 0.2 0.4 0.6 0.8 1 Proposed DLCT Reference (a) Real part of T1 of F1 -4 -3 -2 -1 0 1 2 3 4 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 Proposed DLCT Reference (b) Imaginary part of T1 of F1 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -2 -1 0 1 2 Proposed DLCT Reference (c) Real part of T2 of F2 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -2 -1 0 1 2 Proposed DLCT Reference (d) Imaginary part of T2 of F2 -6 -4 -2 0 2 4 6 -0.5 0 0.5 1 1.5 Proposed DLCT Reference (e) Real part of T3 of F3 -6 -4 -2 0 2 4 6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 Proposed DLCT Reference (f ) Imaginary part of T3 of F3 -4 -3 -2 -1 0 1 2 3 4 -0.5 0 0.5 Proposed DLCT Reference (g) Real part of T4 of F4 -4 -3 -2 -1 0 1 2 3 4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 Proposed DLCT Reference (h) Imaginary part of T4 of F4 Figure 1: Comparison of the proposed DLCT of functions with the reference. 18