Discrete linear canonical transform computation by adaptive method

Feng Zhang; Ran Tao; Yue Wang

doi:10.1364/OE.21.018138

1. Introduction

The linear canonical transform (LCT) is a three-parameter linear integral transform that describes the effect of quadratic phase systems on a wavefield [1,2]. The LCT generalizes many optical transforms, such as the Fourier transform (FT), the fractional Fourier transform (FRFT) [1], the Fresnel transform [3,4], the scaling operation. These integral transforms are of great importance in electromagnetic, acoustic, and other wave propagation problems where they represent the solution of the wave equation under a variety of circumstances [5]. In optics, the LCT can describe a broad class of optical systems including thin lenses, sections of free space in the Fresnel approximation, section of quadratic graded-index media, and linear quadratic phase systems [1]. The LCT itself has found many applications in phase reconstruction, filter design, graded-index media analysis, signal synthesis, pattern recognition, time-frequency analysis, optimal filtering, holographic three-dimensional television and many others [5–8]. Understanding the LCT can help to gain more insight into its special cases and to carry the knowledge gained from one subject to others [3].

Discrete counterpart of a continuous transform is very important for computing the samples of the continuous transform. The well-known discrete Fourier transform (DFT) is commonly used to obtain the samples of the FT, and is thus widely applied in signal and image processing [9,10]. Therefore, after the continuous LCT has been introduced, the definition and the implementation of the discrete LCT have been derived [5,11,12]. These methods are used to approximate the continuous LCT in the same sense that the DFT approximates the continuous FT.

However, the above works to compute the discrete LCT are generally limited to implementing it on a data block and thus not to compute the discrete LCT of the stream data which updates with every new sample. Besides, due to the floating point arithmetic in software and hardware realizations, noise, for example come from roundoff errors, can propagate and accumulate over time.

In this paper, we present the discrete LCT computation by the adaptive method. The approaches for computing the block discrete LCT and the sliding discrete LCT via the LMS algorithm are derived, and the implementation structures of these approaches by the adaptive filter system are considered. Since the adaptive algorithm can automatically adjust for possible errors, the appearing error in the computation process of the discrete LCT is attenuated when it propagates over time. Furthermore, the proposed computation approaches for the block discrete LCT and the sliding discrete LCT have the inherent parallel structures and is suitable for efficient very large scale integration (VLSI) implementations, which implies that the processing speed is essentially independent of the size of the data block $N$ .

2. Linear canonical transform and LMS adaptive filter

2.1 Linear canonical transform

The LCT of a signal $f (t)$ with parameters $(a, b, c, d)$ is defined as [4]

f_{M} (u) = L^{M} [f] (u) = {\begin{matrix} \sqrt{\frac{1}{j 2 π b}} \cdot \exp (j \frac{d}{2 b} u^{2}) \int_{- \infty}^{\infty} f (t) \exp (j \frac{a}{2 b} t^{2} - j \frac{1}{b} u t) d t, \begin{matrix} b \neq 0 \end{matrix} \\ \begin{matrix} \sqrt{d} \exp (j \frac{c d}{2} u^{2}) f (d u), & \begin{matrix}  \end{matrix} & b = 0 \end{matrix} \end{matrix}

where

M = [\begin{matrix} a & b; & c & d \end{matrix}]

is the parameter matrix and

\det (M) = a d - b c = 1

. Two successive LCTs with matrices

M_{1}

and

M_{2}

is another LCT with the matrix

M_{3} = M_{2} M_{1}

, and consequently the inverse LCT is given by the LCT with parameters

(d, - b, - c, a)

. Note that when

b = 0

the LCT becomes a chirp multiplication. Therefore, we only consider the case of

b \neq 0

, and without loss of generality we assume

b > 0

in this paper. When

[\begin{matrix} a & b; & c & d \end{matrix}] = [\begin{matrix} \cos α & \sin α; & - \sin α & \cos α \end{matrix}]

, the LCT reduces to the FRFT. For further details about properties of the LCT, readers can refer to [1–7,11–14].

The LCT of the ideal uniform sampled signal with sampling interval $Δ t$ is [13]

\begin{matrix} {\tilde{f}}_{M} (u) = L^{M} [\sum_{n = - \infty}^{\infty} f (n Δ t) δ (t - n Δ t)] (u) \\ = \sqrt{\frac{1}{j 2 π b}} \cdot \exp (j \frac{d}{2 b} u^{2}) \sum_{n = - \infty}^{\infty} f (n Δ t) \exp (j \frac{a}{2 b} n^{2} Δ t^{2} - j \frac{1}{b} u n Δ t) . \\ = \frac{1}{Δ t} \exp (j \frac{d}{2 b} u^{2}) \sum_{n = - \infty}^{\infty} f_{M} (u - n \frac{2 π b}{Δ t}) \exp [- j \frac{d}{2 b} {(u - n \frac{2 π b}{Δ t})}^{2}] \end{matrix}

Equation (2) shows that the LCT of the uniform sampled signal replicates

f_{M} (u)

with a period of

2 π b / Δ t

.

By sampling uniformly $N$ points in the replicated period $2 π b / Δ t$ in the LCT domain, the sampling interval $Δ u$ in the LCT domain is $Δ u = 2 π b / (N Δ t)$ . Then the discrete LCT of $f [n]$ can been defined as [11,12]

f_{M} [m] = \sqrt{1 / N} \exp (j \frac{d}{2 b} m^{2} Δ u^{2}) \sum_{n \in 〈 N 〉} f [n] \exp (j \frac{a}{2 b} n^{2} Δ t^{2} - j \frac{2 π}{N} m n)

where

〈 N 〉

denotes the interval

[0, 1, \dots, N - 1]

without generality. The conventional time-bandwidth product

(1 / Δ f) \cdot (1 / Δ t) = (2 π / Δ Ω) \cdot (1 / Δ t) = N

is the minimum number of samples to identify a signal whose energy are confined to

(1 / Δ f)

and

(1 / Δ t)

in the time domain and the frequency domain respectively. Likewise, the product

(2 π b / Δ u) \cdot (b / Δ t) \cdot (1 / b) = N

is the minimum number of samples to identify a signal whose energy are confined to

(2 π b / Δ u)

and

(b / Δ t)

in the time domain and the LCT domain respectively [5,12]. This minimum sample number

N

is also referred as the number of freedom degrees.

2.2 LMS adaptive filter

Adaptive filters are widely used in a variety of applications, including noise cancellation, system modeling, linear prediction filtering, and spectral estimation [15]. An adaptive filter has an adaptation process that automatically seeks an optimal impulse response by adjustable weights [16]. The model of the adaptive filter is shown in Fig. 1, where $X [n]$ is the input signal vector, $W [n]$ is the weight vector on the $n th$ iteration of the algorithm, $d [n]$ is the desired response which is also named reference signal, $y [n] = X^{T} [n] W [n]$ is the estimate output, and $e [n] = d [n] - y [n]$ is the error signal.

Fig. 1 The model of the adaptive filter.

Download Full Size | PDF

The objective of the adaptation process is to minimize the average total error power [15,16]

\begin{matrix} ξ = E {{| e [n] |}^{2}} = E {e [n] e^{*} [n]} \\ = E {(d [n] - X^{T} [n] W [n]) (d^{*} [n] - X^{H} [n] W^{*} [n])} \end{matrix}

where

E {\cdot}

denotes statistical expectation. It is shown that the performance function

ξ

(also called the cost function) is a convex quadratic function of weights.

The steepest descent method uses gradients of the performance surface for seeking its minimum, where each change in the weight vector proportional to the negative of the gradient vector

W [n + 1] = W [n] + μ (- \nabla [n])

where

μ

is a step parameter.

The most popular algorithm in adaptive signal processing is the least-mean-square (LMS) algorithm [16]. This algorithm is an implementation of steepest descent using the estimated gradients $\hat{\nabla} [n]$ , which takes the instantaneous gradients of $\hat{ξ} = {| e [n] |}^{2} = e [n] e^{*} [n]$ . The algorithm is given by

\begin{matrix} W [n + 1] = W [n] + μ (- \hat{\nabla} [n]) \\ = W [n] + 2 μ e [n] X^{*} [n] \\ = (I - 2 μ X^{*} [n] X^{T} [n]) W [n] + 2 μ d [n] X^{*} [n] \end{matrix}

The LMS adaptive algorithm has been applied to adaptive process modeling, adaptive noise canceling, adaptive line enhancer, etc [15–18].

3. Block discrete LCT computation by adaptive filter

In [19], the relationship between the DFT and the least-mean-square (LMS) algorithm has been established, which is used to compute the DFT via the adaptive LMS algorithm. Due to the adaptive algorithm adjustment, this computation method for the DFT has a built-in error cancellation mechanism. Moreover, this technique has the inherent parallel structure and is very suitable for efficient VLSI implementation, which implies that the processing speed is independent of the size of the data block $N$ .

Discrete LCT computation is the basis for applications of the continuous LCT, and we propose the discrete LCT computation using the adaptive LMS algorithm. In this section, we present the block-based discrete LCT computation by an adaptive filter, which has the parallel structure that is suitable for efficient VLSI implementation.

In order to implement the computation of the discrete LCT by the adaptive LMS algorithm, the discrete signal to be transformed is served as the desired response $d [n]$ of the adaptive filter system, and the input signal vector $X [n]$ to the system at time index $n$ is constructed as

\begin{array}{l} X [n] = \frac{1}{\sqrt{N}} \exp (- j \frac{a}{2 b} n^{2} Δ t^{2}) \times [\begin{matrix} 1 & \exp (- j \frac{d}{2 b} \cdot 1^{2} \cdot Δ u^{2} + j \frac{2 π}{N} \cdot 1 \cdot n) \end{matrix} \\ \begin{matrix}  \end{matrix} {\begin{matrix} \exp (- j \frac{d}{2 b} \cdot 2^{2} \cdot Δ u^{2} + j \frac{2 π}{N} \cdot 2 \cdot n) & \dots & \exp (- j \frac{d}{2 b} \cdot {(N - 1)}^{2} \cdot Δ u^{2} + j \frac{2 π}{N} \cdot (N - 1) \cdot n) \end{matrix}]}^{T} \end{array}

The choice of

X [n]

can be motivated as follows. The result of the continuous adaptive process is the weight vector

W [n]

, which is expected to be associated with the discrete LCT of

d [n]

. Note that the discrete LCT can be regarded as the expansion coefficient of transformed signal in term of the harmonic series with quadratic phase terms, which can be seen from the discrete transform kernel. This is the generalized results of the FT case that the DFT of a discrete signal is the expansion coefficient in terms of the discrete harmonic series

\exp (j \cdot k \frac{2 π}{N} n)

to the discrete signal. Therefore, by choosing the vector

X [n]

in (7), the weighted sum of harmonic series with quadratic phase terms can well match the desired signal

d [n]

, which can transform

d [n]

into its discrete LCT by means of the adaptive LMS algorithm.

For an adaptive filter, an adaptation process can automatically seek an optimal impulse response by adjusting the weight vector $W [n]$ . In addition to the input signal vector $X [n]$ , the desired response $d [n]$ is supplied to the adaptive filter during the adaptation process. In many applications, considerable ingenuity on $d [n]$ is required to obtain a better performance for an adaptive process. For example in the application of system modeling, the desired response $d [n]$ is the output of the unknown system to be modeled; in the application of interference canceling, the desired response $d [n]$ is the delayed version of the input signal; in the application of channel equalization, the desired response $d [n]$ is the output of the decision device. Compared to these applications, the adaptive filter used here looks somewhat special that rather than the output $y [n]$ the weight vector $W [n]$ is the cared result that gives the “expected” discrete LCT transform with certain parameters $(a, b, c, d)$ by means of the LMS algorithm. It should be pointed that the proposed method is used to compute the discrete LCT with certain known parameters $(a, b, c, d)$ which gives the certain input vector $X [n]$ to the adaptive filter.

The filter weight vector $W [n]$ is adjusted with the LMS algorithm which is given in (6). Let the initial weight vector $W [0]$ be set to zero. We obtain

W [1] = (I - 2 μ X^{*} [0] X^{T} [0]) W [0] + 2 μ d [0] X^{*} [0] = 2 μ d [0] X^{*} [0],

and

\begin{matrix} W [2] = (I - 2 μ X^{*} [1] X^{T} [1]) W [1] + 2 μ d [1] X^{*} [1] \\ = - 4 μ^{2} d [0] X^{*} [1] X^{T} [1] X^{*} [0] + 2 μ (d [0] X^{*} [0] + d [1] X^{*} [1]) . \end{matrix}

Using the relation

X^{T} [1] X^{*} [0] = \frac{1}{N} \exp [- j \frac{a}{2 b} (1^{2} - 0^{2}) Δ t^{2}] \sum_{k = 0}^{N - 1} \exp [j \frac{2 π}{N} \cdot (1 - 0) \cdot k] = 0,

we have

W [2] = 2 μ (d [0] X^{*} [0] + d [1] X^{*} [1]) .

Following the same procedure, we can obtain the weight vector

W [n] = 2 μ \sum_{m = 0}^{n - 1} d [m] X^{*} [m], n = 1, \dots, N .

For the case that

n = N

, it gives

W [N] = \frac{2 μ}{\sqrt{N}} \cdot [\begin{matrix} \sum_{m = 0}^{N - 1} d [m] \exp (j \frac{a}{2 b} \cdot m^{2} \cdot Δ t^{2}) \\ \exp (j \frac{d}{2 b} \cdot 1^{2} \cdot Δ u^{2}) \sum_{m = 0}^{N - 1} d [m] \exp (j \frac{a}{2 b} \cdot m^{2} \cdot Δ t^{2} - j \frac{2 π}{N} \cdot 1 \cdot m) \\ \exp (j \frac{d}{2 b} \cdot 2^{2} \cdot Δ u^{2}) \sum_{m = 0}^{N - 1} d [m] \exp (j \frac{a}{2 b} \cdot m^{2} \cdot Δ t^{2} - j \frac{2 π}{N} \cdot 2 \cdot m) \\ ⋮ \\ \exp (j \frac{d}{2 b} \cdot {(N - 1)}^{2} \cdot Δ u^{2}) \sum_{m = 0}^{N - 1} d [m] \exp (j \frac{a}{2 b} \cdot m^{2} \cdot Δ t^{2} - j \frac{2 π}{N} \cdot (N - 1) \cdot m) \end{matrix}]

We can clearly see that the element of

W [N]

comprise the values of the discrete LCT of the signal block

{[d [0], d [1], \dots, d [N - 1]]}^{T}

, except for a scale factor of

2 μ

.

Equation (12) is only applicable for $W [n]$ with $n = 1, \dots, N$ . Thus, we need to compute $W [N + 1]$ again

\begin{matrix} W [N + 1] = (I - 2 μ X^{*} [N] X^{T} [N]) W [N] + 2 μ d [N] X^{*} [N] \\ = 2 μ \sum_{m = 0}^{N} d [m] X^{*} [m] - 4 μ^{2} X^{*} [N] (X^{T} [N] \sum_{m = 0}^{N - 1} d [m] X^{*} [m]) . \end{matrix}

Note that for the input signal vector

X [n + N] \neq X [n]

. However, the products

X^{T} [N] X^{*} [l]

for

l = 0, 1, \dots, (N - 1)

are

X^{T} [N] X^{*} [l] = {\begin{matrix} \exp (- j \frac{a}{2 b} N^{2} Δ t^{2}), \begin{matrix} l = 0 \end{matrix} \\ 0, \begin{matrix} l = 1, \dots, (N - 1) \end{matrix} \end{matrix} .

Substituting (15) into (14), we can obtain

W [N + 1] = 2 μ \sum_{m = 1}^{N} d [m] X^{*} [m] - 4 μ^{2} X^{*} [N] d [0] \exp (- j \frac{a}{2 b} N^{2} Δ t^{2}) + 2 μ d [0] X^{*} [0] .

Combining the relationship between

X [N]

and

X [0]

that

X [N] = X [0] \exp (- j \frac{a}{2 b} N^{2} Δ t^{2})

, we can obtain

W [N + 1] = 2 μ \sum_{m = 1}^{N} d [m] X^{*} [m] + 2 μ (1 - 2 μ) d [0] X^{*} [0] .

Similarly, using the relationship

X^{T} (N + 1) = X^{T} (1) \exp {- j \frac{a}{2 b} [{(N + 1)}^{2} - 1^{2}] Δ t^{2}}

and the product

X^{T} [N + 1] X^{*} [1] = \exp {- j \frac{a}{2 b} [{(N + 1)}^{2} - 1^{2}] Δ t^{2}}

, we can obtain the weight vector at the time

(N + 2)

W [N + 2] = 2 μ \sum_{m = 2}^{N + 1} d [m] X^{*} [m] + 2 μ (1 - 2 μ) (X^{*} [1] d [1] + X^{*} [0] d [0]) .

Along this line, we can obtain

W [n] = 2 μ \sum_{m = n - N}^{n - 1} d [m] X^{*} [m] + 2 μ (1 - 2 μ) \sum_{m = 0}^{n - N - 1} d [m] X^{*} [m], n = (N + 1), \dots, 2 N .

Moreover, we can obtain a completely general formula for

n \geq 1

with the initial condition

W (0) = 0

. That is

\begin{array}{l} W [n] = 2 μ \sum_{m = n - N}^{n - 1} d [m] X^{*} [m] + 2 μ (1 - 2 μ) \sum_{m = n - 2 N}^{n - N - 1} d [m] X^{*} [m] + \\ 2 μ)^{2} \sum_{m = n - 3 N}^{n - 2 N - 1} d [m] X^{*} [m] + \dots . \end{array}

Note that the objective using the LMS adaptive algorithm here is to compute the discrete LCT. Except the first term, every term of the series (20) has a multiplicative factor of $(1 - 2 μ)$ before the summation. Therefore, if the step size is chosen as $μ = 1 / 2$ then the series reduces to its first term only, regardless of the time index $n$ . In fact, with this particular choice of the step size that $μ = 1 / 2$ , the discrete LCT can be computed by means of the adaptive LMS algorithm which is shown as follows. At the time $n = l N$ where $l$ is an integer the weight vector becomes

\begin{matrix} W [l N] = \sum_{m = l N - N}^{l N - 1} d [m] X^{*} [m] \\ = \frac{1}{\sqrt{N}} [\begin{matrix} \sum_{m = l N - N}^{l N - 1} d [m] \exp (j \frac{a}{2 b} \cdot m^{2} \cdot Δ t^{2}) \\ \exp (j \frac{d}{2 b} \cdot 1^{2} \cdot Δ u^{2}) \sum_{m = l N - N}^{l N - 1} d [m] \exp (j \frac{a}{2 b} \cdot m^{2} \cdot Δ t^{2} - j \frac{2 π}{N} \cdot 1 \cdot m) \\ \exp (j \frac{d}{2 b} \cdot 2^{2} \cdot Δ u^{2}) \sum_{m = l N - N}^{l N - 1} d [m] \exp (j \frac{a}{2 b} \cdot m^{2} \cdot Δ t^{2} - j \frac{2 π}{N} \cdot 2 \cdot m) \\ ⋮ \\ \exp (j \frac{d}{2 b} \cdot {(N - 1)}^{2} \cdot Δ u^{2}) \sum_{m = l N - N}^{l N - 1} d [m] \exp (j \frac{a}{2 b} \cdot m^{2} \cdot Δ t^{2} - j \frac{2 π}{N} \cdot (N - 1) \cdot m) \end{matrix}] . \end{matrix}

Although Eq. (21) is very similar to the discrete LCT of the input signal block

{[d [l N - N], d [l N - N + 1], \dots, d [l N - 1]]}^{T}

at the first look, they are different indeed. That is to say we can’t use (21) to compute the block discrete LCT directly. From (21) the weights

w_{k} [l N]

are

\begin{matrix} w_{k} [l N] = \frac{1}{\sqrt{N}} \exp (j \frac{d}{2 b} \cdot k^{2} \cdot Δ u^{2}) \sum_{m = l N - N}^{l N - 1} d [m] \exp (j \frac{a}{2 b} \cdot m^{2} \cdot Δ t^{2} - j \frac{2 π}{N} \cdot m \cdot k) \\ = \frac{1}{\sqrt{N}} \exp (j \frac{d}{2 b} \cdot k^{2} \cdot Δ u^{2}) \sum_{m = l N - N}^{l N - 1} d [m] \exp {j \frac{a}{2 b} \cdot {[m - (l N - N)]}^{2} \cdot Δ t^{2} - j \frac{2 π}{N} \cdot m \cdot k} \times \\ \begin{matrix}  \end{matrix} \exp {j \frac{a}{2 b} \cdot [- {(l N - N)}^{2} + 2 m (l N - N)] \cdot Δ t^{2}} \\ = \frac{1}{\sqrt{N}} \exp (j \frac{d}{2 b} \cdot k^{2} \cdot Δ u^{2}) \sum_{m = 0}^{N - 1} d [m + l N - N] \exp (j \frac{a}{2 b} \cdot m^{2} \cdot Δ t^{2} - j \frac{2 π}{N} \cdot m \cdot k) \times \\ \begin{matrix}  \end{matrix} \exp {j \frac{a}{2 b} \cdot [2 m (l N - N) + {(l N - N)}^{2}] \cdot Δ t^{2}} . \end{matrix}

Therefore, if we want to compute the block discrete LCT of the input signal block

{[d [l N - N], d [l N - N + 1], \dots, d [l N - 1]]}^{T}

, we need to multiply a phase factor

\exp {j \frac{a}{2 b} \cdot [- {(l N - N)}^{2} + 2 n (l N - N)] \cdot Δ t^{2}}

with

d [n]

from the time

n = l N - N

to

n = l N - 1

. Thus, the computation of the block discrete LCT by using the LMS adaptive filtering system is shown in Fig. 2. From Fig. 2 we can see that compared to directly compute the discrete LCT, the discrete LCT computation by the adaptive filtering technique has the inherent parallel structure that makes it suitable for efficient VLSI implementation. Moreover, this method can automatically adjust for possible errors as it propagates over time, which is not shared by other method for the computation of the discrete LCT.

Fig. 2 Block discrete LCT computation by LMS adaptive filter.

Download Full Size | PDF

4. Sliding discrete LCT computation by adaptive filter

We have shown the block discrete LCT computation by adaptive filter system. This method is suitable for implementing discrete LCT on a block by block basis. In some application scenarios, discrete time signals need to be transformed continuously in the LCT domain. At each instant, the $N$ most recent samples of the input signal are transformed in the LCT domain, which is referred as the sliding discrete LCT.

For the case when the signal is a stream of samples, the sliding discrete LCT vector is

\begin{array}{l} LCT [n] = [\begin{matrix} {LCT}_{0} [n] \\ {LCT}_{1} [n] \\ {LCT}_{2} [n] \\ ⋮ \\ {LCT}_{N - 1} [n] \end{matrix}] = \\ \frac{1}{\sqrt{N}} [\begin{matrix} \sum_{m = 0}^{N - 1} d [n - (N - 1) + m] \exp (j \frac{a}{2 b} \cdot m^{2} \cdot Δ t^{2}) \\ \exp (j \frac{d}{2 b} \cdot 1^{2} \cdot Δ u^{2}) \sum_{m = 0}^{N - 1} d [n - (N - 1) + m] \exp (j \frac{a}{2 b} \cdot m^{2} \cdot Δ t^{2} - j \frac{2 π}{N} \cdot 1 \cdot m) \\ \exp (j \frac{d}{2 b} \cdot 2^{2} \cdot Δ u^{2}) \sum_{m = 0}^{N - 1} d [n - (N - 1) + m] \exp (j \frac{a}{2 b} \cdot m^{2} \cdot Δ t^{2} - j \frac{2 π}{N} \cdot 2 \cdot m) \\ ⋮ \\ \exp (j \frac{d}{2 b} \cdot {(N - 1)}^{2} \cdot Δ u^{2}) \sum_{m = 0}^{N - 1} d [n - (N - 1) + m] \exp (j \frac{a}{2 b} \cdot m^{2} \cdot Δ t^{2} - j \frac{2 π}{N} \cdot (N - 1) \cdot m) \end{matrix}] . \end{array}

When

μ = 1 / 2

, from (20) we have

\begin{matrix} W [n] = \sum_{m = n - N}^{n - 1} d [m] X^{*} [m] \\ = \frac{1}{\sqrt{N}} [\begin{matrix} \sum_{m = n - N}^{n - 1} d [m] \exp (j \frac{a}{2 b} \cdot m^{2} \cdot Δ t^{2}) \\ \exp (j \frac{d}{2 b} \cdot 1^{2} \cdot Δ u^{2}) \sum_{m = n - N}^{n - 1} d [m] \exp (j \frac{a}{2 b} \cdot m^{2} \cdot Δ t^{2} - j \frac{2 π}{N} \cdot 1 \cdot m) \\ \exp (j \frac{d}{2 b} \cdot 2^{2} \cdot Δ u^{2}) \sum_{m = n - N}^{n - 1} d [m] \exp (j \frac{a}{2 b} \cdot m^{2} \cdot Δ t^{2} - j \frac{2 π}{N} \cdot 2 \cdot m) \\ ⋮ \\ \exp (j \frac{d}{2 b} \cdot {(N - 1)}^{2} \cdot Δ u^{2}) \sum_{m = n - N}^{n - 1} d [m] \exp (j \frac{a}{2 b} \cdot m^{2} \cdot Δ t^{2} - j \frac{2 π}{N} \cdot (N - 1) \cdot m) \end{matrix}] . \end{matrix}

The weights

w_{k} [n]

can be rewritten as

\begin{matrix} w_{k} [n] = \frac{1}{\sqrt{N}} \exp (j \frac{d}{2 b} \cdot k^{2} \cdot Δ u^{2}) \sum_{m = n - N}^{n - 1} d [m] \exp (j \frac{a}{2 b} \cdot m^{2} \cdot Δ t^{2} - j \frac{2 π}{N} \cdot m \cdot k) \\ = \frac{1}{\sqrt{N}} \exp (j \frac{d}{2 b} \cdot k^{2} \cdot Δ u^{2}) \sum_{m = 0}^{N - 1} d [m + n - N] \times \\ \begin{matrix} \begin{matrix}  \end{matrix} \end{matrix} \exp [j \frac{a}{2 b} \cdot {(m + n - N)}^{2} \cdot Δ t^{2} - j \frac{2 π}{N} \cdot (m + n - N) \cdot k] \\ = \frac{1}{\sqrt{N}} \exp (j \frac{d}{2 b} \cdot k^{2} \cdot Δ u^{2} - j \frac{2 π}{N} \cdot n \cdot k) \sum_{m = 0}^{N - 1} d [m + n - N] \times . \\ \begin{matrix}  \end{matrix} \exp [j \frac{a}{2 b} \cdot {(m + n - N)}^{2} \cdot Δ t^{2} - j \frac{2 π}{N} \cdot m \cdot k] \\ = \frac{1}{\sqrt{N}} \exp (j \frac{d}{2 b} \cdot k^{2} \cdot Δ u^{2} - j \frac{2 π}{N} \cdot n \cdot k + j \frac{a}{2 b} \cdot {(n - N)}^{2} \cdot Δ t^{2}) \sum_{m = 0}^{N - 1} d [m + n - N] \times \\ \begin{matrix}  \end{matrix} \exp {j \frac{a}{2 b} \cdot [m^{2} \cdot Δ t^{2} + 2 m (n - N) \cdot Δ t^{2}] - j \frac{2 π}{N} \cdot m \cdot k} \end{matrix}

From (23) we have

{LCT}_{k} [n - 1] = \frac{1}{\sqrt{N}} \exp (j \frac{d}{2 b} \cdot k^{2} \cdot Δ u^{2}) \sum_{m = 0}^{N - 1} d [n - N + m] \exp (j \frac{a}{2 b} \cdot m^{2} \cdot Δ t^{2} - j \frac{2 π}{N} \cdot k \cdot m) .

Comparing (25) and (26) we can see that in order to complete the sliding discrete LCT computation of

{LCT}_{k} [n - 1]

by weights

w_{k} [n]

, the weights

w_{k} [n]

at time index

n

need to multiply the corresponding phase factor

\exp [- j \frac{a}{2 b} \cdot {(n - N)}^{2} \cdot Δ t^{2} + j \frac{2 π}{N} \cdot n \cdot k]

and the input samples

d [n]

need to multiply a time-variant phase factor

\exp [j \frac{a}{2 b} \cdot 2 m (n - N) \cdot Δ t^{2}]

,

m = 0, 1, \dots, N - 1

, which can be complemented by a time-variant phase network.

Figure 3(a) shows the configuration of the computation of sliding discrete LCT by the adaptive method. Figure 3(b) shows the details of the time-variant phase network, where $z^{- D}$ denotes $D$ times delay, $D = 0, 1, \dots, N - 1$ .

Fig. 3 (a) Sliding discrete LCT computation by LMS adaptive filter and (b) details of time-variant phase network.

Download Full Size | PDF

5. Simulation results

In this Section, numerical examples of discrete LCT computation by the proposed adaptive LMS algorithm are demonstrated. The original continuous signal is a Gaussian function $f (t) = \exp [- {(t - t_{0})}^{2} / 2]$ where $t_{0} = 4$ . The discrete signal $f [n]$ is sampled from the continuous signal $f (t)$ with sampling interval $Δ t = 0.04$ , that is $f [n] = f (n Δ t)$ . The continuous signal $f (t)$ and the discrete signal $f [n]$ are shown in Fig. 4. Since $\exp (- 8) \approx 0$ , signals $f (t)$ and $f [n]$ are approximated as finite signals in the interval $[0, 8]$ , and thus can be numerically simulated in a computer.

Fig. 4 (a) The continuous signal and (b) its discrete version.

Download Full Size | PDF

First, we consider the discrete LCT of $f [n]$ with parameters $(a, b, c, d) = (3, 1, 1 / 2, 1 / 2)$ . The discrete LCTs of $f [n]$ computed by the definition (3) directly and by the proposed adaptive method are shown in Fig. 5. The amplitude and real part results of the discrete LCT are plotted in Fig. 5(a) and Fig. 5(b), respectively. Also, the continuous LCT of original continuous signal $f (t)$ is plotted in Fig. 5. From Fig. 5, we can see that the discrete LCT of $f [n]$ computed by the proposed adaptive method is exactly equal to the discrete LCT computed by the definition (3), which are well approximated as the continuous LCT of $f (t)$ .

Fig. 5 (a) The amplitude and (b) the real part results of the continuous LCT and the discrete LCT with parameters $(a, b, c, d) = (3, 1, 1 / 2, 1 / 2)$ .

Download Full Size | PDF

Then, we consider the discrete LCT of $f [n]$ with parameters $(a, b, c, d) = (\cos α, \sin α, - \sin α, \cos α)$ by $α = 3 π / 20$ , which is essentially the discrete FRFT of $f [n]$ . The discrete LCTs of $f [n]$ computed by the definition (3) directly and by the proposed adaptive method are shown in Fig. 6, where the amplitude and real part of the discrete LCT are plotted in Fig. 6(a) and Fig. 6(b), respectively. Also, the continuous LCT of original continuous signal $f (t)$ is plotted in Fig. 6. For a Gaussian function $\exp (- t^{2} / 2)$ , its FRFT has the analytic closed-form $\exp (- u^{2} / 2)$ . Then according to the shift property of the FRFT, the FRFT of signal $f (t) = \exp [- {(t - t_{0})}^{2} / 2]$ is $f_{M} (u) = \exp [- {(u - t_{0} \cos α)}^{2} / 2] \exp (j t_{0}^{2} \sin α \cos α / 2 - j u t_{0} \sin α)$ . The position of the maximum amplitude value of continuous FRFT is $t_{0} \cos α = 3.564$ theoretically. From Fig. 6 it is shown that the discrete LCT of $f [n]$ computed by the proposed adaptive method is exactly equal to the discrete LCT computed by the definition (3). Besides, from Fig. 6(a), the position of the maximum amplitude value of the discrete LCT is 3.566, which implies the well approximation to the continuous LCT.

Fig. 6 (a) The amplitude and (b) the real part results of the continuous LCT and the discrete LCT with parameters $(a, b, c, d) = (\cos α, \sin α, - \sin α, \cos α)$ by $α = 3 π / 20$ .

Download Full Size | PDF

6. Conclusion

In this paper, the computation of the block discrete LCT and the sliding discrete LCT by LMS adaptive filter system has been proposed. We have shown that with a particular choice of the step to be $μ = \frac{1}{2}$ , the discrete LCT of a signal block with length $N$ can be computed by means of the adaptive LMS algorithm. Besides, the sliding discrete LCT computation by the adaptive system has been shown. The proposed computation methods are robust to the errors, and the inherent parallel structure makes the proposed computation techniques very suitable for efficient VLSI implementation.

Acknowledgments

This work was supported in part by the National Science Foundation of China for Distinguished Young Scholars under Grant No. 60625104 and by the National Natural Science Foundation of China under Grant No. 61201354, and by the Basic Science Foundation of Beijing Institute of Technology under Grant 20120542005.

References and links

1. H. M. Ozaktas, Z. Zalevsky, and M. A. Kutay, The Fractional Fourier Transform with Applications in Optics and Signal Processing (Wiley, 2001).

2. J. J. Healy and J. T. Sheridan, “Cases where the linear canonical transform of a signal has compact support or is band-limited,” Opt. Lett. 33(3), 228–230 (2008). [CrossRef] [PubMed]

3. A. Stern, “Why is the linear canonical transform so little known?” Proc. AIP860,225–234 (2006). [CrossRef]

4. D. F. V. James and G. S. Agarwal, “The generalized Fresnel transform and its applications to optics,” Opt. Commun. 126(4-6), 207–212 (1996). [CrossRef]

5. A. Koc, H. M. Ozaktas, C. Candan, and M. A. Kutay, “Digital computation of linear canonical transforms,” IEEE Trans. Signal Process. 56(6), 2383–2394 (2008). [CrossRef]

6. S. C. Pei and J. J. Ding, “Eigenfunctions of linear canonical transform,” IEEE Trans. Signal Process. 50(1), 11–26 (2002). [CrossRef]

7. B. Deng, R. Tao, and Y. Wang, “Convolution theorems for the linear canonical transform and their applications,” Science in China Ser. F 49(5), 592–603 (2006). [CrossRef]

8. L. Onural, A. Gotchev, H. M. Ozaktas, and E. Stoykova, “A survey of signal processing problems and tools in holographic three-dimensional television,” IEEE Trans. Circ. Syst. Video Tech. 17(11), 1631–1646 (2007). [CrossRef]

9. A. V. Oppenheim and R. W. Schafer, Digital Signal Processing (Prentice-Hall, 1975)

10. L. R. Rabiner and B. Gold, Theory and Application of Digital Signal Processing (Prentice-Hall, 1975).

11. S. C. Pei and J. J. Ding, “Closed-form discrete fractional and affine Fourier transforms,” IEEE Trans. Signal Process. 48(5), 1338–1353 (2000). [CrossRef]

12. F. S. Oktem and H. M. Ozaktas, “Exact relation between continuous and discrete linear canonical transforms,” IEEE Signal Process. Lett. 16(8), 727–730 (2009). [CrossRef]

13. J. Zhao, R. Tao, and Y. Wang, “Sampling rate conversion for linear canonical transform,” Signal Process. 88(11), 2825–2832 (2008). [CrossRef]

14. J. Zhao, R. Tao, Y. L. Li, and Y. Wang, “Uncertainty principles for linear canonical transform,” IEEE Trans. Signal Process. 57(7), 2856–2858 (2009). [CrossRef]

15. B. Widrow, J. R. Glover Jr, J. M. McCool, J. Kaunitz, C. S. Williams, R. H. Hearn, J. R. Zeidler, E. Dong Jr, and R. C. Goodlin, “Adaptive noise canceling: Principles and applications,” Proc. IEEE 63(12), 1692–1716 (1975). [CrossRef]

16. B. Widrow, J. M. McCool, M. G. Larimore, and C. R. Johnson Jr., “Stationary and nonstationary learning characteristics of the LMS adaptive filter,” Proc. IEEE 64(8), 1151–1162 (1976). [CrossRef]

17. B. Widrow and S. D. Stearns, Adaptive Signal Processing (Prentice-Hall, 1985).

18. B. Widrow, J. M. McCool, and M. Ball, “The complex LMS algorithm,” Proc. IEEE 63(4), 719–720 (1975). [CrossRef]

19. B. Widrow, P. Baudrenghien, M. Vetterli, and P. F. Titchener, “Fundamental relations between the LMS algorithm and the DFT,” IEEE Trans. Circ. Syst. 34(7), 814–820 (1987). [CrossRef]

Discrete linear canonical transform computation by adaptive method

Abstract

1. Introduction

2. Linear canonical transform and LMS adaptive filter

2.1 Linear canonical transform

2.2 LMS adaptive filter

3. Block discrete LCT computation by adaptive filter

4. Sliding discrete LCT computation by adaptive filter

5. Simulation results

6. Conclusion

Acknowledgments

References and links

Cited By

Figures (6)

Equations (26)

Optics Express