Fast Discrete Linear Canonical Transform Based on CM-CC-CM Decomposition and FFT

# Fast Discrete Linear Canonical Transform Based on CM-CC-CM Decomposition and FFT

## Abstract

In this paper, a discrete LCT (DLCT) irrelevant to the sampling periods and without oversampling operation is developed. This DLCT is based on the well-known CM-CC-CM decomposition, that is, implemented by two discrete chirp multiplications (CMs) and one discrete chirp convolution (CC). This decomposition doesn’t use any scaling operation which will change the sampling period or cause the interpolation error. Compared with previous works, DLCT calculated by direct summation and DLCT based on center discrete dilated Hermite functions (CDDHFs), the proposed method implemented by FFTs has much lower computational complexity. The relation between the proposed DLCT and the continuous LCT is also derived to approximate the samples of the continuous LCT. Simulation results show that the proposed method somewhat outperforms the CDDHFs-based method in the approximation accuracy. Besides, the proposed method has approximate additivity property with error as small as the CDDHFs-based method. Most importantly, the proposed method has perfect reversibility, which doesn’t hold in many existing DLCTs. With this property, it is unnecessary to develop the inverse DLCT additionally because it can be replaced by the forward DLCT.

{keywords}

ABCD transform, affine Fourier transform, fractional Fourier transform, linear canonical transform, quadratic-phase integrals.

## I Introduction

The linear canonical transform (LCT), first introduced in [1, 2], is a parameterized general linear integral transform with three degrees of freedom. The LCT unifies a variety of transforms from the well-known Fourier transform (FT), fractional Fourier transform (FRFT) and Fresnel transform (also known as chirp convolution (CC)) to simple operations such as scaling and chirp multiplication (CM) [3, 4, 5]. The LCT is an important tool in optics because the paraxial light propagation through a first-order optical system can be modeled by the LCT [6, 7, 3]. Besides, as a generalization of the transforms mentioned above, the LCT could be more useful and attractive in many signal processing applications including filter design, radar system analysis, signal synthesis, time-frequency analysis, phase reconstruction, pattern recognition, graded index media analysis, encryption and modulation [8, 9, 10, 11, 12, 13, 14]. In some papers, the LCT is also called Collins formula [1], affine Fourier transform [15], almost-Fresnel transformations [16], generalized Fresnel transforms[17], ABCD transforms [18] or quadratic-phase integrals [3].

In this paper, the definition of the LCT with four parameters [3, 19, 14] is used:

 X(u)≜OM\textmdLCT{x(t)} =⎧⎪ ⎪⎨⎪ ⎪⎩√1jb∞∫−∞ej2π(d2bu2−1but+a2bt2)x(t)dt, b≠0√d ejπcdu2x(du), b=0 (1)

where denotes the LCT operator and the parameter matrix is defined as

Although there are four parameters , the degree of freedom is three according to the constraint in (2). The LCT reduces to the FT when , and becomes the FRFT when . (Note that there is a constant phase difference between the LCT and FT/FRFT. One can refer to [5] for detailed description of the relations between the LCT and its special cases.) The LCT has two important and useful properties: reversibility and additivity [5, 19]. The reversibility property allows one to realize the inverse LCT (ILCT) with parameter matrix by the forward LCT with parameter matrix :

 OM\textmdILCTΔ=[OM\textmdLCT]−1=OM−1\textmdLCT. (3)

The additivity property of the LCT is given by

 OM1\textmdLCTOM2\textmdLCT=OM1×M2\textmdLCT. (4)

It implies that the cascade of several LCTs with parameter matrices can be replaced by only one LCT using parameter matrix . In other words, the LCT can be decomposed into a cascade of multiple LCTs. The reversibility property is a special case of the additivity property when . Besides, it is worth noting that the LCT is not commutative in most cases, i.e. , because .

Consider a set of sampled data with sampling period , i.e. . As , the sampled output of the LCT is given by

 X[k]=√1jb∑nej2π(d2bk2Δ2u−1bknΔuΔt+a2bn2Δ2t)x[n]Δt. (5)

This direct summation method is inefficient because of its high computational complexity. In [20], the authors proposed some conditions for the sampling periods such that (5) becomes an unitary discrete transform, however, no low-complexity implementation method is proposed. Taking the benefit of the additivity property mentioned in (4), the complexity can be reduced by decomposing the LCT into a sequence of simpler operations, such as scaling, chirp multiplication (CM), chirp convolution (CC) and FRFT [21, 22, 23, 24, 19]. In[21], two kinds of decompositions are proposed, decomposing the LCT into CMs and CCs. But there is no further discussion on the digital implementation based on the decompositions. Most of the other kinds of decompositions are introduced in [24]. The first method in [24] decomposes the LCT into scaling, CMs and Fourier transforms. The second method is based on the well-known Iwasawa decomposition [25, 26], i.e. one FRFT, one scaling and one CM. It follows that the digital computation of the LCT can be realized by discrete CM, discrete Fourier transform (DFT) and/or discrete FRFT (DFRFT). It has been shown in [10, 24] that the Wigner distribution function (WDF) of the LCT output is a linearly affine distorted version of the input. The bandwidth-time (BT) product of the output would be larger than the input. To avoid aliasing effect, oversampling is utilized to increase the number of samples (i.e. reduce the sampling period). The two methods in [24] require smaller oversampling rate than the direct computation in (5), and thus can further lower the complexity. In [19], a DLCT with no oversampling involved is proposed, which is also based on the Iwasawa decomposition. Since the Hermite functions are the eigenfunctions of the FRFT [5], the DFRFT can be implemented by an orthonormal basis of discrete Hermite functions. Since Iwasawa decomposition involves a scaling operation, the scalable discrete Hermite functions, called center discrete dilated Hermite functions (CDDHFs) [27], are used to implement the DLCT. The oversampling operation is not used in this DLCT. So the length of output signal will remain the same as the input. Besides, simulation results in [19] show that this DLCT has smaller error in the reversibility and additivity properties in comparison with Koç’s methods [24]. However, the main disadvantage is the high computational complexity on the generation of CDDHFs.

Actually, Koç’s work [24] is mainly focused on determining the sampling periods and number of samples so that the continuous LCT can be recovered from its samples. In this paper, the focus is on the development of a discrete LCT (DLCT) which is irrelevant to the sampling periods and doesn’t involve oversampling. Considering the FT, i.e. a special case of the LCT, it is expected that the DLCT can reduce to the DFT in this special case. Sampling periods and aliasing effect are not the concerns of the DFT, and thus also not of the DLCT. Since scaling operation will change the sampling period [24] or cause the interpolation error, the CM-CC-CM decomposition [21, 3] that involves no scaling is adopted in this paper. The CM-CC-CM decomposition is also used in the implementations of many other transforms such as FRFT [28, 3], complex LCT [29, 30] and gyrator transform [31]. The proposed DLCT based on CM-CC-CM decomposition consists of three discrete CMs, one DFT and one IDFT. In comparison with the DLCT based on CDDHFs [19], which is also irrelevant to sampling periods, the proposed DLCT has much less computational complexity. Besides, simulation results show that in some cases, the proposed DLCT can approximate the samples of the continuous LCT with higher accuracy and yield smaller error in the additivity property. Most importantly, the proposed DLCT has “perfect” reversibility and thus is more useful in some applications such as encryption/decryption.

This paper is organized as follows. Section II gives a review of the DLCT proposed by Pei and Lai [19] for comparison with the newly proposed DLCT. In Section III, we develop a DLCT based on the CM-CC-CM decomposition. Some special cases of the proposed DLCT are discussed. We also investigate the relation between the proposed DLCT and the continuous LCT in Section IV. In Section V, the proposed DLCT will be compared with the previous work [19]. Finally, conclusions are made in Section VI.

## Ii Review of DLCT Based on Center Discrete Dilated Hermite Functions [19]

The goal of this paper is on the development of a discrete LCT (DLCT) which is irrelevant to the sampling periods and without oversampling procedure. We first give a brief review of the DLCT in [19], which will be used for comparison with the newly proposed DLCT. In [19], the DLCT is based on the Iwasawa decomposition [25, 26]:

 M=[abcd]=[10ξ1][σ00σ−1][cosαsinα−sinαcosα]. (6)

The first matrix corresponds to chirp multiplication (CM) with chirp rate . The second matrix corresponds to scaling operation with scaling parameter . The last one corresponds to FRFT with fractional angle that satisfies and .

It has been known that Hermite functions (HFs) are the eigenfunctions of the FRFT [32, 33]. From (6), if the input is first expanded by the HFs, the LCT output can be expressed by the scaled HFs multiplied by a chirp function. In order to develop a unitary DLCT, orthonormal discrete HFs (DHFs) and orthonormal scalable DHFs are necessary. Therefore, Pei and Lai utilize the center discrete dilated Hermite functions (CDDHFs) [27], which are orthonormal and can approximate the samples of the scaled HFs. Let denote the CDDHF of order with scaling parameter . Then, is the undilated DHF of order . The DFRFT and discrete scaling are designed satisfying

 O(cosα,sinα;−sinα,cosα)\textmdDLCT {Φ\textmdCDDHFp;1[n]}=e−j(p+12)α ⋅Φ\textmdCDDHFp;1[n], (7) O(σ,0;\textmd0,σ−1)\textmdDLCT {Φ\textmdCDDHFp;1[n]}=Φ\textmdCDDHFp;σ[n], (8)

respectively. Denote as an orthonormal matrix with th column being . From (7) and (8), the matrix-computational expressions of the discrete CM, discrete scaling and DFRFT in (6) are , and , respectively, where and are diagonal matrices consisting of the eigenvalues in (7) and chirp samples, respectively:

 [Vα]p,p=e−j(p+12)α,[Cξ]k,k=ejπNξk2. (9)

Therefore, the CDDHFs-based DLCT is given by

 X=OM\textmdDLCT{x} =Cξ(EσET1)(E1VαET1)x =CξEσVαET1x, (10)

where and are vectors consisting of and the DLCT output , respectively.

This DLCT has three main disadvantages. First, simulation results show that this DLCT can approximate the continuous LCT well; however, the sampling periods are restricted to that is used when generating the CDDHFs. Second, it doesn’t satisfy the reversibility property perfectly, even though the error is smaller than other previous works. The last one is the relatively high computational complexity, even higher than the direct computation in (5). Besides, the generation of CDDHFs (i.e. ) is based on eigendecomposition, which is time-consuming, and precomputing is impractical because different is used for different (different parameter matrix ).

## Iii DLCT Based on CM-CC-CM Decomposition

In this paper, the CM-CC-CM decomposition [21, 3] is utilized. The digital computation of the LCT based on this decomposition has been mentioned by Koç et al. in one brief paragraph of [24]. However, their work concentrates on the oversampling operation in each step so that the number of samples is sufficient for recovering the corresponding continuous signal. The development of DLCT, which is irrelevant to the sampling periods and doesn’t involve oversampling, for discrete data is the main object of our interest. In short, Koç’s work is like digitally computing the FT with careful attention to the sampling periods, while our work is much like the development of DFT.

Although a variety of decompositions have been presented in [24], we adopt the CM-CC-CM decomposition because it doesn’t have scaling operation. In [24], scaling merely changes the sampling period and reinterprets the same samples with the scaled period; that is, the output of is .

### Iii-a Formulation of the DLCT For B≠0

We directly develop the DLCT in the discrete domain. Firstly, consider the DLCT based on direct summation:

 X (11) \textmdi.e.  X[k] =√1jBNN/2−1∑n=−N/2ej2πN(D2Bk2−1Bkn+A2Bn2)x[n], (12)

where and are vectors consisting of and the DLCT output , respectively; is the DLCT kernel of size ; and is the parameter matrix defined as . The range of in the summation of (12), i.e. , is replaced by if is odd. Note that would be different from the parameters used in continuous LCT. Later we will show that depends on and the sampling period . It is obvious that (12) reduces to the DFT multiplied by a constant phase when . If we want the inverse DLCT (IDLCT) given by

 x =OM−1\textmdDLCT{X}=K†x, (13) \textmdi.e.  x[n] =√1−jBNN/2−1∑k=−N/2ej2πN(−A2Bn2+1Bkn−D2Bk2)X[k], (14)

it is required that and

Unfortunately, in other cases, this definition usually violates the reversibility property, i.e. . Besides, the direct summation is quite inefficient because it involves complex multiplications.

Accordingly, we introduce the CM-CC-CM decomposition:

 M=[ABCD]=[10D−1B1][1B01][10A−1B1]. (16)

These three matrices from left to right represent chirp multiplication (CM) with chirp rate , chirp convolution (CC) with parameter , and again CM with chirp rate , respectively. The CC can be further decomposed into IDFT-CM-DFT, i.e.

 [1B01]=[0−110][10−B1][01−10]. (17)

Then, we can develop a DLCT completely composed of DFT, IDFT and discrete CMs. The three steps of computation of the proposed DLCT are shown below:

 \textmd(CM) x1[n]=ejπNA−1Bn2x[n], (18) \textmd(CC) X1[k]=1NN/2−1∑m=−N/2e−jπNBm2 ⋅⎛⎝N/2−1∑n=−N/2x1[n]e−j2πNmn⎞⎠ej2πNmk, (19) \textmd(CM) X[k]=ejπND−1Bk2X1[k]. (20)

Denote and as the DFT and IDFT matrices, and as a diagonal matrix where the -th element is . Then, the matrix form of the DLCT is given by

 X=OM\textmdDLCT{x}=CD−1BF†C−BFCA−1Bx. (21)

Since the DFT/IDFT can be implemented by fast Fourier transform (FFT), the computational complexity is apparently lower than that of the direct summation method (12). Replacing in (21) by leads to the IDLCT based on CM-CC-CM decomposition:

 x=OM−1\textmdDLCT{X}=C−A−1BF†CBFC−D−1BX. (22)

It is obvious that the proposed DLCT satisfies the reversibility property because

 x =C−A−1BF†CBFC−D−1B(CD−1BF†C−BFCA−1Bx) =OM−1\textmdDLCT{OM\textmdDLCT{x}}. (23)

### Iii-B Formulation of the DLCT For B=0

The definition in (21) is invalid when . Fortunately, if , one has and because . In this case, the following two kinds of decompositions are considered:

 M=[A0CD] =[01−10][−C−DA0] =[01−10][101D1][1−D01][10C+1D1], (24)

and

 M=[A0CD] =[0A−DC][0−110] =[10C−1A1][1A01][10−1A1][0−110]. (25)

The matrix forms of the DLCTs based on (III-B) and (III-B) are given by

 X =O(A,0;C,D)\textmdDLCT{x}=√−j FC1DF†CDFCC+1Dx, (26) X =O(A,0;C,D)\textmdDLCT{x}=√j CC−1AF†C−AFC−1AF†x, (27)

respectively. Note that (or ) is the phase difference between LCT and FT (or inverse FT).

The reversibility property for is given by

 x=O(D,0;−C,A)\textmdDLCT{O(A,0;C,D)\textmdDLCT{x}}. (28)

If both the forward and inverse transforms use the same DLCT form in (26) or (27), the reversibility property doesn’t holds due to no intermediate cancellation in between:

 x ≠√−j FC1AF†CAFC−C+1A {√−j FC1DF†CDFCC+1Dx} (29) ≠√j C−C−1DF†C−DFC−1DF† {√j CC−1AF†C−AFC−1AF†x}. (30)

However, if we let the forward and inverse transforms use different DLCT forms, reversibility can be achieved because and in between:

 x =√j C−C−1DF†C−DFC−1DF† {√−j FC1DF†CDFCC+1Dx} (31) =√−j FC1AF†CAFC−C+1A {√j CC−1AF†C−AFC−1AF†x}. (32)

Accordingly, we make the assumption that (26) is used when and (27) is used when :

 O(A,0;C,D)\textmdDLCT{x} =⎧⎪⎨⎪⎩√−j FC1DF†CDFCC+1Dx,\textmdfor |A|>|D|√j CC−1AF†C−AFC−1AF†x,\textmdfor |A|<|D|. (33)

(Alternatively, one can choose using another assumption that (26) is used when and (27) is used when .) For example, consider and then . The DLCT with uses the form for in (III-B) because . The IDLCT, i.e. DLCT with , uses the form for in (III-B) because . The reversibility holds obviously because

 x =O(0.5,0;−1,2)\textmdDLCT{O(2,0;1,0.5)\textmdDLCT{x}} =√j C−4F†C−0.5FC−2F†(√−j FC2F†C0.5FC4x). (34)

### Iii-C Special Cases of the DLCT

In this subsection, we discuss some special cases of the proposed DLCT, including DFRFT, discrete Fresnel transform and discrete scaling operation.

#### Discrete fractional Fourier transform

In [5], it has been shown that the FRFT of fractional angle is the special case of the LCT of with some phase difference:

 Oα\textmdFRFT{x(t)}=ejα2O(cosα,sinα;−sinα,cosα)\textmdLCT{x(t)}. (35)

Discard the case that , i.e. , because it is just an identity operation. Then, the DFRFT can be implemented by the proposed DLCT defined in (21) with the same phase difference . That is,

 Oα\textmdDFRFT{x} Δ=ejα2O(cosα,sinα;−sinα,cosα)\textmdDLCT{x} =ejα2Ccosα−1sinαF†C−sinαFCcosα−1sinαx =ejα2C−tanα2F†C−sinαFC−tanα2x. (36)

In [28], the authors proposed two digital computation methods for the FRFT. The second method is similar to the direct summation method in (12) with and . The first method presents similar notion as that in (III-C1), i.e. decomposing the FRFT into CM-CC-CM. However, the difference between them is the implementation of the CC. In [28], the CC is digitally implemented by a discrete linear convolution and is recommended using FFTs. So the inputs of the convolution have to be padded with zeros before passing through FFTs. The CC in (III-C1) is digitally implemented by two FFTs and one discrete CM, like a discrete circular convolution without zero-padding.

#### Discrete Fresnel transform

The 1-D Fresnel transform with wavelength and propagation distance , denoted by , is also a special case of the LCT where [5], i.e. a CC,

 Oλ,z\textmdFresnel{x(t)}=ejπzλO(1,λz;0,1)\textmdLCT{x(t)}. (37)

Accordingly, we can design the discrete Fresnel transform as follows using the proposed DLCT defined in (21):

 Oλ,z\textmdDFresnel{x} Δ=ejπzλO(1,λz;0,1)\textmdDLCT{x} =ejπzλF†C−λzFx. (38)

#### Discrete scaling

When , the LCT reduces to the scaling operation with scaling parameter , denoted by :

 Oσ\textmdScal{x(t)}=x(t/σ)=O(σ,0;σ−1,0)\textmdLCT{x(t)}. (39)

Recall the proposed DLCT for shown in (III-B). The discrete scaling operation is given by

 Oσ\textmdDScal{x}=⎧⎨⎩√−jFCσF†C1σFCσx,|σ|>1√jC−1σF†C−σFC−1σF†x,|σ|<1. (40)

As mentioned in (III-B)-(32), the above two kinds of decompositions are used for perfect reversibility.

#### Other discrete operations

The FT, inverse FT and CM are also the special cases of the LCT with parameter matrix being , and , respectively. The discrete versions of these operations are simply the DFT (), IDFT () and discrete CM (), respectively, without the need of the proposed DLCT.

## Iv Relation Between Proposed DLCT and Continuous LCT

In this section, we discuss the connections between the proposed DLCT and continuous LCT, including derivation of the DLCT from the continuous LCT, relation between of the DLCT and of the continuous LCT, and oversampling for the DLCT to approximate the samples of the continuous LCT. Some related works regarding sampling and oversampling of the LCT include [34, 23, 35, 24, 36]. In this paper, we discuss sampling and oversampling from the point of view of CM-CC-CM decomposition.

### Iv-a Derivation of Proposed DLCT From Continuous LCT

Since the proposed DLCT is based on the CM-CC-CM decomposition, consider the expression of the continuous LCT based on the same decomposition:

 \textmd(CM) x1(t)=ejπa−1bt2x(t), (41) \textmd(CC) X1(u)=∞∫−∞e−jπbf2⎛⎜⎝∞∫−∞x1(t)e−j2πftdt⎞⎟⎠ej2πfudf, (42) \textmd(CM) X(u)=ejπd−1bu2X1(u). (43)

Let and denote the FTs of and , respectively:

 ˆx1(f)=F{x1(t)}\textmdandˆX1(f)=F{X1(u)}, (44)

where denotes the FT operation. From (44) and (42), it follows that

 ˆX1(f)=e−jπbf2ˆx1(f). (45)

Sample and with the same sampling period . Replacing and in (44) by their discrete samples and leads to

 ˆx1,1Δ(f)Δ=∞∑l=−∞ˆx1(f−lΔ)=Δ∞∑n=−∞x1(nΔ)e−j2πfnΔ, (46)
 ˆX1,1Δ(f)Δ=∞∑l=−∞ˆX1(f−lΔ)=Δ∞∑k=−∞X1(kΔ)e−j2πfkΔ. (47)

To further simplify the above computations, we sample and by sampling period where is some positive integer. Then, with , (46) becomes

 ˆx1,1Δ(mNΔ) =Δ∞∑n=−∞x1(nΔ)e−j2πNmn =ΔN/2−1∑n=−N/2x1,NΔ(nΔ)e−j2πNmn, (48)

where is the periodic summation of with period :

 x1,NΔ(t)Δ=∞∑l=−∞x1(t−lNΔ), (49)

Similarly, sampling (47) with leads to

 ˆX1,1Δ(mNΔ) =ΔN/2−1∑k=−N/2X1,NΔ(kΔ)e−j2πNmk, (50)

where is the periodic summation of :

 X1,NΔ(u)Δ=∞∑l=−∞X1(u−lNΔ). (51)

The inverse transform of (50) is given by

 X1,NΔ(kΔ)=1NΔN/2−1∑m=−N/2ˆX1,1Δ(mNΔ)ej2πNmk. (52)

Finally, a DLCT can be developed by the following five steps (the combination of the 2nd to 4th steps corresponds to the CC procedure):

 x1,NΔ(nΔ) =ejπNa−1bNΔ2n2x(nΔ), (53) ˆx1,1Δ(mNΔ) =ΔN/2−1∑n=−N/2x1,NΔ(nΔ)e−j2πNmn, (54) ˆX1,1Δ(mNΔ) =e−jπNb1NΔ2m2ˆx1,1Δ(mNΔ), (55)
 X1,NΔ(kΔ) =1NΔN/2−1∑m=−N/2ˆX1,1Δ(mNΔ)ej2πNmk, (56) X[k] =ejπNd−1bNΔ2k2X1,NΔ(kΔ), (57)

where . The assumptions in (53), (55) and (57) may be inconsistent with the relations in (41), (45) and (43). If so, the output of the DLCT will be different from the sampled output of the continuous LCT. In next subsection, we will discuss the conditions such that the above DLCT can produce an accurate approximation to .

Compared with (18)-(20), the above DLCT is equivalent to the proposed DLCT if , and . It implies that the relation between the parameter matrix of the proposed DLCT and the parameter matrix of the continuous LCT is given by

 [ABCD]=[abNΔ2cNΔ2d], (58)

where when .

### Iv-B Oversampling for Approximating the Samples of Continuous LCT

In (III-A), (III-B) and (32), it has been proved that the proposed DLCT is always reversible. However, the output of the proposed DLCT may be different from the sampled output of the continuous LCT because of some aliasing and overlapping problems. Recall the DLCT in (53)-(57). If we want , (53), (55) and (57) have to be consistent with (41), (45) and (43), respectively. That is, for , we need

 x1,NΔ(nΔ) =∞∑l=−∞x1(nΔ−lNΔ)=x1(nΔ), (59) ˆx1,1Δ(mNΔ) =∞∑