Research Progress on Discretization of Linear Canonical Transform

Yannan Sun,Bingzhao Li,Ran Tao

Abstract:Linear canonical transformation (LCT) is a generalization of the Fourier transform and fractional Fourier transform.The recent research has shown that the LCT is widely used in signal processing and applied mathematics,and the discretization of the LCT becomes vital for the applications of LCT.Based on the development of discretization LCT,a review of important research progress and current situation is presented,which can help researchers to further understand the discretization of LCT and can promote its engineering application.Meanwhile,the connection among different discretization algorithms and the future research are given.

Keywords:linear canonical transform(LCT);discrete linear canonical transform;sampling;Wigner-Ville distribution;fast algorithm

1 Introduction

Linear canonical transform (LCT) is a class of linear integral transform with three parameters,which is the generalization of the Fourier transform (FT),the fractional Fourier transform(FRFT).It was first introduced in the 1970s by Collins and Moshinsky Quesne [1,2].The detail introduction of the LCT is given in [3].The definition of the LCT is given as [4]

Eq.(1) is the LCT in the sense of angular frequency.Through variable substitution,the other definition of LCT can be obtained,and the transformation kernel is as

In general,the definition of LCT with three parameters is given by takingα=d/b,β=1/b,γ=a/b

These definitions forms are different,but they can transform into each other.The LCT has some important and useful properties:reversibility,additively and unitary.The additive property implies that the LCT can be decomposed into a series of multiple LCTs.The unitary property is that ILCT is the conjugate transpose of LCT with parameter matrix,namelyKA−1(t,u)=

In recent years,the LCT theories have been widely concerned [5−42].Scholars at home and abroad have done a lot of research in discretization method and fast algorithm [13−21],filtering[12],sampling [22−27],convolution theory [28−30],and so on [31−41].From the current research results,the basic theory of LCT is constantly improving.Compared with Fourier transform and fractional Fourier transform (FRFT),there are many aspects to be further studied,such as the discretization and fast algorithm of LCT[13,41,42].In 2000,Pei obtained the closed form definition of discrete linear canonical transform(DLCT) for the first time by sampling the time domain and LCT domain [13].After that,many discretization methods and fast algorithms of LCT have been proposed.For example,in 2005,Sheridan proposed a basis decomposition fast algorithm based on the idea of fast Fourier transform (FFT) [41].Since then,the algorithm based on LCT parameter matrix factorization has been widely concerned [42].For different types of signals and different factorization forms,various discretization and fast calculation methods are proposed.At present,some review articles have been given on the discretization methods of FRFT [43,44].However,there is still no a review on LCT discretization methods and fast algorithms.This paper reviews the important research progress and current situation of LCT discretization algorithm [45−56].

This paper is organized as follows.The different definitions of direct discrete LCT are summarized in Section 2.Section 3,the fast DLCT algorithm based on different decomposition bases is studied.In Section 4,a fast LCT algorithm based on LCT parameter matrix factorization is discussed.In Section 5,other discretization and fast implementation methods of DLCT are introduced.Conclusions and future research are given in Section 6.

2 Direct Discrete LCT

There are two kinds of direct discrete LCT:the one is to discretize the variablestandudirectly,and change the integral into summation.The other is to obtain the discrete time LCT(DTLCT),and then obtain DLCT.These two discretization methods are different,but they are equivalent under certain conditions.Pei and Sheridan have proposed the following definition of DLCT by using these two methods,namely,Pei method and Sheridan method.

2.1 Pei Method

Pei used the first method to obtain the following DLCT in 2000 [13]

where ∆tand ∆uis sampling interval in time domain and LCT domain,m=−M/2,−M/2+1,···,M/2−1,NandMare sampling points in time domain and LCT domain,respectively.In general,this direct discretization method didn’t satisfy the unitary property.However,Eq.(4) is unitary for.By normalizing the transformation matrix,the DLCT with unitary and reversible is obtained

The corresponding inverse discrete linear canonical transformation (IDLCT) is

Eq.(5) is not true forb=0.Whenb=0,dis integer,the variableuis discretized directly,

In order to approach continuous LCT,the sampling interval of this kind of DLCT should satisfy the sampling theorem

Asdis not integer,we will use LCT parameter matrix decomposition to get the discrete form,which will be analyzed and discussed in detail in Section 4.2.After this kind of algorithm was proposed,Sheridan proposed DLCT by using the second method.

2.2 Sheridan Method

Hennelly [41]proposed the DTLCT as

and discretization of variableu

where−N/2 ≤m≤N/2−1,the discrete transformation kernel is

ForM=N,this two DLCT in [13]and [41]are essentially equivalent,which required the sampling interval to satisfy ∆t∆u=2π|b|/Nand∆t∆u=1/(N|β|) respectively.In [51],it is shown that DLCT has unitary property for∆u=q2π|b|/(N∆t),qandNare prime.For DLCT with unitary property,this property is very important in engineering applications,such as phase modulation.Therefore,the conclusion of [51]will be very important for the DLCT.

The definition of DLCT in [41]is constrained by the extent of input,output range,and output sampling rate.To overcome this limitation,the following three improved DLCT are proposed in [52],namely,scale-preserving DLCT,scale variation DLCT,arbitrarily scaled DLCT.Next,the existing fast DLCT was introduced.

3 Basis Decomposition DLCT

In recent years,some fast DLCT algorithms based on basis decomposition have been proposed using the time-shift and frequency-shift properties,chirp periodic properties and symmetry properties of discrete transform kernel and so on [53−56].

3.1 Frequency Decimation Algorithm

The radix-2 frequency decimation algorithm was proposed in [54]

Then,Eq.(9) is simplified as the following two cases based on the parity ofm

In order to make the length of decomposition more flexible,a general frequency domain decimation algorithm was proposed in [56].This algorithm only required that the length of the input data is a composite number,and the length of the split base is a factor ofN.It is more flexible.In addition,the discrete data is input sequentially in the decomposition process,but the output sequence is extracted at a certain interval.For the opposite case,the time domain d ecimation algorithms are proposed.

3.2 Time Decimation Algorithm

A radix-2 time extraction algorithm was proposed in [57].Firstly,the input discrete sequence was decomposed into two parts,one part contains even discrete pointsx(2n∆t) and the other contains odd discrete pointsx((2n+1)∆t),

where−N/2 ≤m≤−1.Therefore,twoN/2-point DLCT could be used to obtain a DLCT with a length ofN.Similarly,the sequence ofN/2 points can be repeatedly decomposed and iterated until it is decomposed into 2-point DLCT.

In order to improve the flexibility of the algorithm,[55]further proposed a general basis time extraction algorithm,where the length of input data was required to be composite.This is similar to the frequency extraction algorithm,and the length of the basis can be selected according to the actual needs.From the current research situation,many literatures have proposed different basis decomposition algorithms [58],but there was no optimal selection method,which needs further research.In addition,in many references,using the additive property of LCT and the existing algorithm of LCT special case,the operator decomposition algorithm was given.

4 Operator Decomposition Algorithm

The LCT operator is decomposed into the product of special operators with fast algorithms,such as Fourier operator,fractional Fourier operator,Fresnel operator,chirp product,scaling operator,etc.,and the fast algorithm of LCT is obtained based on the ones of these special operators.This type of algorithm is called operator decomposition algorithm [59−63].It focused on the sufficiency of sampling in the process of decomposition to ensure that the discrete points can recover the continuous transformation well.The determination of sampling rate and sampling points is the key of the algorithm.The calculation of signal wide bandwidth product (SBP) is given in [42,63].In the following,the existing fast LCT algorithms will be summarized in detail.

4.1 Matrix Factorization Algorithm for Time Frequency with Any Shape

The most widely FFT can be used to calculate fractional Fourier transform [64],FST.Therefore,the LCT can be decomposed into the transform including FT to study the discrete algorithm [65].Combined with the existing discretization algorithm of FST and FRFT,the LCT operator was decomposed into the combination of FT,ST,CM,FST and FRFT operators in [42],and the corresponding fast algorithm was proposed.Based on the decomposition of FST and FRFT,[42]proposed the decomposition of LCT

wheremcos(pπ/2)=a;mλqsin(pπ/2)=b.Selectingfto makea/λf −sin(pπ/2)/mqλ=c.

According to [42],LCT can be decomposed into ten forms.Different decompositions had different interpolation or extraction multiples,and the corresponding computational complexity was different.When the time-frequency region is rectangular,Ozaktas proposed two matrix factorization DLCT algorithms for different transformation parameters in [63].

4.2 Matrix Factorization Algorithm for Rectangular Time Frequency Domain

In [63],Ozaktas normalized the signal before discretization,which limited the time and frequency range of the analyzed signal to the interval [−W/2,W/2],and the number of sampling points isN0=W2.Based on the parameter values,[63]proposed two different decomposition of LCT.The Method I was to decompose the LCT into a combination of scale transform,Fourier transform and chirp product,namely

The first matrix from right to left is scale transformation,which only changed the sampling interval of time domain and frequency domain.After this transformation,the time width was|β|W,and frequency width wasW/|β|,andN1=N0.The second matrix was CM,the time width and frequency width of the signal after this transformation wereW/|β|+W|γ/β|and|β|W,andN2=(1+|γ|)N0.The third matrix corresponded to FT,which did not require an increase sampling points.After the last CM,the SBP of the signal wasN3=(1+|γ|+|α|(1+|γ|)2/|β|2)W2.IfN3≤2N0,then interpolation and decimation were not needed.IfN3>2N0,the following inequality was established

wherekwas the smallest integer that made the inequality hold.Therefore,in the discrete process,the result of the previous transformation was used as the input of the new transformation,and it was necessary to judge whether the number of sampling points was sufficient.

We have also added the computational of each stage to obtain the overall computational complexity of the decomposition as

whereI(x,y) stands for the cost to interpolatexsamples by a factor ofyto obtainxysamples andS(N) stands for the cost of the Scaling operation onNsamples.

From the above decomposition,it can be concluded that only CM changes sampling points.Therefore,so as to reduce interpolation or extraction multiple,CM should be included as little as possible in the decomposition,and appear as late as possible in the order of matrix decomposition.

The second method is based on the following decomposition involving the FRFT,Scaling,and Chirp multiplication:

It was a special case of the widely known Iwasawa decomposition.When the algorithm was implemented,FRFT can be calculated by using the existing discrete algorithm of FRFT [66],and then the result can be used as the input to perform the scaling operation corresponding to the second matrix,and then the transformation result can be used as the input of the next transformation.Because the next transform is chirp product,it will change the signal bandwidth product (SBP) of the signal.To ensure the sufficiency of sampling,it was necessary to judge whether it was necessary to interpolate or extract according to the SBP.

In practical application,the appropriate decomposition form can be selected according to the different parameters to achieve fast calculation.Compared with the algorithm in [42]and[63],the former which can be regarded as the theoretical basis of the later focused on decomposition theory without considering the computational complexity,and the algorithm did not have any restriction on the shape of signal timefrequency region.However,[63]normalized the signal firstly,and considered the influence of the change of signal SBP,and analyzed the computational complexity of the algorithm.Although the interpolation times and computational complexity of these decomposition methods were different,the discrete transformation results can well approximate the case of continuous transformation.

In [15],Pei proposed DLCT by combining FRFT,scaling transformation and chirp product’s extensible Hermite function as characteristic function.Its decomposition method was the same as the method II in [58],but it did not need over sampling.This kind of discrete algorithm has also some good properties,such as unitary property,reversible property and additive property.It is the discrete algorithm that satisfies the most LCT properties.However,this algorithm has no closed form and the computational complexity wasO(N3).Therefore,it is very limited in practical application.Although the fast DLCT algorithm in [42,63]can well approximate the continuous LCT.The algorithm was generally irreversible,and did not satisfy the additive propsampling.A fast LCT algorithm for pure discrete date was presented.The DLCT ofx(n) was defined as erty,and needed interpolation or extraction in the process of calculation.The algorithm in [15]overcomes the shortcomings of [42,63],but it had no closed form and fast computational efficiency.In order to further overcome these shortage,[16]proposed a DLCT by CM-CC-CM decomposition [67,68],which independent of sampling interval and does not need over-

whereA′=[a′,b′;c′,d′],a′d′−b′c′=1.According to [42,63],CC can be further decomposed as following

Eq.(15) can be rewritten as

whereb′≠0

FandF −1represented discrete Fourier transform and inverse discrete Fourier transform respectively.Forb′=0,it can be transformed into the case ofb′≠0 by matrix decomposition.This kind of algorithm consisted of two FFT algorithm and three Chirp products,the total calculation was 3N+Nlog2(N) complex number times.

When matrixA′is transformed intoA′−1,the above algorithm was transformed into inverse DLCT

The relation between transformation parametersa,b,c,danda′,b′,c′,d′was

where ∆twas sampling interval.Fora,b,c,dwere equal toa′,b′,c′,d′respectively.

The DLCT was independent of sampling interval and had reversible and additive properties.When the transformation parameter matrix was taken as a special case,a fast algorithm for the special transformation can be obtained.

4.3 Matrix Factorization DLCT Algorithm Based on Sampling Theorem in LCT Domain

The fast DLCT in Sections 4.1 and 4.2 required the signal was band limited in time domain and traditional frequency domain,and the sampling should meet the Nyquist Shannon.However,the unitary DLCT proposed in [13]required the signal to be band limited in time domain and LCT domain,and the sampling interval was determined by the sampling theorem of LCT.The DLCT in [13]can also be implemented by matrix decomposition.This algorithm did not require dimensional normalization,but the algorithm did not hold forb=0.In this case,dis a non-integer,the LCT can be decomposed as

The corresponding DLCT was

where ∆u=|a|∆tN/M.This kind of algorithm did not research the exact relationship between LCT and DLCT,and how to determine the number of sampling points.

In [45],Oktem and Ozaktas gave the exact relationship between them as following

which were the periodically replicated functions of a continuous-time signalx(t) and its LCT;N=BW/|b|is the number of sampling,WandBwere arbitrary,were sampling interval.It was also worth noting that the functions (18) and (19) were Chirp-periodic.Eq.(17)was the generalization of the exact relation between continuous and discrete FTs.If a signal was aperiodic,most of the energy was concentrated in time domain [−W/2,W/2]and LCT domain [−B/2,B/2],the DLCT of a signal was

Ifx(t) andXA(u) were 0 in given interval outside,then Eq.(20) was exactly equal.As we all know,it was impossible for a signal to be bandlimited in both time domain and transform domain.Therefore,only Chirp periodic signals can have exact equality between DLCT and LCT.

5 Other DLCT

In addition to the above DLCT,there were other types of algorithms.Sun and Li has proposed fast LCT for local spectra with flexible resolution ability [69−71].In [72],Campos et al.proposed a fast DLCT algorithm by using the same decomposition method as [13],which independent of sampling interval and has unitary property.In [72,73],the approximate zeros of Hermite polynomials were used as the sampling points of transformation to realize the fast algorithm of DLCT.This algorithm can be well applied to image edge detection.In [74],Li proposed the characteristic discrete algorithm based on the characteristic function of LCT for|b|≤2,and gave the definition spectral decomposition which similar to FRFT[75].By discretizing the definition of spectral decomposition,the corresponding DLCT was obtained.This algorithm was simple to implement and easy to understand,but it had high computational complexity.In [76],Wei proposed the random discrete linear canonical transform by randomizing the kernel transform matrix of the discrete linear canonical.It has a greater degree of randomness because of the randomization in terms of both eigenvectors and eigenvalues.This method can be well applied to image encryption and information security,but it has the computational complexityN2.In [77],the DLCT using the adaptive least-mean-square algorithm was presented which had the characteristics of parallel computing and could be implemented effectively by VLSI.In recent years,some researchers have studied the complex linear canonical transformation [78,79].In [79],the algorithm of CM-CC-CC decomposition was proposed by using the method of parameter matrix decomposition.Unlike the previous decomposition [16,63],these matrix elements were complex numbers.In[80],a fast algorithm independent of sampling interval and sampling points was proposed by using FFT algorithm,which was widely used at present.In [81],Pei proposed four kinds of fast algorithms for complex linear canonical transform with high accuracy based on Gabor transform,and studied the discretization method of inverse Bargmann transform.

Although there were many LCT discretization methods,these algorithms have their own limitations.These algorithms do not satisfy the additive property except the algorithm in [15,16].For different types of signals,the fast implementation methods were selected.Therefore,it was meaningless to compare these algorithms.At present,there was no universal algorithm like FFT and FRFT,so there were still many aspects in the research of LCT discretization algorithm and fast implementation method.

6 Conclusion and Future Research

In this paper,we present the discretization of the LCT developed in the last nearly two decades,as well as their applications.The discretization included discrete-time LCT,LCS,discrete LCT.In order to achieve the digital computation of one continuous transform,we need an accurate and efficient discretization of this transform which is consistent with its continuous counterpart.After the LCT has been derived,many researchers have tried to derive their discrete counterpart.Because of the fast oscillation of the LCT kernel,the discretization of the LCT is more difficult than the conventional FT.Until now,no unified definition has been obtained.Currently the DLCT based on the Operator decomposition method has attracted more and more attention.

There are several works for future research:1) acquire fast LCT with satisfying more properties of continuous LCT;2) derive nonuniform fast LCT;3) establish the fast of the short-time LCT.We are confident that the above problems will be settled soon.