1

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

arXiv:1709.06222v1 [cs.IT] 19 Sep 2017

Soo-Chang Pei, Life Fellow, IEEE, and Shih-Gu Huang

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-CCCM 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.
Index Terms—ABCD transform, affine Fourier transform,
fractional Fourier transform, linear canonical transform,
quadratic-phase integrals.

I. I NTRODUCTION
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]–[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 [3], [6], [7]. 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]–[14]. In some
Copyright (c) 2015 IEEE. Personal use of this material is permitted.
However, permission to use this material for any other purposes must be
obtained from the IEEE by sending a request to pubs-permissions@ieee.org.
This work was supported by the Ministry of Science and Technology,
Taiwan, under Contracts MOST 104-2221-E-002-096-MY3 and MOST 1042221-E-002-006.
S. C. Pei is with the Department of Electrical Engineering & Graduate
Institute of Communication Engineering, National Taiwan University, Taipei
10617, Taiwan (e-mail: peisc@ntu.edu.tw).
S.-G. Huang is with the Graduate Institute of Communication Engineering, National Taiwan University, Taipei 10617, Taiwan (e-mail:
d98942023@ntu.edu.tw).

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], [14], [19] is used:
M
X(u) , OLCT
{x(t)}
q
∞
R
d 2
1
a 2
1

ej2π( 2b u − b ut+ 2b t ) x(t) dt, b 6= 0
jb
= √ −∞
2

d ejπ cd u x(du),
b=0

(1)

M
where OLCT
denotes the LCT operator and the 2×2 parameter
matrix M is defined as


a b
M = (a, b; c, d) =
and ad − bc = 1.
(2)
c d

Although there are four parameters (a, b; c, d), the degree of
freedom is three according to the constraint in (2). The LCT
reduces to the FT when M = (0, 1; −1, 0), and becomes the
FRFT when M = (cos α, sin α; − sin α, cos α). (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 M by the forward
LCT with parameter matrix M−1 :
∆  M −1
M−1
M
= OLCT
.
(3)
= OLCT
OILCT
The additivity property of the LCT is given by
M1
M1 ×M2
M2
OLCT
.
= OLCT
OLCT

(4)

It implies that the cascade of several LCTs with parameter
matrices M1 , M2 , . . . , Mk can be replaced by only one LCT
using parameter matrix M = M1 × M2 × · · · × Mk . 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 M1 = M−1
2 . Besides, it is worth
noting that the LCT is not commutative in most cases, i.e.
M1 M2
M2 M1
OLCT
OLCT 6= OLCT
OLCT , because M1 M2 6= M2 M1 .
Consider a set of sampled data x[n] with sampling period
∆t , i.e. x[n] = x(n∆t ). As b 6= 0, the sampled output of the
LCT X[k] = X(k∆u ) is given by
r
a
d 2 2
1 X j2π( 2b
k ∆u − 1b kn∆u ∆t + 2b
n2 ∆2t )
x[n]∆t .
e
X[k] =
jb n
(5)

2

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 [19], [21]–[24]. 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 wellknown 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 Koç’s
methods [24]. However, the main disadvantage is the high
computational complexity on the generation of CDDHFs.
Actually, Koç’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 [3], [21] 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 [3], [28], 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. R EVIEW OF DLCT BASED ON C ENTER D ISCRETE
D ILATED H ERMITE F UNCTIONS [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]:

 



a b
1 0 σ 0
cos α sin α
M=
=
.
(6)
c d
ξ 1 0 σ −1 − sin α cos α
The first matrix corresponds to chirp multiplication (CM)
with chirp rate ξ = (ac + bd)/(a2 + b2 ). The second matrix
corresponds
to scaling operation with scaling parameter σ =
√
a2 + b2 . The last one corresponds to FRFT with fractional
angle α that satisfies cos α = a/σ and sin α = b/σ.
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 ΦCDDHF
[n] denote the
p;σ
CDDHF of order p with scaling parameter σ. Then, ΦCDDHF
[n]
p;1
is the undilated DHF of order p. The DFRFT and discrete
scaling are designed satisfying
1
(cos α,sin α;− sin α,cos α)  CDDHF
ODLCT
Φp;1 [n] = e−j (p+ 2 )α
(σ,0; 0,σ

ODLCT

−1

)

· ΦCDDHF
[n], (7)
p;1

 CDDHF
Φp;1 [n] = ΦCDDHF
[n],
p;σ

(8)

respectively. Denote Eσ as an N ×N orthonormal matrix with
pth column being ΦCDDHF
[n]. From (7) and (8), the matrixp;σ
computational expressions of the discrete CM, discrete scaling
and DFRFT in (6) are Cξ , Eσ ET1 and E1 Vα ET1 , respectively,
where Vα and Cξ are N × N diagonal matrices consisting of

3

the eigenvalues in (7) and chirp samples, respectively:
[Vα ]p,p = e−j (

p+ 12

)α ,

[Cξ ]k,k = ej

π
N

ξk

2

.

Therefore, the CDDHFs-based DLCT is given by


M
X = ODLCT
{x} = Cξ Eσ ET1 E1 Vα ET1 x
= Cξ Eσ Vα ET1 x,

(9)

(10)

where x and X are N × 1 vectors consisting of x[n] and the
DLCT output X[k], respectively.
This DLCT has three main disadvantages. First, simulation
results show that this DLCT can approximate the continuous
LCT well; however,
p the sampling periods are restricted to
1/N that is used when generating the
∆t = ∆u =
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. Eσ ) is based on
eigendecomposition, which is time-consuming, and precomputing Eσ is impractical because different Eσ is used for
different σ (different parameter matrix M).
III. DLCT BASED

ON

CM-CC-CM D ECOMPOSITION

In this paper, the CM-CC-CM decomposition [3], [21] is
utilized. The digital computation of the LCT based on this
decomposition has been mentioned by Koç 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, Koç’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 x[n] = x(n∆) is X[k] =
x[k] = x(kd∆).

if n is odd. Note that (A, B; C, D) would be different from
the parameters (a, b; c, d) used in continuous LCT. Later we
will show that (A, B; C, D) depends on (a, b; c, d) and the
sampling period ∆. It is obvious
√that (12) reduces to the DFT
multiplied by a constant phase −j when M = (0, 1; −1, 0).
If we want the inverse DLCT (IDLCT) given by
M
x = ODLCT
{X} = K† x,
(13)
r
N/2−1
X
2
A
1
D 2
2π
1
ej N (− 2B n + B kn− 2B k ) X[k],
i.e. x[n] =
−jBN
−1

k=−N/2

(14)

it is required that M−1 = (D, −B; −C, A) and
AD − BC = 1.

Unfortunately, in other cases, this definitionusually violates
M−1
M
the reversibility property, i.e. x[n] 6= ODLCT
ODLCT
{x[n]} .
Besides, the direct summation is quite inefficient because it
involves N 2 complex multiplications.
Accordingly, we introduce the CM-CC-CM decomposition:

 



1 0 1 B
1 0
A B
M=
= D−1
.
(16)
C D
1 0 1 A−1
1
B
B
These three matrices from left to right represent chirp multiplication (CM) with chirp rate (D−1)/B, chirp convolution (CC)
with parameter B, and again CM with chirp rate (A − 1)/B,
respectively. The CC can be further decomposed into IDFTCM-DFT, i.e.

 



1 B
0 −1
1 0
0 1
=
.
(17)
0 1
1 0
−B 1 −1 0
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:
π A−1 2
B n

(CM) x1 [n] = ej N
(CC) X1 [k] =

1
N

We directly develop the DLCT in the discrete domain.
Firstly, consider the DLCT based on direct summation:
M
X = ODLCT
{x} = Kx,
(11)
r
N/2−1
X
2
2π
D 2
1
A
1
ej N ( 2B k − B kn+ 2B n ) x[n],
i.e. X[k] =
jBN
n=−N/2

(12)

where x and X are N × 1 vectors consisting of x[n] and the
DLCT output X[k], respectively; K is the DLCT kernel of
size N × N ; and M is the 2 × 2 parameter matrix defined as
M = (A, B; C, D). The range of n in the summation of (12),
i.e. [−N/2, N/2 − 1], is replaced by [−(N − 1)/2, (N − 1)/2]

x[n],

(18)

N/2−1

X

π

e−j N Bm

2

m=−N/2



(CM) X[k] = e

A. Formulation of the DLCT For B 6= 0

(15)

·

π
jN

N/2−1

X

x1 [n]e−j

n=−N/2

D−1 2
B k

X1 [k].

2π
N mn



 ej 2π
N mk , (19)

(20)

Denote F and F† as the DFT and IDFT matrices, and Cξ as
a diagonal matrix where the (n, n)-th element is [Cξ ]n,n =
2
π
ej N ξn . Then, the matrix form of the DLCT is given by
M
X = ODLCT
{x} = C D−1 F† C−B FC A−1 x.
B

B

(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 M
in (21) by M−1 leads to the IDLCT based on CM-CC-CM
decomposition:
M
x = ODLCT
{X} = C− A−1 F† CB FC− D−1 X.
−1

B

B

(22)

It is obvious that the proposed DLCT satisfies the reversibility

4

property because
x = C− A−1 F† CB FC− D−1
B
B
 M
M−1
= ODLCT ODLCT {x} .



C D−1 F† C−B FC A−1 x
B

B

(23)

The definition in (21) is invalid when B = 0. Fortunately, if
B = 0, one has A 6= 0 and D 6= 0 because AD−BC = AD =
1. In this case, the following two kinds of decompositions are
considered:

 


A 0
0 1 −C −D
M=
=
C D
−1 0
A
0





1 0
0 1 1 0 1 −D
, (24)
=
1
C+1
1 0 1
−1 0 D
1
D
M=



 

A 0
0 A 0
=
C D
−D C 1


1 0 1
= C−1
1 0
A

(A,0;C,D)

{x}
ODLCT
(√
−j FC D1 F† CD FC C+1 x, for |A| > |D|
D
= √
. (33)
j C C−1 F† C−A FC− A1 F† x, for |A| < |D|
A

B. Formulation of the DLCT For B = 0

and

|A| > |D| and (27) is used when |A| < |D|:


−1
0



A
1 0 0 −1
. (25)
1
1 −A
1 1 0

The matrix forms of the DLCTs based on (24) and (25) are
given by
p
(A,0;C,D)
{x} = −j FC D1 F† CD FC C+1 x, (26)
X = ODLCT
D
p
(A,0;C,D)
{x} = j C C−1 F† C−A FC− A1 F† x, (27)
X = ODLCT
A
√
√
respectively. Note that −j (or j) is the phase difference
between LCT and FT (or inverse FT).
The reversibility property for B = 0 is given by
n
o
(D,0;−C,A)
(A,0;C,D)
x = ODLCT
ODLCT
{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:
p
x 6= −j FC A1 F† CA FC −C+1
np A
o
−j FC D1 F† CD FC C+1 x
(29)
D
p
6= j C −C−1 F† C−D FC− D1 F†
D
np
o
j C C−1 F† C−A FC− A1 F† x . (30)
A

However, if we let the forward and inverse transforms use
different DLCT forms, reversibility can be achieved because
F† F = FF† = I and Cξ C−ξ = I in between:
p
x = j C −C−1 F† C−D FC− D1 F†
D
np
o
−j FC D1 F† CD FC C+1 x
(31)
D
p
= −j FC A1 F† CA FC −C+1
A
o
np
j C C−1 F† C−A FC− A1 F† x . (32)
A

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

(Alternatively, one can choose using another assumption that
(26) is used when |A| < |D| and (27) is used when |A| >
|D|.) For example, consider (A, B; C, D) = (2, 0; 1, 0.5)
and then (A, B; C, D)−1 = (0.5, 0; −1, 2). The DLCT with
(2, 0; 1, 0.5) uses the form for |A| > |D| in (33) because
2 > 0.5. The IDLCT, i.e. DLCT with (0.5, 0; −1, 2), uses the
form for |A| < |D| in (33) because 0.5 < 2. The reversibility
holds obviously because
n
o
(0.5,0;−1,2)
(2,0;1,0.5)
x = ODLCT
ODLCT
{x}
p

p
−j FC2 F† C0.5 FC4 x .
= j C−4 F† C−0.5 FC−2 F†
(34)

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.
1) Discrete fractional Fourier transform: In [5], it has been
shown that the FRFT of fractional angle α is the special case
of the LCT of M = (cos α, sin α; − sin α, cos α) with some
phase difference:
α

(cos α,sin α;− sin α,cos α)

α
{x(t)} = ej 2 OLCT
OFRFT

{x(t)} .

(35)

Discard the case that α = 0, i.e. M = (1, 0; 0, 1), 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 ej 2 . That is,
∆

α

(cos α,sin α;− sin α,cos α)

α
ODFRFT
{x} = ej 2 ODLCT
α
2

=e

j

=e

jα
2

†

{x}

C cos α−1 F C− sin α FC cos α−1 x
sin α

sin α

C− tan α2 F† C− sin α FC− tan α2 x.

(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 (A, B; C, D) =
(cos α, sin α; − sin α, cos α) and N = (2∆x )2 . The first
method presents similar notion as that in (36), 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 (36) is digitally implemented by two FFTs and
one discrete CM, like a discrete circular convolution without
zero-padding.
2) Discrete Fresnel transform: The 1-D Fresnel transform
with wavelength λ and propagation distance z, denoted by
λ,z
OFresnel
, is also a special case of the LCT where M =

5

(1, λz; 0, 1) [5], i.e. a CC,
λ,z
OFresnel

{x(t)} = ej

πz
λ

(1,λz;0,1)
OLCT

{x(t)} .

(37)

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

λ,z
ODFresnel
{x} = ej

=e

πz
λ

j πz
λ

(1,λz;0,1)

ODLCT
†

{x}

F C−λz Fx.

(38)

3) Discrete scaling: When M = (σ, 0; σ −1 , 0), the LCT
reduces to the scaling operation with scaling parameter σ,
σ
denoted by OScal
:
σ
OScal

{x(t)} = x(t/σ) =

(σ,0;σ−1 ,0)
OLCT
{x(t)} .

(39)

Recall the proposed DLCT for B = 0 shown in (33). The
σ
discrete scaling operation ODScal
is given by
( √
−j FCσ F† C σ1 FCσ x, |σ| > 1
σ
ODScal
{x} = √
. (40)
j C− σ1 F† C−σ FC− σ1 F† x, |σ| < 1
As mentioned in (29)-(32), the above two kinds of decompositions are used for perfect reversibility.
4) Other discrete operations: The FT, inverse FT and CM
are also the special cases of the LCT with parameter matrix
being (0, 1; −1, 0), (0, −1; 1, 0) and (1, 0; ξ, 1), respectively.
The discrete versions of these operations are simply the DFT
(F), IDFT (F† ) and discrete CM (Cξ ), respectively, without
the need of the proposed DLCT.

where F denotes the FT operation. From (44) and (42), it
follows that
b1 (f ) = e−jπbf 2 x
X
b1 (f ).

Sample x1 (t) and X1 (u) with the same sampling period ∆t =
∆u = ∆. Replacing x1 (t) and X1 (u) in (44) by their discrete
samples x1 (n∆) and X1 (k∆) leads to


∞
∞
X
X
l
∆
1
=∆
x1 (n∆)e−j2πf n∆ ,
x
b1, ∆ (f ) =
x
b1 f −
∆
n=−∞
l=−∞

(46)

∆
b 1 (f ) =
X
1, ∆

∞
X



∞
X
b1 f − l = ∆
X
X1 (k∆)e−j2πf k∆ .
∆

l=−∞

k=−∞

(47)

To further simplify the above computations, we sample
b 1 (f ) by sampling period ∆f = 1 where
x
b1, ∆1 (f ) and X
1, ∆
N∆
N is some positive integer. Then, with f = Nm∆ , (46) becomes
x
b1, ∆1

∞
 m 
X
2π
=∆
x1 (n∆)e−j N mn
N∆
n=−∞
N/2−1

X

=∆

2π

x1,N ∆ (n∆) e−j N mn ,

where x1,N ∆ (t) is the periodic summation of x1 (t) with
period N ∆:
x1,N ∆ (t) =

AND

∞
X

l=−∞

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
(A, B; C, D) of the DLCT and (a, b; c, d) 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 [23], [24],
[34]–[36]. In this paper, we discuss sampling and oversampling from the point of view of CM-CC-CM decomposition.

(48)

n=−N/2

∆

IV. R ELATION B ETWEEN P ROPOSED DLCT
C ONTINUOUS LCT

(45)

x1 (t − lN ∆),

Similarly, sampling (47) with f =
b1, 1
X
∆

 m 
=∆
N∆

m
N∆

(49)

leads to

N/2−1

X

2π

X1,N ∆ (k∆) e−j N mk ,

(50)

k=−N/2

where X1,N ∆ (u) is the periodic summation of X1 (u):
∆

X1,N ∆ (u) =

∞
X

l=−∞

X1 (u − lN ∆).

(51)

The inverse transform of (50) is given by
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:
a−1 2

(41)
(CM) x1 (t) = ejπ b t x(t),
 ∞

∞
Z
Z
2
(CC) X1 (u) =
e−jπbf 
x1 (t)e−j2πf t dtej2πf u df,
−∞

(CM) X(u) = e

jπ

−∞

d−1 2
b u

X1 (u).

(42)

(43)

b1 (f ) denote the FTs of x1 (t) and X1 (u),
Let x
b1 (f ) and X
respectively:
x
b1 (f ) = F {x1 (t)}

b1 (f ) = F {X1 (u)} ,
and X

(44)

1
X1,N ∆ (k∆) =
N∆

N/2−1

X

m=−N/2

b 1
X
1, ∆

 m  2π
ej N mk .
N∆

(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):
2 2
π a−1
b N∆ n

x1,N ∆ (n∆) = ej N
x
b1, ∆1

b1, 1
X
∆

x(n∆),

N/2−1
 m 
X
2π
x1,N ∆ (n∆)e−j N mn ,
=∆
N∆
n=−N/2
 m 
 m 
π
b N 1∆2 m2
−j N
=e
x
b1, ∆1
,
N∆
N∆

(53)
(54)
(55)

6

1
X1,N ∆ (k∆) =
N∆
X[k] = e

π
jN

N/2−1

X

m=−N/2

d−1
2 2
b N∆ k

b1, 1
X
∆

 m  2π
ej N mk ,
N∆

X1,N ∆ (k∆) ,

Compared with (18)-(20), the above DLCT is equivalent to
A−1
1
2
the proposed DLCT if a−1
b N∆ =
B , b N ∆2 = B and
D−1
d−1
2
b N∆ =
B . It implies that the relation between the
parameter matrix (A, B; C, D) of the proposed DLCT and the
parameter matrix (a, b; c, d) of the continuous LCT is given
by

 

b
A B
a
2
N
∆
=
,
(58)
C D
cN ∆2 d
p
where (A, B; C, D) = (a, b; c, d) when ∆ = 1/N .

B. Oversampling for Approximating the Samples of Continuous LCT
In (23), (31) 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 X[k] ≈ X(k∆), (53), (55) and (57) have to be
consistent with (41), (45) and (43), respectively. That is, for
− N2 ≤ n, m, k ≤ N2 − 1, we need
∞
X

 m 
=
N∆

l=−∞
∞
X

X1,N ∆ (k∆) =

l=−∞
∞
X

x
b1, ∆1

b1, 1
X
∆

l=−∞

x1 (n∆ − lN ∆) = x1 (n∆),
x
b1



m
l
−
N∆ ∆



≈x
b1

(59)

 m 
, (60)
N∆



∞


 m 
X
m
l
b1 m , (61)
b
=
≈X
X1
−
N∆
N∆ ∆
N∆
l=−∞

f

(57)

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

x1,N ∆ (n∆) =

x1 (t)
¯ a−1 ¯
¯
¯T + F f
b

x(t)

(56)

X1 (k∆ − lN ∆) = X1 (k∆) .

(62)

The above four conditions imply that x1 (t) and X1 (t) have
to be time-limited to [− N2∆ , N2∆ ] and approximately band1
1
, 2∆
]. (A time-limited signal cannot be comlimited to [− 2∆
pletely band-limbed.) According to (41), x1 (t) has the same
time duration as x(t). And according to (45), x1 (t) and X1 (t)
have the same bandwidth. Therefore, N and ∆ need to satisfy
the following conditions such that the output of the DLCT can

F
t

t

T
T
(b)

(a)

X1 (u)
¯ a−1 ¯
¯T + F f
¯
b

X (u)
f

|c|T + |d|F

u

u

|a|T + |b|F

|a|T + |b|F

(c)

(d)

Fig. 1.
Effect of each step of the LCT with M = (a, b; c, d) in timefrequency plane: (a) the original signal x(t), (b) x1 (t) i.e. x(t) after the CM
step, (c) X1 (u) i.e. x1 (t) after the CC step and (d) X(u) i.e. X1 (u) after the
CM step. T and F are the time duration and bandwidth of x(t), respectively.

approximate the sampled output of the continuous LCT:
N∆
,
(63)
2
1
∆
x
b1 (f ) = F {x1 (t)} ≈ 0 for |f | ≥
,
(64)
2∆
N∆
X1 (u) = 0 for |u| ≥
.
(65)
2
It can be found that the bandwidths of x(t) and X(u) are not
taken into account. However, if one wants to reconstruct X(u)
from the DLCT output X[k], it is required that
x(t) = 0 for |t| ≥

1
.
(66)
2∆
Plus, if one wants to recover x(t) from x[n] after some inverse
DLCT work, sampling period ∆ should also satisfy
F {X(u)} ≈ 0 for |f | ≥

F {x(t)} ≈ 0 for |f | ≥

1
.
2∆

(67)

Consider that the input signal x(t) has time duration T and
bandwidth F , as shown in Fig. 1(a). In [10], [24], it has been
shown that the relation of the LCT to the Wigner distribution
function (WDF) is given by
WY (u, v) = Wy (du − bv, −cu + av) ,

(68)

where WY and Wy are the WDFs of Y (u) and y(t), respectively, and Y (u) is the LCT of y(t) with parameter matrix
M = (a, b; c, d). The CM with chirp rate (a − 1)/b, i.e. M =
(1, 0; (a − 1)/b, 1), will lead to shearing in frequency domain.
T + F while
Therefore, the bandwidth of x1 (t) become a−1
b
the time duration remains the same as x(t), as shown in
Fig. 1(b). In contrast, the CC, i.e. M = (1, b; 0, 1), will lead

7

x(t)
f
|c0 |T0 + |b0 c0 + 1|F0

x0 (t)

Then, according to (68), the time duration T and bandwidth
F of x(t) are given by

f

p2

T = T0 + |b0 |F0 ,
F = |c0 |T0 + |b0 c0 + 1|F0 .

F0

p1
t

p3
p4

t

(a)

M0
Since x(t) = OLCT
{x0 (t)}, the relation of X(u) and x0 (t) is
M·M0
M
X(u) = OLCT
{x(t)} = OLCT
{x0 (t)}.

T0

T0 + |b0 |F0
(b)

(73)
(74)

Assume M1 = M · M0 , i.e.
 



b0
a b ∆ a b 1
.
M1 = 1 1 =
c1 d1
c d c0 b 0 c0 + 1

(75)

(76)

Fig. 2.
Time-frequency distributions: (a) the original signal x(t); and
(b) x0 (t) i.e. the LCT output of x(t) with parameter matrix M0 =
(1, b0 ; c0 , b0 c0 + 1). The M0 is well-designed such that the time-frequency
distribution of x0 (t) is a rectangle with time duration T0 and bandwidth F0 .

In Fig. 1, replacing x(t) by x0 (t) and (a, b; c, d) by
(a1 , b1 ; c1 , d1 ), the time durations and bandwidths of x1 (t),
X1 (u) and X(u) can be expressed in terms of T0 and F0 .
Then, (69), (70) and (71) become

to shearing in time domain and change the time duration
to |a|T + |b|F , as shown in Fig. 1(c). Fig. 1(d) shows that
the CM changes the bandwidth again to |c|T + |d|F . With
these information of the time durations and bandwidths of
x(t), x1 (t), X1 (u) and X(u), one can determine the sampling
period ∆ and the number of samples N . According to (64),
∆ needs to satisfy

1
a1 − 1
T0 + F0 ,
(77)
≥
∆
b1


a1 − 1
1
T0 + F0 , |c1 |T0 + |d1 |F0 , F , (78)
≥ max
∆
b1
1
max {T, |a1 |T0 + |b1 |F0 } ,
(79)
N≥
∆
respectively, where T and F are given in (73) and (74). Note
that the constraint (78) instead of (77) is used when one wants
to reconstruct X(u) from X[k] and x(t) from x[n]. It is
possible that F is equal to or larger than the right-hand sides
of (77) and (78) when the input x(t) has larger bandwidth
than the intermediate result x1 (t) after the CM step and the
final LCT output X(u). If so, oversampling is unnecessary.
We can summarize the results in (70) and (78) simply as
follows. Let F , F1 and Fout denote the bandwidths of the input
signal x(t), intermediate result x1 (t) after the CM step and
the LCT output X(u), respectively.
• If F < F1 or F < Fout , then oversampling is required
due to LCT bandwidth expansion.
• If F ≥ F1 and F ≥ Fout , then oversampling is not
required due to LCT bandwidth compression.
If one wants to calculate the samples of the continuous LCT
by the proposed DLCT, once the sampling period ∆ and
the number of samples N are determined, the parameter
matrix of the proposed DLCT can be obtained from (58), i.e.
(A, B; C, D) = (a, b/(N ∆2 ); cN ∆2 , d).
Note that oversampling for decreasing ∆ and zero-padding
for increasing N are classified as data preprocessors, not parts
of the proposed DLCT, because they are totally unnecessary
in some applications such as data encryption and decryption.

1
a−1
≥
T + F.
∆
b

(69)

If the additional conditions (66) and (67) are also considered,
we need


a−1
1
≥ max
T + F, |c|T + |d|F, F
∆
b


a−1
T + F, |c|T + |d|F .
(70)
= max
b
Once ∆ is determined, N can be determined based on (63)
and (65):
1
max {T, |a|T + |b|F } .
(71)
∆
Except when a = 1, no mater which one of (69) and (70) is
adopted, oversampling is required if x(t) is originally sampled
based on its bandwidth, i.e. ∆ = 1/F .
In fact, oversampling would be unnecessary if more timefrequency information of x(t) is provided. Consider that the
time-frequency distribution of x(t) is within some fundamental parallelogram. For example, assume the time-frequency
distribution of x(t) is as shown in Fig. 2(a). And assume
x0 (t) is the LCT output of x(t) with parameter matrix
M0 = (1, b0 ; c0 , b0 c0 + 1). If M0 is well-designed, the timefrequency distribution of x0 (t) will be a rectangle with time
duration T0 and bandwidth F0 , as shown in Fig. 2(b). Given
the position information of the four vertices p1 , p2 , p3 , p4 of
the fundamental parallelogram in Fig. 2(a), say p1 = (t1 , f1 )
and p2 = (t2 , f2 ), we have
N≥

T 0 = t1 − t2 ,
F0 = f1 + f2 − c0 (t1 + t2 ),

f1 − f2
,
T0
t1 + t2
b0 =
.
F0

c0 =

(72)

V. C OMPARISONS B ETWEEN P ROPOSED DLCT
CDDHF S - BASED DLCT [19]

AND

In this section, the proposed DLCT will be compared with
the previous work, CDDHFs-based DLCT in [19], which is
also irrelevant to the sampling periods and without oversampling operation. In the following, the comparisons in
computational complexity, accuracy of the approximation to
sampled continuous LCT, additivity property and reversibility
property are presented.

8

TABLE I
[19] AND PROPOSED DLCT
Complexity1

Direct sum

Kx

4N 2

CDDHFs-based2

Cξ Eσ Vα ET
1x

4N 2 + 8N

Proposed (B=0)

F† C

C D−1
12N + 4N log2 N
−B FC A−1 x
B
√ B
†
−jFC 1 F CD FC C+1 x
D
D
12N + 6N log2 N
√
jC C−1 F† C−A FC− 1 F† x
A

A

1. List the number of real multiplications.
2. The complexity of matrix eigendecomposition for Eσ and E1 is not
included.

A. Computational Complexity

NMSE

10

Proposed DLCT

0

10

ξ
σ
α

−10

Matrix form

Proposed (B6=0)

CDDHFs−based DLCT

0

10

C OMPLEXITY OF DIRECT SUMMATION METHOD , CDDHF S - BASED DLCT

ξ

3

ξ

−10

10

2

ξ

1

−20

−20

10

10

−30

−30

10

10
0.5

1

1.5

2

2.5

3

0.5

α, σ or ξ
(a)

1

1.5

2

ξ1, ξ2 or ξ3

2.5

3

(b)

Fig. 3.
Normalized mean-square error (NMSE) of the approximation to
the sampled LCT versus the parameters of (a) CDDHFs-based DLCT and (b)
2
proposed DLCT. The input is sampled Gaussian function g1 [n] = e−2πn /N .
When one of the parameters varies from 0.1 to 3, the other two are fixed to
0.1 (σ to 1.1). (a) shows that the scaling parameter σ dominates the accuracy
of CDDHFs-based DLCT, while (b) shows that the accuracy of the proposed
method depends on the first and second chirp rates, i.e. ξ1 and ξ2 .
0

B. Accuracy of Approximation to Sampled LCT
In this subsection, we will examine the errors between the
DLCTs and the sampled output of the continuous LCT. In
[19], it has been mentioned that the continuous LCT of the
2
Gaussian function gs (t) = e−2πst is given by
r
ad−1+j2sbd 2 2
1
π u
(a,b;c,d)
{g
(t)}
=
e 2πsb2 −jπab
GM
(u)
=
O
,
s
s
LCT
a + j2bs
(80)
wherepb 6= 0. We adopt the sampling periods ∆t = ∆u =
∆ = 1/N which are suitable for both the proposed DLCT
and the CDDHFs-based DLCT [19]. Then, the sampled input
and sampled output of (80) are given by




k
n
M
M
√
, Gs [k] = Gs
,
(81)
gs [n] = gs √
N
N

10

−10

NMSE

Recall the matrix form of the DLCT based on direct
summation (11), i.e. X = Kx. It involves N 2 complex multiplications, and thus the computational complexity is 4N 2 real
multiplications. Next, consider the CDDHFs-based DLCT [19]
with matrix form shown in (10), i.e. X = Cξ Eσ Vα ET1 x. Each
of the real matrices E1 and Eσ yields 2N 2 real multiplications,
while each of the complex diagonal matrices Cξ and Vα leads
to N complex multiplications. Therefore, the computation
totally contains 4N 2 + 8N real multiplications. Regarding the
proposed DLCT, firstly consider the B 6= 0 case with matrix
form given in (21): X = C D−1 F† C−B FC A−1 x. There
B
B
are three discrete CMs which totally require 3N complex
multiplications. The DFT (F) and IDFT (F† ) can be computed
by the FFT along with (N/2) log2 N complex multiplications.
Accordingly, the complexity of the proposed DLCT for B 6= 0
is 12N + 4N log2 N real multiplications. The proposed DLCT
for B = 0 given in (26) and (27) contains one more DFT (or
IDFT) than the B 6= 0 case. It follows that the complexity
becomes 12N + 6N log2 N real multiplications. At last, we
summarize the comparison of the computational complexity
in TABLE I.

10

σ=1
σ=1.2
σ=0.7
σ=1.8

−20

10

−30

10

0

20

40

60

80

100

order p

Fig. 4.
Normalized mean-square error (NMSE) between the CDDHFs
(i.e.ΦCDDHF
) and the samples of scaled HFs for N = 101 and scaling
p;σ
parameter σ = 0.7, 1, 1.2, 1.8. There are 101 CDDHFs, i.e. ΦCDDHF
of order
p;σ
p = 0, 1, . . . , 100, which form an orthonormal set. As σ is farther from 1,
the error is larger.

respectively. The accuracy is measured by the normalized
mean-square error (NMSE) defined as
P M
2
M
Gs [k] − ODLCT
{gs [n]}
.
(82)
NMSE = k
P M 2
|Gs [k]|
k

Note that (A, B; C, D) in the proposed DLCT is equal
p to
(a, b; c, d) in the CDDHFs-based DLCT because ∆ = 1/N
and (58).
The CDDHFs-based DLCT can also be controlled by three
independent parameters: DFRFT angle α, scaling parameter
σ and chirp rate ξ, as shown in TABLE I. To examine the
influence of α on the accuracy, let α vary from 0.1 to 3 while
σ and ξ are fixed to 1.1 and 0.1, respectively. (If σ and ξ
are fixed to 1 and 0, respectively, b will become zero and not
suitable for (80).) Similarly, α and ξ are both fixed to 0.1 when
σ varies from 0.1 to 3. And α and σ are fixed to 0.1 and 1.1,
respectively, when ξ varies from 0.1 to 3. The NMSEs versus
α, σ and ξ are depicted in Fig. 3(a). Here, gs [n] with s = 1
and N = 101 is used as the input. We can find out that the
accuracy mainly depends on σ. It is because the approximation
errors of CDDHFs significantly increase as σ is farther from
1, as shown in Fig. 4. In contrast, the DFRFT angle α has very
little influence on the accuracy. It is because the DFRFT uses
only the undilated DHFs, i.e. CDDHFs with σ = 1, which
have good approximation. The chirp rate ξ in the last step

9

5

g0.2[n]

g [n]
1

5

5

gπ/2[n]

0

Approximation error (s = 0.2)

10

−10

0

0

−5
−5

−5
−5

−5
−5

10

−20

−20

10

5

0
(a)

5

G M1[k]
0.2

5

M2
G1 [k]

5

0

−5
−5

0
(b)

5

0

0
(d)

5

0
(c)

0
(e)

5

−5
−5

2

2

CDDHFs−based
Proposed
−30

10

G M3[k]
π/2

−30

50

100

150

200

10

50

(a)

100

150

200

(b)

Fig. 6. Normalized mean-square errors (NMSEs) of the approximation to
the sampled LCT for 200 different M’s. The NMSEs are sorted in ascending
order. The input signals are sampled Gaussian functions with (a) s = 0.2
i.e. g0.2 [n] and (b) s = π/2 i.e. gπ/2 [n]. The parameters in each M are
uniformly distributed random numbers on the interval (−2, 2).
0
(f)

5

Fig. 5.
The time-frequency distributions (TFDs) of the sampled inputs
(sampled Gaussian functions) gs [n] and the sampled LCT outputs GM
s [k].
In (a) and (d), s = 0.2 and M1 = (0, 4; −1/4, 2) are used. In (b) and (e),
s = 1 and M
 2 = (0.44, −0.08; 4.8,
 1.4) are used. In (c) and (f), s = π/2
0.8 √
√ , 7.65
√
and M3 = √
are used.
, 0.8 ; 5.15
2

10

5

0

−5
−5

Approximation error (s = π/2)

−10

10

NMSE

0

0

10

2

rameters in M are uniformly distributed random numbers on
the interval (−2, 2). And we sort the 200 data of NMSEs in
ascending order. The results using g0.2 [n] and gπ/2 [n] as the
inputs are depicted in Fig. 6(a) and (b), respectively. Generally,
the proposed method has somewhat higher accuracy than the
CDDHFs-based method.
C. Additivity Property

won’t introduce any error. The reason will be discussed later
together with the proposed DLCT.
Consider the proposed DLCT with B 6= 0 because
p (80) is
valid only when b 6= 0 (b = B because ∆ = 1/N ). As
shown in TABLE I, the proposed DLCT is actually controlled
by three independent chirp rates, i.e. ξ1 = A−1
B , ξ2 = −B and
ξ3 = D−1
B . Again, we examine the influence of each chirp rate
on the accuracy by fixing the other two to 0.1. The NMSEs
versus ξ1 , ξ2 ans ξ3 from 0.1 to 3 are shown in Fig. 3(b).
As mentioned in Sec. IV-B, CMs will produce shearing in
frequency domain while CC will produce shearing in time
domain. Therefore, larger values of ξ1 and ξ2 will lead to
larger bandwidths of x1 (t) and longer time duration of X1 (u),
cause larger errors in the two necessary conditions (64) and
(65) followed by lower accurate DLCT output X[k]. Although
the last CM with chirp rate ξ3 will also yield frequency domain
shearing, it won’t affect the accuracy of X[k] but the accuracy
of recovering X(u) from X[k].
Recall (63)-(67) for the proposed DLCT. The output X[k]
without aliasing effect doesn’t implies high accuracy of
X[k], and vice versa. This statement is also true for the
CDDHFs-based DLCT. Three examples are presented. The
time-frequency distributions (TFDs) of the sampled Gaussian
functions with s = 0.2, 1 and π/2 (i.e. g0.2 [n], g1 [n] and
gπ/2 [n]) are shown in Fig. 5(a)-(c). The TFDs of the sampled
LCT outputs with some parameter matrices M1 , M2 and M3
M2
M3
1
(i.e. GM
0.2 [k], G1 [k] and Gπ/2 [k]) are shown in Fig. 5(d)-(f).
M2
1
There’s no aliasing effect in GM
0.2 [k] and G1 [k]. However,
M1
the NMSE of G0.2 [k] calculated by the CDDHFs-based DLCT
2
is 1.9 × 10−4 . The NMSE of GM
1 [k] calculated by the
−4
3
proposed DLCT is up to 8.6 × 10 . In contrast, GM
π/2 [k] has
serious aliasing effect, but both the CDDHFs-based and the
proposed DLCTs have high accuracy – NMSEs below 10−15 .
At last, the proposed DLCT is compared with the CDDHFsbased DLCT by 200 simulation runs. In each run, the pa-

The additivity property for DLCTs is defined as
M1
M2
M1 ×M2
ODLCT
ODLCT
= ODLCT
.

(83)

In the following, we examine the NMSE of the additivity
property for the CDDHFs-based and proposed methods:
2
P M1 ×M2
M1
M2
ODLCT
{hi [n]} − ODLCT
ODLCT
{hi [n]}
NMSE = k
.
2
P M1 ×M2
ODLCT
{hi [n]}
k

(84)

Four kinds of discrete p
signals h1 [n], h2 [n], h3 [n] and h4 [n]
with sampling periods 1/N are shown in Fig. 7(a)-(d) and
described as follows:
− π n2 −j π n2
• h1 [n] = e
h N  N with
 N =128;

i
1 2
n
√n − 1
• h2 [n] = 2 cos 2π √
+
j
sin
π
e− N n
N
N
with N = 101;
• h3 [n] is a binary sequence with N = 280;
• h4 [n] is a trapezoidal-shaped function with N = 201.
The DFTs of these four signals are shown in Fig. 7(e)-(h),
respectively.
Again, let the parameters in M1 and M2 be uniformly
distributed random numbers on the interval (−2, 2), and obtain
the additivity NMSE in (84) for 200 simulation runs. The
NMSEs sorted in ascending order using h1 [n], h2 [n], h3 [n]
and h4 [n] as the input are plotted in Fig. 8(a)-(d), respectively.
These four examples reveal that the proposed DLCT has
performance similar to the CDDHFs-based DLCT in the
additivity property. Besides, if N is large enough so that
the energy is well concentrated around the origin of timefrequency plane, such as h1 [n], “approximate” additivity can
be achieved. When N is not large enough, one way to reduce
the NMSE of additivity is to limit the values of M1 and M2 .
Choosing all the chirp rates sufficiently small, approximate

10
h2[n]

h [n]
1

2

0.8
0.6
0.4
0.2
0
−0.2

h3[n]

1

h4[n]

1

1
0.5

0

0.5

−1
−5

0
time
(a)

5

−4

−2

DFT of h [n]

4

0

−5

Real
Imag

0.6
0.4

1

0.2

0.5

0

0

−0.2

−0.5
5

0
time
(c)

0

5

−5

DFT of h [n]

2

1.5

0
frequency
(e)

2

DFT of h [n]

1

−5

0
time
(b)

0
time
(d)

5

DFT of h [n]

3

4

4

4

3
2

2
1

0

0
−4

−2

0
2
frequency
(f)

4

−5

0
frequency
(g)

5

−5

0
frequency
(h)

5

p
Fig. 7. The time domain and frequency domain plots of four kinds of discrete signals h1 [n], h2 [n], h3 [n] and h4 [n] with sampling periods 1/N . Both
h1 [n] and h2 [n] are approximately time-limited and band-limited, and h1 [n] has energy more concentrated around the origin then h2 [n]. Both h3 [n] and
h4 [n] are time-limited but not band-limited, and h3 [n] has a larger bandwidth than h4 [n].
0

10

Additivity error for h1[n]

0

10

Additivity error for h2[n]

Additivity of proposed DLCT

0

10

−2

10

−4

10

−5

10

NMSE

NMSE

10
−10

CDDHFs−based
Proposed

−20

10

−6

10

−8

10

−10

50

100

150

200

10

50

(a)
0

10

100

150

200

(b)

Additivity error for h3[n]

0

10

ξ ∈[−2.5, 2.5]
ξ ∈[−1.5, 1.5]

−10

10

ξ ∈[−0.5, 0.5]

Additivity error for h4[n]

−12

10

20

40

60

80

100

120

140

160

180

200

NMSE

Simulation runs

−5

Fig. 9. Normalized mean-square errors (NMSEs) of the additivity property of
M1 ×M2
M2
M1
and ODLCT
, ODLCT
the proposed DLCT with all the chirp rates in ODLCT
are uniformly distributed random numbers within [−2.5, 2.5], [−1.5, 1.5] or
[−0.5, 0.5]. The NMSEs are obtained from 200 simulation runs and sorted
in ascending order.

−5

10

10
CDDHFs−based
Proposed

−10

10

−10

50

100

(c)

150

200

10

50

100

150

200

D. Reversibility Property

(d)

Fig. 8. Normalized mean-square errors (NMSEs) of the additivity property
for 200 different sets of M1 and M2 . The NMSEs are sorted in ascending
order. The input signals are (a) h1 [n], (b) h2 [n], (c) h3 [n] and (d) h4 [n]
depicted in Fig. 7. The parameters in M1 and M2 are uniformly distributed
random numbers on the interval (−2, 2).

Next, we examine the NMSE of the reversibility property:
2
P
M−1 M
hi [n] − ODLCT
ODLCT {hi [n]}
NMSE = n
.
(85)
P
|hi [n]|2
n

additivity will be achieved. For example, use h4 [n] as the
input. With all the chirp rates being uniformly distributed
random numbers within [−2.5, 2.5], [−1.5, 1.5] or [−0.5, 0.5],
the NMSEs of the additivity property from 200 simulation
runs are plotted in Fig. 9. It can be found that the NMSE of
ξ ∈ [−0.5, 0.5] is about 105 ∼ 106 times smaller than the
NMSE of ξ ∈ [−2.5, 2.5].

Again, let the parameters in M be uniformly distributed
random numbers on the interval (−2, 2). The NMSEs resulting
from 200 simulation runs are sorted in ascending order and
displayed in Fig. 10. The CDDHFs-based method doesn’t
satisfy the reversibility property perfectly. Although the proposed DLCT doesn’t has perfect additivity, it satisfies the
reversibility property perfectly. In Fig. 10, all the NMSEs of
the proposed DLCT are below 10−25 and numerically verify
the proofs. With the reversibility property, it is unnecessary
to develop the inverse DLCT additionally because it can be
realized by the forward DLCT with M−1 .
At the end, comparisons between the CDDHFs-based DLCT
and the proposed DLCT are summarized in TABLE II.

11

0

10

Reversibility error for h1[n]

−10

TABLE II

Reversibility error for h2[n]

C OMPARISONS BETWEEN THE CDDHF S - BASED DLCT [19] AND
PROPOSED DLCT

−10

10

NMSE

0

10

10
CDDHFs−based
Proposed

−20

Sampling periods

−20

10

10

Complexity
−30

10

0

10

−30

50

100

150

200

100

150

(a)

(b)
Reversibility error for h4[n]

0

10

−10

NMSE

50

Reversibility error for h3[n]

200

−10

10

10

CDDHFs−based −20
10
Proposed

−20

10

−30

10

10

−30

50

100

150

200

10

(c)

50

100

150

200

(d)

Fig. 10.
Normalized mean-square errors (NMSEs) of the reversibility
property for 200 different M’s. The NMSEs are sorted in ascending order.
The input signals are (a) h1 [n], (b) h2 [n], (c) h3 [n] and (d) h4 [n] depicted in
Fig. 7. The parameters in each M are uniformly distributed random numbers
on the interval (−2, 2).

VI. C ONCLUSION
In this paper, we develop a discrete LCT (DLCT) which
is irrelevant to the sampling periods and doesn’t involve oversampling operation. The proposed DLCT is based on the wellknown CM-CC-CM decomposition, which decomposes the
LCT to two chirp multiplications (CMs) and one chirp convolution (CC). One advantage of this decomposition over many
other decompositions is no scaling operation involved because
scaling operation will change the sampling period or introduce
interpolation error. The CM-CC-CM decomposition is invalid
for B = 0. Accordingly, we modify the decomposition and the
proposed DLCT fit for the B = 0 case. We also investigate
special cases of the proposed DLCT. The proposed DLCT can
be implemented by three discrete CMs and two FFTs (three
for B = 0), which yield lower computational complexity than
the previous works, DLCT calculated by direct summation
and DLCT based on center discrete dilated Hermite functions
(CDDHFs) [19]. The relation between the proposed DLCT
and the continuous LCT is also derived to approximate the
samples of the continuous LCT. Compared with the CDDHFsbased method, the proposed method has somewhat higher
approximation accuracy. Besides, simulation results show that
approximate additivity property can be achieved with error as
small as the CDDHFs-based method. Most importantly, the
proposed method has perfect reversibility, which is proved
mathematically and by numerical examples. With the reversibility property, the inverse transform of the proposed
DLCT can be realized by the forward DLCT.
R EFERENCES
[1] S. A. Collins, Jr., “Lens-system diffraction integral written in terms of
matrix optics,” JOSA, vol. 60, no. 9, pp. 1168–1177, 1970.

CDDHFs-based
p
∆x = ∆u = 1/N
O(N 2 )

Proposed
∆x = ∆u
O(N log2 N )

Accuracy

Worse

Better

Additivity

Approximate

Approximate

Reversibility

Approximate

Perfect

[2] M. Moshinsky and C. Quesne, “Linear canonical transformations and
their unitary representations,” Journal of Mathematical Physics, vol. 12,
p. 1772, 1971.
[3] H. M. Ozaktas, M. A. Kutay, and Z. Zalevsky, The fractional Fourier
transform with applications in optics and signal processing. New York:
Wiley, 2001.
[4] J.-J. Ding, “Research of fractional Fourier transform and linear canonical
transform,” Ph.D. dissertation, Ph. D. Thesis, National Taiwan University, 2001.
[5] S. C. Pei and J.-J. Ding, “Eigenfunctions of linear canonical transform,”
IEEE Transactions on Signal Processing, vol. 50, no. 1, pp. 11–26, 2002.
[6] M. Nazarathy and J. Shamir, “First-order opticsa canonical operator
representation: lossless systems,” JOSA, vol. 72, no. 3, pp. 356–364,
1982.
[7] M. J. Bastiaans, “Propagation laws for the second-order moments of
the wigner distribution function in first-order optical systems,” Optik,
vol. 82, no. 4, pp. 173–181, 1989.
[8] B. Barshan, M. A. Kutay, and H. M. Ozaktas, “Optimal filtering with
linear canonical transformations,” Optics communications, vol. 135, no.
1-3, pp. 32–36, 1997.
[9] S. C. Pei and J.-J. Ding, “Simplified fractional Fourier transforms,” JOSA
A, vol. 17, no. 12, pp. 2355–2367, 2000.
[10] S. C. Pei and J.-J. Ding, “Relations between fractional operations and
time-frequency distributions, and their applications,” IEEE Transactions
on Signal Processing, vol. 49, no. 8, pp. 1638–1655, 2001.
[11] M. J. Bastiaans and K. B. Wolf, “Phase reconstruction from intensity
measurements in linear systems,” JOSA A, vol. 20, no. 6, pp. 1046–1049,
2003.
[12] B. M. Hennelly and J. T. Sheridan, “Optical encryption and the space
bandwidth product,” Optics communications, vol. 247, no. 4, pp. 291–
305, 2005.
[13] K. K. Sharma and S. D. Joshi, “Signal separation using linear canonical
and fractional Fourier transforms,” Optics communications, vol. 265,
no. 2, pp. 454–460, 2006.
[14] S. C. Pei and S.-G. Huang, “Reversible joint hilbert and linear canonical
transform without distortion,” IEEE transactions on signal processing,
vol. 61, no. 17-20, pp. 4768–4781, 2013.
[15] S. Abe and J. T. Sheridan, “Optical operations on wave functions as the
abelian subgroups of the special affine Fourier transformation,” Optics
letters, vol. 19, no. 22, pp. 1801–1803, 1994.
[16] S. Abe and J. T. Sheridan, “Almost-Fourier and almost-Fresnel transformations,” Optics communications, vol. 113, no. 4, pp. 385–388, 1995.
[17] D. F. James and G. S. Agarwal, “The generalized Fresnel transform and
its application to optics,” Optics communications, vol. 126, no. 4-6, pp.
207–212, 1996.
[18] L. M. Bernardo, “Abcd matrix formalism of fractional Fourier optics,”
Optical Engineering, vol. 35, p. 732, 1996.
[19] S. C. Pei and Y.-C. Lai, “Discrete linear canonical transforms based on
dilated hermite functions,” JOSA A, vol. 28, no. 8, pp. 1695–1708, 2011.
[20] L. Zhao, J. J. Healy, and J. T. Sheridan, “Unitary discrete linear canonical
transform: analysis and application,” Applied optics, vol. 52, no. 7, pp.
C30–C36, 2013.
[21] A. Papoulis, Signal analysis. McGraw-Hill, 1977, vol. 191.
[22] B. M. Hennelly and J. T. Sheridan, “Fast numerical algorithm for the
linear canonical transform,” JOSA A, vol. 22, no. 5, pp. 928–937, 2005.
[23] H. M. Ozaktas, A. Koç, I. Sari, and M. A. Kutay, “Efficient computation
of quadratic-phase integrals in optics,” Optics letters, vol. 31, no. 1, pp.
35–37, 2006.
[24] A. Koç, H. M. Ozaktas, C. Candan, and M. A. Kutay, “Digital computation of linear canonical transforms,” IEEE Transactions on Signal
Processing, vol. 56, no. 6, pp. 2383–2394, 2008.
[25] K. B. Wolf, Geometric optics on phase space. Springer, 2004.

12

[26] R. Simon and N. Mukunda, “Iwasawa decomposition in first-order
optics: universal treatment of shape-invariant propagation for coherent
and partially coherent beams,” JOSA A, vol. 15, no. 8, pp. 2146–2155,
1998.
[27] S.-C. Pei and Y.-C. Lai, “Signal scaling by centered discrete dilated
hermite functions,” IEEE Transactions on Signal Processing, vol. 60,
no. 1, pp. 498–503, 2012.
[28] H. M. Ozaktas, O. Arikan, M. A. Kutay, and G. Bozdagt, “Digital
computation of the fractional Fourier transform,” Signal Processing,
IEEE Transactions on, vol. 44, no. 9, pp. 2141–2150, 1996.
[29] A. Koç, H. M. Ozaktas, and L. Hesselink, “Fast and accurate algorithms
for quadratic phase integrals in optics and signal processing,” in SPIE
Defense, Security, and Sensing. International Society for Optics and
Photonics, 2011, pp. 804 304–804 304.
[30] H. M. Ozaktas and A. Koç, “Fast and accurate linear canonical transform
algorithms,” in Signal Processing and Communications Applications
Conference (SIU), 2015 23th. IEEE, 2015, pp. 1409–1412.
[31] S.-C. Pei, S.-G. Huang, and J.-J. Ding, “Discrete gyrator transforms:
Computational algorithms and applications,” Signal Processing, IEEE
Transactions on, vol. 63, no. 16, pp. 4207–4222, Aug 2015.
[32] A. Papoulis, The Fourier integral and its applications. McGraw-Hill,
1962.
[33] V. Namias, “The fractional order Fourier transform and its application to
quantum mechanics,” J. Inst. Math. Appl., vol. 25, no. 3, pp. 241–265,
1980.
[34] A. Stern, “Sampling of linear canonical transformed signals,” Signal
Processing, vol. 86, no. 7, pp. 1421–1425, 2006.
[35] B.-Z. Li, R. Tao, and Y. Wang, “New sampling formulae related to linear
canonical transform,” Signal Processing, vol. 87, no. 5, pp. 983–990,
2007.
[36] J. J. Healy and J. T. Sheridan, “Sampling and discretization of the linear
canonical transform,” Signal Processing, vol. 89, no. 4, pp. 641–648,
2009.

Soo-Chang Pei (SM’89-F’00-LF’15) was born in
Soo-Auo, Taiwan, China on February 20, 1949. He
received the B. S. degree from National Taiwan University in 1970 and the M. S. and Ph. D. degree from
the University of California, Santa Barbara in 1972
and 1975 respectively, all in electrical engineering.
He was an engineering officer in the Chinese
Navy Shipyard from 1970 to 1971. From 1971 to
1975, he was a research assistant at the University of
California, Santa Barbara. He was the Professor and
Chairman in the EE department of Tatung Institute
of Technology and National Taiwan University, from 1981 to 1983 and 1995 to
1998, respectively. Presently, he is the Professor of EE department at National
Taiwan University. His research interests include digital signal processing, image processing, optical information processing, and laser holography. Dr. Pei
received National Sun Yet- Sen Academic Achievement Award in Engineering
in 1984, the Distinguished Research Award from the National Science Council
from 1990-1998, outstanding Electrical Engineering Professor Award from
the Chinese Institute of Electrical Engineering in 1998, and the Academic
Achievement Award in Engineering from the Ministry of Education in 1998,
the IEEE Fellow in 2000 for contributions to the development of digital
eigenfilter design, color image coding and signal compression, and to electrical
engineering education in Taiwan, the Pan Wen-Yuan Distinguished Research
Award in 2002, and the National Chair Professor Award from Ministry of
Education in 2002 and 2008. The IEEE Life Fellow in 2015 for recognition
of the years of royal membership and support of the activities of IEEE. He
has been President of the Chinese Image Processing and Pattern Recognition
Society in Taiwan from 1996-1998.
Dr. Pei is IEEE Life Fellow and a member of Eta Keppa Nu and the Optical
Society of America.

Shih-Gu Huang was born in Taiwan in 1984. He
received the B.S. degree in electrical engineering
and the M.S. degree in communications engineering
from National Tsing Hua University, Hsinchu, Taiwan, in 2007 and 2009, respectively. He is currently
working toward the Ph.D. degree in the Graduate
Institute of Communication Engineering, National
Taiwan University, Taipei, Taiwan. His research
interests include digital signal processing, timefrequency analysis, fractional Fourier transform, and
linear canonical transform.

