# DFRFT: A Classified Review of Recent Methods with Its Application

Authors: Ashutosh Kumar Singh; Rajiv Saxena
Journal: Journal of Engineering
 (2013)

---

## Abstract

In the literature, there are various algorithms available for computing the discrete fractional Fourier transform (DFRFT). In this paper, all the existing methods are reviewed, classified into four categories, and subsequently compared to find out the best alternative from the view point of minimal computational error, computational complexity, transform features, and additional features like security. Subsequently, the correlation theorem of FRFT has been utilized to remove significantly the Doppler shift caused due to motion of receiver in the DSB-SC AM signal. Finally, the role of DFRFT has been investigated in the area of steganography.


---

## Body

## 1. Introduction

Due to the inadequate performance of Fourier transform for the analysis of nonstationary signal, normally encountered in communication systems when either transmitter or receiver is moving, fractional Fourier transform (FRFT) can be utilized. The concept of FRFT was first introduced by N. Wiener in 1929. However, the FRFT was recognized as a “transform method” by mathematical bodies after the work of Victor Namias in 1980 in which the concept of FRFT had been introduced by considering fractional power of eigenfunctions of the ordinary FT. The mathematical description of FRFT was given by McBride and Keer in 1987. And finally, a general definition of FRFT for all classes of signals was given by Cariolaro et al. To compute the FRFT of any signal its discrete version was needed which initiates the work for defining discrete FRFT (DFRFT). However, due to availability of different approaches for evaluating DFRFT, a detailed study is needed in the context of suitable definition of DFRFT. This study along with comparison between the existing methods has been included in Section 3.This paper includes five sections. The basic concept of fractional Fourier transform is described in Section2. The detailed study of discrete FRFT (DFRFT) along with the comparative analysis between the different methods of evaluating DFRFT is included in Section 3. Subsequently, the utility of FRFT in DSB-SC signal detection under noisy condition and the application of DFRFT in steganography are being dealt in Section 4.  Finally, conclusions are made in Section 5.


## 2. Fractional Fourier Transform (FRFT)

The FRFT is a mathematical tool which maps a signal from one domain to other in time-frequency plane (e.g., FRFT becomes FT when it maps a signal in time domain to frequency domain), shown in Figure1. The FRFT of a signal x(t) with angle parameter α, represented by Xα(u) is defined as [1]:
(1)Fα[x(t)]=Xα(u)=1-jcotα2π∫-∞∞x(t)e(j/2){(t2+u2)cotα-2tucscα}dt
The angle parameter associated with the FRFT is given as:
(2)α=aπ2.
When α=0, the FRFT becomes identity operator, and for α=π/2(a=1), the FRFT becomes Fourier transform (FT). Similarly, for α=π(a=2) FRFT behaves as an inverse operator, and for α=3π/2(a=3), the expression converts into an inverse FT, as shown in Figure 1.Figure 1
FRFT in time-frequency plane.In signal processing, cross-correlation is a measure of similarity of two waveforms as a function of a time lag applied to one of them. The cross-correlation theorem for the FRFT is given as follows [2]. Let rxy(τ) be defined as weighted cross-correlation of two functions x(t) and y(t), and Xα(u), Yα(u), and Rxy,α(u), be defined as FRFT of x(t), y(t), and rxy(τ), respectively. Then,
(3)rxy(τ)=∫-∞∞x(t)y(t+τ)ejt(τ+t)cotαdt⟷Rxy,α(u)=2π1-jcotαe-(j/2)u2cotαXα(-u)Yα(u).
This cross-correlation theorem of FRFT is utilized for detection of DSB-SC signal when the signal encounters a Doppler shift due to motion of transmitter or receiver.


## 3. Discrete Fractional Fourier Transform (DFRFT)

Many researchers have tried to come up with some method of evaluating DFRFT. On this journey, many methods and algorithms had come into existence. But none got acclamation due to deficiency of either not satisfying some of the prime properties that its continuous type possessses or nonexistence of a closed-form expression. Previously, many classifications have already been documented [3]. However, in this section, all the approaches available in the literature are included and divided into four classes on the basis of their methodology of evaluation. These classes are(i)
sampling-based method(ii)
linear combination method(iii)
eigenvector-based method(iv)
weighted summation-based method.


### 3.1. Sampling-Based Method

In this approach, by considering the samples of the kernel matrix under some predefined constraint, two different algorithms are available in the literature and they are named as sampling type DFRFT and closed-form type DFRFT. These two methods are described below:


#### 3.1.1. Sampling Type DFRFT

The sampling theorem for the FRFT of band-limited and time-limited signals is followed from Shannon’s sampling theorem [4–9]. With this approach, the simplest way to derive the DFRFT is sampling the continuous FRFT and computing it directly from the samples, but by sampling the continuous FRFT directly, the resultant discrete transform obtained will lose many important properties. The most serious problem with this type DFRFT will be noncompliance of unitary and reversibility properties. In addition, the DFRFT obtained by direct sampling of the FRFT lacks closed-form expressions and is nonadditive, which affects its application domain. In order to maintain some of the FRFT properties, a type of DFRFT was derived as a special case of the continuous FRFT as reported in [10]. Specifically, the input function was assumed as a periodic, equally spaced impulse train. Since this type of DFRFT is a special case of continuous FRFT, many properties of the FRFT exist and have the fast algorithm. However, this type of DFRFT cannot be defined for all values of α due to various imposed constraints. At the same time, Ozaktas et al. [11] had proposed two innovative approaches for obtaining the DFRFT through the sampling of the FRFT. Both methods were based on the idea of manipulating the expression of the FRFT by sampling appropriately.In the first method, the DFRFT of a signalx(t) can be obtained by multiplying the signal with a chirp function as followes:
(4)x′(t)=e-jπt2tan(α/2)x(t).For this chirp multiplication, the bandwidth and time-bandwidth product ofx′(t) can be as large as twice of x(t). Thus, the sampling of x′(t) will be done at intervals of (1/2Δx). If the samples of x(t) are spaced at (1/Δx), then this needs to be interpolated before multiplying it by the samples of the chirp function to obtain the desired samples of x′(t).The next step was to convolve the signalx′(t) with another chirp function. To perform this convolution, the chirp signal was replaced by its band-limited version because x′(t) is assumed as a band-limited signal. 										(5)x′′(t)=Aα∫-∞∞x(τ)ejπ(t-τ)2cscαdτ.The band-limited version of chirp function can be expressed in the form of Fresnel integral. Finally, again performing a chirp multiplication, the DFRFT of signalx(t) can be obtained as:
(6)ℑ[x(t)]=Xα=e-jπt2tan(α/2)x′′(t).The convolution operation in the above equation can be achieved by samplingx(t) and performing the convolution using the fast Fourier transform (FFT). Overall, the procedure starts with N samples spaced at (1/Δx), which uniquely characterizes the function x(t) and returns the same for Xα. For x and Xα, denote column vectors with N elements containing the samples of x(t) and its DFRFT, respectively. The above procedure can be given in a matrix notation as:
(7)Xα=DΛHΛJx.Here,D and J are matrices representing the decimation and interpolation operations, Λ is a diagonal matrix that corresponds to chirp multiplication, and H corresponds to the convolution operation. The above method had allowed obtaining the samples of the ath transform in terms of the samples of the original function.First method includes the analysis of one step in terms of Fresnel integral, which requires large number of computation time particularly in inverse DFRFT. To remove this discrepancy, a second method was introduced by Ozaktas et al. [11]. In this method, one assumption was taken that the Wigner distribution of function x(t) could be zero outside a circle of diameter Δu centered at the origin and by limiting the order “a” to the interval 0.5≤|a|≤1.5. Simultaneously, the amount of vertical shear in Wigner space resulting from the chirp modulation is bounded by Δu/2. The FRFT can be written as:
(8)ℑ[x(t)]=Xα(u)=Aαejπu2cotα∫-∞∞e-jπtucscα{ejπt2cotαx(t)}dt.
Thus, ejπt2cotαx(t) can be written by using Shannon’s interpolation formula as:
(9)ejπt2cotαx(t)=∑n=-NNejπ(n/2Δu)2cotαx(n2Δu)sinc(2Δu(t-n2Δu)),
where N=(Δu)2; putting (9) in (8) and interchanging the order of integration and summation,
(10)Xα(u)=Aαejπu2cotα∑n=-NNejπ(n/2Δu)2cotαx(n2Δu)×∫-∞∞e-jπtucscαsinc(2Δu(t-n2Δu))dt.
And after some algebraic manipulation, (10) can be restructured as:
(11)Xα(m2Δu)=Aα2Δuejπ(cotα-cscα)(m/2Δu)2×∑n=-NNejπ(m-n/2Δu)2cscαejπ(cotα-cscα)(n/2Δu)2×x(n2Δu).This summation can be interpreted as the convolution ofejπ(n/2Δu)2cscα and the chirp modulated function ejπ(cotα-cscα)(n/2Δu)2x(n/2Δu). This convolution can be computed by fast Fourier transform (FFT). And the overall methodology can be written in matrix form as:
(12)Xα=DKαJx,
where the Kα matrix was given as:
(13)Kα(m,n)=Aα2Δuejπ{cotα(m/2Δu)2-2cscα(mn/(2Δu)2)+cotα(n/2Δu)2}.In the same article, Ozaktas et al. [11] gave alternate methods to calculate the DFRFT by introducing scaling process. All these methods resulted in the same time complexity. However, in practice, these alternate methods are not preferable because they require coordinate scaling. Scaling is not an advantageous process if the original data has already been sampled. This is because scaling requires additional interpolations and so forth, which will require additional computation.However, in both cases, the Wigner distribution of the signal was confined to a circle of diameterΔu around the origin. Thus, there might be several discrete fractional Fourier transform matrices, some of which were more elegantly expressible than the others, which gave the same result within the accuracy of this approximation. All these matrices would yield results that were in increasingly better agreement as the signal energy contained in the circle was increasingly closer to the total energy.


#### 3.1.2. Closed-Form Type DFRFT

Pei and Ding [12] had derived a different type of DFRFT which was neat in concept and reversible but lacks the additivity property. Due to the orientation of practical usage, two types of DFRFT were derived which are different in parameterization. The parameters of the first type have directly linked to the continuous FRFT and suit the applications of computing the continuous FRFT. On the other hand, the second type has the simpler parameters set and allows more elegant expression for the operator kernels and more is suitable for other applications of DFRFT, such as the filter design, pattern recognition, and phase retrieval. Considering the expression of continuous FRFT:
(14)Fα(u)=1-jcotα2πej(u2/2)cotα∫-∞∞e-jtucscαej(t2/2)cotαf(t)dt.
And to derive the DFRFT, sampling of f(t) and Fα(u) with the interval Δt and Δu needs to be performed as:
(15)y(n)=f(n*Δt),Yα(m)=Fα(m*Δu),
where n=-N,…,N and m=-M,…,M. By using this sampling, FRFT can be converted [13] as:
(16)Yα(m)=1-jcotα2πΔtej(m2Δu2/2)cotα×∑n=-NNe-jnmΔtΔucscαej(n2Δt2/2)cotαy(n).The above can be written as:(17)Yα(m)=∑n=-NNFα(m,n)y(n),
where
(18)Fα(m,n)=1-jcotα2πΔtej(m2Δu2/2)cotα×e-jnmΔtΔucscαej(n2Δt2/2)cotα.
And Yα(m) can be reversible only if the inverse transform of Fα(m,n), for all M≥N, would be Hermitian (conjugate and transpose), given [12] as:
(19)y(n)=∑m=-MMFα*(m,n)Yα(m)forM≥N.From (17) and (19),
(20)y(n)=∑m=-MM∑k=-NNFα*(m,n)Fα(m,k)y(k).After some manipulation,(21)y(n)=Δt22π|sinα|∑m=-MM∑k=-NN×ej((k2-n2)Δt2/2)cotαejm(n-k)ΔtΔucscαy(k).The summation form in (21) should be equal to δ(n-k), then
(22)ΔtΔu=S2πsinα2m+1,
where |S| is some integer prime to (2m+1); in this case, (18) becomes
(23)Fα(m,n)=1-jcotα2πΔtej(m2Δu2/2)cotαe-j(S2πnm/(2m+1))×ej(n2Δt2/2)cotα,∑m=-MM∑k=-NNFα*(m,n)Fα(m,k)y(k)=2M+12π|sinα|Δt2y(n).
After normalization, the transform matrix [12] can be given as:
(24)Fα(m,n)=|sinα|-jsgn(sinα)cosα2M+1ej(m2Δu2/2)cotα×e-j(sgn(sinα)2πnm/(2M+1))ej(n2Δt2/2)cotα.
The definition of DFRFT is further divided into two parts [12]. The first part was defined for “sinα>0” as:
(25)Yα(m)=sinα-jcosα2M+1ej(m2Δu2/2)cotα×∑n=-NNe-j(2πnm/(2M+1))ej(n2Δt2/2)cotαy(n).
And the second part was defined for “sinα<0” as:
(26)Yα(m)=-sinα+jcosα2M+1ej(m2Δu2/2)cotα×∑n=-NNej(2πnm/(2M+1))ej(n2Δt2/2)cotαy(n).
This definition of DFRFT possesses reversibility property and the periodicity property:
(27)Yα(-m)=Yα+π(m),Yα(m)=Yα+2π(m).The DFRFT of type 1 has a very important advantage in terms of its efficiency in calculating and implementing the DFRFT [12]. Due to two chirp multiplications and one FFT required for the implementation of this type of DFRFT, the total number of the multiplication operations required was “2P+(P/2)log2P,” for length of the output as P=2M+1.Subsequently, another type of DFRFT, named as type 2, was introduced by modifying the definitions obtained for discrete affine Fourier transform (DAFT) and presented as:(28)Yp(M)=12M+1ej(m2p/2)×∑n=-NNe±j(2πnm/(2M+1))ej(n2p/2)y(n),∀M≥N.For the DFRFT of type 2, the parameter “p” is used to control the variation of the chirp in frequency domain [13]. The DFRFT of type 2 also needed “2P+(P/2)log2P” multiplication operations. It does not follow the additivity property but observes the convertibility operation.These definitions of DFRFT are better to previous ones, that is, sampling type, linear combination type, and eigen-function type, as it require less number of multiplications, having less complexity. Also it has a closed form expression being rich in properties.


### 3.2. Linear Combination Type DFRFT

In [14–17], the discrete fractional Fourier transform was derived by using the linear combination of identity operation (F0), discrete Fourier transform (Fπ/2), time inverse operation (Fπ), and inverse discrete Fourier transform (F3π/2). In [15], the concept of DFRFT was given by linear combination of Lagrange interpolation polynomial of degree 3. For any function x(v)(29)x(v)=∑k=03ϕk(v)g(λk),
where
(30)λk=ejkπ/2,for0≤k≤3.Andϕk(v) is the Lagrange interpolation polynomial of degree 3 taking value “one” at v=λk and “zero” otherwise. It follows that
(31)ϕk(v)=14∑i=03uie-jkiπ/2.
Finally, by taking the principle nth-root of eigenvalue λk, the fractional transform can be given as:
(32)Fn=∑i=03Fiαi(n),
where
(33)αi(n)=14∑k=14e-jk(n-i)π/2,for0≤i≤3.Subsequently, in [14, 15], it had also been shown that the operator defined by (17) is unitary; satisfies angle additivity and angle multiplicity properties, and the operator is periodic with a fundamental period of four. One analogous approach of this method was presented in [17], in which the DFRFT was calculated by first chirping the signal in time domain, then taking its Fourier transform, and finally chirping again in frequency domain. Depending on the character of the chirped signal, the embedded FT can be one of the two known classes in discrete case—discrete-time and periodic discrete-time.For discrete-time DFRFT, the expression for DFRFT of signalx(t) was given as:
(34)Xα(f)=∑k=-∞∞Tx(kT)Kα(f,kT).
And its inverse transform as:
(35)x(kT)=∫0FpXα(f)Kα*(f,kT)df,
where Fp=2/Tcscα and Kα(t,u) is the kernel of FRFT. Similarly, for periodic discrete-time signal x(t) with period Tp, the expression for DFRFT was given as:
(36)Xα(nf)=∑k=-∞∞Tx(kT)Kα(nF,kT).
And its inverse transform as:
(37)x(kT)=∑n=0N-1FXα(nF)Kα*(nF,kT),
where Tp=NT and Fp=NF=2/Tcscα.In this type of DFRFT, the transform matrix is orthogonal and satisfies the additivity property along with the reversibility property. However, the main problem with this method is that the transform result is not matched to the continuous FRFT. Besides, it performed very similarly to the original Fourier transform or the identity operation but lost the important characteristic of fractionalization.


### 3.3. Eigenvector Type DFRFT

In this classification, the basic approach of eigenvalues and eigenvectors of the kernel matrix of the FRFT is considered to define the discrete FRFT. Here, the eigenvalues of the discrete Fourier transform (DFT) are considered first and then by obtaining the powers of these eigenvectors of DFT matrix, the corresponding eigenvalues and eigenvectors of DFRFT matrix are derived. This classification includes three different approaches, named as, eigenvector decomposition method, random type method, andV¯ method DFRFT.


#### 3.3.1. Eigenvector Decomposition Type DFRFT

The authors had derived another type of discrete fractional Fourier transform as reported in [13, 18–25] by searching the eigenvectors and eigenvalues of the DFT matrix followed by computing the fractional power of the DFT matrix based on Hermite function. This type of DFRFT worked very similarly to the continuous FRFT, and it fulfills the properties of orthogonality, additivity, and reversibility. The fractional power of matrix could be calculated from its eigendecomposition and the power of eigenvalues. Unfortunately, there exist two types of ambiguity in deciding the fractional power of the DFT kernel matrix.(i)
Ambiguity in Deciding the Fractional Powers of Eigenvalues.The square roots of unity are “1” and “–1” from elementary mathematics. This indicates that there exists root ambiguity in deciding the fractional power of eigenvalues.(ii)
Ambiguity in Deciding the Eigenvectors of the DFT Kernel Matrix.The DFT eigenvectors constitute four major eigensubspaces; therefore, the choices for the DFT eigenvectors to construct the DFRFT kernel are multiple and not unique.Because of the above mentioned ambiguities, several DFRFT kernel matrices were possible which could obey the rotational properties. The idea for developing the DFRFT is to find the discrete form of the following:(38)Kα(t,u)=∑n=0∞e-jnαHn(t)Hn(u),
where α indicates the rotation angle associated with FRFT of signal and Hn(t) represents the nth-order normalized Hermite function with unit variance. The nth-order normalized Hermite function with variance σ is defined as:
(39)Hn(t)=12nn!σπhn(tσ)e-(t2/2σ2),
where, hn(·) represents the nth-order Hermite polynomial. However, the normalized Hermite functions with unitary variance Hn(·) indicate the eigenfunction of the FRFT. With the help of Hermite functions, the FRFT of a signal can be computed [24] as:
(40)Xα(u)=∫-∞∞x(t)Kα(t,u)dt=∑n=0∞Hn(u)e-jnα{∫-∞∞x(t)Hn(t)dt}.The FRFT can be interpreted as a weighting summation of Hermite functions, as given in (40). The weighting coefficients are obtained by multiplying the phase term e-jnα and the inner product of the input signal with the corresponding Hermite functions [24]. The development of this DFRFT was based on the eigendecomposition of the DFT kernel. The eigenvalues of DFT kernel “F” can be given as {1,-j,-1,j}, and corresponding to each eigenvalue one eigenvector can be evaluated. The continuous FRFT has a Hermite function with unitary variance as its eigenfunction. The corresponding eigenfunction property for the DFT can be given as:
(41)F2α/π[u^n]=e-jnαu^n,
where u^n represents the eigenvector of DFT corresponding to the nth-order discrete Hermite functions. In order to retain the eigenfunction properties, the unit variance Hermite function is sampled with a period Ts=2π/N. In the case of continuous FRFT, the terms of the Hermite functions are summed up from order zero to infinity. However, for the discrete case, only eigenvectors for the DFT Hermite eigenvectors can be added. The selection of the DFT Hermite eigenvectors is usually made from low to high orders, due to small approximation error of the low DFT Hermite eigenvectors [24]. Finally, the transform kernel of DFRFT can be defined as:
(42)F2α/π=U^D2α/πU^T={∑k=0N-1e-jkαu^ku^kT,∀N=4m+1{∑k=0N-2e-jkαu^ku^kT}+e-jNαu^Nu^NT,u^5NT∀N=4m,
where U^=[u^0|u^1|⋯u^N-1|] for N as odd and U^=[u^0|u^1|⋯u^N-2|u^N] for N is even, u^k represents the normalized eigenvector corresponding to the kth-order discrete Hermite function, and D can be defined as follows.ForN odd,
(43)D2α/π=[e-j00e-jα⋱0e-jα(N-1)].
And for N even,
(44)D2α/π=[e-j00⋱e-jα(N-2)0e-jαN]The convergence for the eigenvectors obtained from these matrixes is not so fast for the high-order Hermite functions. To ensure orthogonality of the DFT Hermite eigenvectorsu^k, the Gram-Schmidt algorithm (GSA), or the orthogonal procrustes algorithm (OPA) can be used [24]. The GSA minimizes the errors between the samples of the Hermite functions and the orthogonal DFT Hermite eigenvectors. On the other hand, the OPA minimizes the total errors between those samples. The main difference between the approach proposed in [24] and similar approaches proposed in [21–23] is found in the eigenvectors obtained. The eigenvectors obtained by method reported in [24] were discrete Mathieu function, although the Mathieu functions can converge to Hermite functions as obtained by the methods reported in [21–23]. Given the number of points N and the rotation angle α, the DFT Hermite eigenvectors can be computed followed by determining the eigenvalues of DFRFT. The computation of the DFRFT can be implemented only by a transform kernel matrix multiplication. The complexity of computing the DFRFT is O(N2) same as in the DFT case. By adjusting the rotational angles, the method for implementing the DFRFT of a signal x can be given as:
(45)F2α/π=U^D2α/πU^Tx=∑k=0N-1ake-jkαu^k.The coefficientsak’s are the inner products of signal and eigenvectors, and they can be computed in advance. If rotation angle changed then only the diagonal matrix D2α/π needs to be recomputed. Additionally, the authors [24] investigated the relationship between the FRFT and the DFRFT and found that for a sampling period equal to 2π/N, the DFRFT performs a circular rotation of the signal in the time-frequency plane. However, the DFRFT becomes an elliptical rotation in the continuous time-frequency plane for sampling periods different from 2π/N [24]. Therefore, for these elliptical rotations an angle modification and a postphase compensation in the DFRFT have been required to obtain results similar to the continuous FRFT [24]. This approach has been extended to the so-called multiple-parameter discrete fractional Fourier transform (MPDFRFT) [20, 25]. In fact, the MPDFRFT maintains all desired properties and reduces to the DFRFT when all of its order parameters are the same.A similar approach to [24] has been proposed in [19]. However, authors in [19] believe that the discrete time counterparts of the continuous time Hermite-Gaussians maintained the same properties because these discrete time counterparts exhibited better approximations than the other proposed approaches [18]. In order to resolve this issue about the approximation of Hermite-Gaussian functions, a nearly tridiagonal commuting matrix of the DFT and a corresponding version of the DFRFT were proposed in [13]. Most of the eigenvectors of this proposed nearly tridiagonal matrix result in a good approximation of the continuous Hermite-Gaussian functions by providing a smaller approximation error in comparison to the previous approaches.


#### 3.3.2. Random Type DFRFT

Pei and Hsue [26] had introduced a modified method to calculate the DFRFT based on the random DFT eigenvectors and eigenvalues. In this method, a new commuting matrix with random DFT eigenvectors was constructed. Then, a random discrete fractional Fourier transform (RDFRFT) kernel matrix with random DFT eigenvectors and eigenvalues was developed. The magnitude and phase of the transformed output are found to be random in RDFRFT. This can be considered as a special feature of RDFRFT. Previously, the DFT was randomized to define the discrete random Fourier transform (DRFT) in [27] by taking random powers of eigenvalues of the DFT matrix. However, the eigenvectors of the DRFT were not random and were computed using the same method as proposed in [23]. Although, the phase of the DRFT was random, its magnitude was not random, as shown in [27]. The eigen-decomposition of DFT matrix “F” was given as:
(46)F=∑k=0N-1λkekekT.
The DFRFT with one order parameter “a” was defined by [25]:
(47)Fa=∑k=0N-1λkaekekT.The definition of RDFRFT was based on the definition of multiple-parameter DFRFT (MPDFRFT) [25], given as:
(48)Fa^=∑k=0N-1λkakekekT,
where the 1×N parameter vector a^ was defined as:
(49)a^=[a0,a1,…,aN-1].This order parameter of the MPDFRFT can be taken as independent random numbers to define an MPDFRFT with the following random eigenvalues [26], λkak, k=0,1,…,N-1, and random eigenvalues were given as:
(50)λkak={(e-j2π)ak,forλk=1(e-jπ/2)ak,forλk=-j(e-jπ)ak,forλk=-1(e-j3π/2)ak,forλk=j.The generation of RDFRFT is given in the following section. In this generation of RDFRFT, the DFT-commuting matrix with orthonormal random DFT eigenvectors is used.First defining aN×N real random Matrix “A” being “K” symmetric and generated by A=KAK, where K represents circular reversal matrix, which is given by:
(51)K=F2=F-2=[100JN-1],
where, JN-1 is  the (N-1)×(N-1) reversal matrix, having nonzero entries “ones” placed on the antidiagonal locations. Taking the Ksymmetric part “E” of the random matrix “A” as B={(A+KAK)/2}, with the help of this matrix, a random matrix “C” was generated as C=(B+BT)/2, where “T” denotes the matrix transpose operation, and finally a DFT commuting matrix “D” will be formed [27] as D=C+FCF-1. Subsequently, with the random DFT-commuting matrix “D” and the MPDFRFT, the random DFRFT (RDFRFT) with the 1×N parameter vector a^=[a0,a1,…,aN-1] of a signal “x” can be determined [26] as:
(52)X(a^,D)=FDa^*x={∑k=0N-1λkakrkrkT}*x,
where rk represents the orthonormal random DFT eigenvectors computed from “D” and λk represents the DFT eigenvalue corresponding to rk, given by (50).The RDFRFT kernel matrixFDa^ satisfied various properties of the FRFT like—additivity, FT convertibility, and angle additivity. The computational complexity of the -point RDFRFT was O(N2) same as that of the MPDFRFT, but RDFRFT has an advantage over MPDFRFT in terms of one more free parameter associated with it as parameter vector a^. The magnitude and phase responses of the transformed output, with the RDFRFT, are random because the eigenvectors and eigenvalues of the RDFRFT matrix are random [26], as shown in (52). This aspect establishes its advantageous role in security areas like—cryptography, encryption, and watermarking.


#### 3.3.3. TheV¯ Method DFRFT

In this work [28], authors focus on a different approach to derive eigenvectors of the DFT matrix in order to build the DFRFT matrix that mimics the properties of its continuous counterpart. A straightforward and refined derivation of the eigenvectors of the DFT matrix is proposed without using any commuting matrices, but using only the factored form of a basic property of the DFT matrix stating that four consecutive Fourier transforms is the identity transform
(53)WN4=IN,
where WN4 and IN are the centered DFT (CDFT) and the identity matrices of order N, respectively. We utilize factorization of the CDFT matrix to obtain new matrices, whose columns are eigenvectors of the CDFT matrix. Thereafter, GSA was used to find the orthonormal eigenvectors obtained by factorization of the CDFT matrix. By employing a DFT-shift matrix K, eigenvectors of the ordinary-DFT matrix are obtained by using the CDFT matrix. Later on, the S and T matrices together are also used to boost the performance.In [24], Pei et al. compute the samples of the Hermite-Gauss functions and employ Gram-Schmidt algorithm (GSA) to orthogonalize these samples by projecting on the S matrix, which is equivalent to orthogonalizing the samples of the Hermite-Gauss functions without using the S matrix. Hanna et al. [29] claim to find the eigenvectors of the DFT matrix using a similar method. However, they use the ordinary-DFT matrix as a basis to their work, but Serbes and Durak-Ata [28] employs the CDFT matrix. Additionally, the eigenvectors are nonorthogonal and dissimilar to the samples of Hermite-Gauss functions. Serbes and Durak-Ata [28] work is different from these methods in the sense that the non-orthogonal eigenvectors of the CDFT matrix is determined without utilizing any S matrix or samples of the Hermite-Gauss functions, and then orthogonalization of eigenvectors is obtained by employing the GSA.The DFT maps the signal from[0,N-1] discrete input space to [0,2π] discrete frequency space, whereas the CDFT maps [-(N-1)/2,(N-1)/2] to [-π,π]. This allows CDFT to define even and odd functions such as Hermite-Gauss functions, that is contrary to the ordinary DFT definition. In order to approximate the samples of the continuous FrFT and to imply the rotation property, Hermite-Gauss-like eigenvectors of the DFT matrix have to be obtained. Serbes and Durak-Ata [28] use the CDFT matrix instead of the ordinary DFT matrix, because the Hermite-Gauss-like vectors are eigenvectors of only the CDFT matrix, not the DFT matrix, that is, when a discrete Hermite-Gauss-like eigenvector is transformed using the ordinary DFT matrix, the output is the shifted version of the Hermite-Gauss vector multiplied by a complex sinusoidal.The FRFT operator of orderα can be defined as the αth-power of the ordinary DFT operator WN. Hence, the DFRFT matrix can be expressed by means of its eigenvector decomposition,
(54)WNα=UNΛNαUNT,
where ΛNα is given as—ΛNα=diag(e-j0,e-jπ/2a,…,e-jπ/2a(N-2),e-jπ/2a(N-1)) for the centered-DFRFT.The centered-DFRFT matrix can be obtained by calculating the Hermite-Gaussian-like eigenvectors of DFRFT as followes:(55)WNα=V1¯Λ¯k1αV¯1T+V2¯Λ¯k2αV¯2T+V3¯Λ¯k3αV¯3T+V4¯Λ¯k4αV¯4T,
where
(56)Λ¯k1α=diag(e-j2π(k1-1)a,…,e-j2πa,e-j0),Λ¯k2α=diag(e-j(π+(k2-1)2π)a,…,e-j3πa,e-jπa),Λ¯k3α=diag(e-j(3π/2+(k3-1)2π)a,…,e-j7πa/2,e-j3πa/2),Λ¯k4α=diag(e-j(π/2+(k4-1)2π)a,…,e-j5πa/2,e-jπa/2).The above mentioned method satisfies unitary and angle additivity properties. Also, it reduces to the ordinary DFT whena=1 and approximates to the samples of the continuous FRFT for fractional values of a.Performance of the Serbes and Durak-Ata [28] algorithm can be improved by using it in combination with S and T matrices. Pei et al. [13] uses the S method together with the Grunbaum’s T matrix [30] as a linear combination, in which it has been observed that Hermite-Gauss-like eigenvectors of (S+15T) are more accurate than both S and T alone. Serbes and Durak-Ata [28] proposed that eigenvectors of linear combinations of S, T, and VT matrices as k1S+k1T+k1VT produce more accurate eigenvectors. Pei and Hsue [25] showed that using the linear combination of S and T as (S+15T) produces better results than the S or T alone. Whereas, Serbes and Durak-Ata [28] showed that using (S+30T-7VT) produces a more accurate approximation to the continuous FRFT. Although, better commuting matrix generation method is available in the literature [31–33], but all these methods are explained in the context of DFT and some other transform only, excluding DFRFT. Therefore, the betterment offered by these methods can be utilized in devising the new algorithm of DFRFT based on these commuting matrixes.


### 3.4. Weighted Summation Type DFRFT

Yeh and Pei [34] had developed a new method for DFRFT computation. The idea of the developed method was to compute the DFRFT at any angle by a weighted summation of the DFRFT’s with the special angles. In this method DFRFT computation with odd point of length can be realized by the weighted summation of the DFRFTs in special angles. The special angles were multiples of 2π/N for the odd case and multiples of 2π/(N+1) for even length. Regardless of even- or odd-length cases, the weighting coefficients were obtained from an IDFT operation. The N-point IDFT was needed for odd length, and the (N+1) - point was computed for even length.For discrete signal “x” having odd length “N”, the DFRFT of “x” for rotation angle “α” can be computed as:
(57)Xα=∑n=0N-1Bn,αXn,β∀β=2πN,
where the weighting coefficients Bn,α was given as:
(58)Bn,α=IDFT[e-jkα]∀k=0,1,…,N-1=1N∑n=0N-1e-jkαej(2π/N)nk.
The weighted summation Bn,α has a closed-form solution:
(59)Bn,α={1N1-e-j(N-1)(α-nβ)1-e-j(α-nβ),∀α≠kβδ(n-k),∀α=kβ.Similarly, for discrete signal “x” having even length N, the DFRFT of “x” for rotation angle “α” can be computed by the following expression:
(60)Xα=∑n=0NBn,αXn,β∀β=2πN+1.
The weighting coefficients, Bn,α can be computed as:
(61)Bn,α=IDFT[e-jkα]∀k=0,1,…,N=1N+1∑n=0Ne-jkαej(2π/(N+1))nk.
Similar to the odd length, the weighted summation Bn,α for even case has also a closed-form solution:
(62)Bn,α={1N+11-e-jN(α-nβ)1-e-j(α-nβ)∀α≠kβδ(n-k)∀α=kβThe eigenvector decomposition methods requireN2 storages to store the transform kernel. If the computing angle is changed, the transform kernel is also changed and needs to be recomputed. In this algorithm, the DFRFT of the specified angle was stored. The DFRFT at any angle can be obtained by weighted summation of these specified DFRFTs. The weighted computation in this algorithm still takes O(N2) multiplications; therefore, its computation load is still O(N2). In [34], two implementation methods for the DFRFT computation algorithm were introduced. One was called the parallel method and the other cascade method. Similar to the parallel method, the special angle and numbers of terms of cascade forms are also different for the even and odd cases. Both of the implementation methods can have the advantages over each other for different application. The parallel method was suitable for the signal, whose DFRFTs with special angles were already known. Hence, the computation of the DFRFT will become only a linear combination of the DFRFT’s in special angles. And, the cascade method means that the computation of the DFRFT can be realized from the DFRFT with only one specified angle. If the DFRFT with the specified angle can be computed efficiently, the computation of the DFRFT will become efficient. Therefore, it has been found suitable for VLSI implementation [34].Although, the existence of one well-accepted closed-form definition of DFRFT is not available; however many definitions with their pros and cons are accessible. A comparison of all the proposed classes of DFRFT based on some properties and complexity is shown in Table1. In this Table, yes is included against the method which possesses given property/parameter, otherwise NO is mentioned.Table 1
Comparison of Different Types of DFRFT Algorithms.
Properties/parameters
Sampling-based method
Linear combination-based method
Eigenvector-based method
Weighted summation-based method
Sampling type
Closed-form type
Linearcombination type
Eigenvector decomposition type
Random type
V
¯ Method type
Weighted summation Type
Reversibility
NO
YES
YES
YES
YES
YES
YES
Additivity
NO
NO
YES
YES
YES
YES
YES
Periodicity
YES
YES
YES
YES
YES
YES
YES
Symmetry
YES
YES
NO
YES
YES
YES
YES
Angle additivity
NO
YES
NO
YES
YES
YES
YES
Closed-form expression
YES
YES
YES
YES
YES
YES
YES
Similarity with FRFT
YES
YES
NO
YES
YES
YES
NO
Complexity
O(N log2N)
O(N log2N)
O(N log2N)
O(N2)
O
(
N
2
)
O(N log2N)
O(N2)
Constraints
Many
Fewer
Many
Many
Very less
Many
FewerFrom the above mentioned comparative analysis, it can be recapitulated that RDFRFT provides a better solution for computation of the DFRFT. The computation of RDFRFT requires more time than sampling type, linear combination type, and closed form type; however, it satisfies more properties than any other method with very less number of imposed constraints. Previously, IEEE has conferred to use the Eigenvectors decomposition type DFRFT as the most appropriate definition of DFRFT. In the year 2006, Pei et al. [25] had established the superiority of their method (MPDFRFT) over the earlier reported methods. Subsequently, in the year 2009, the MPDFRFT was modified by Pei et al. [26] and given as RDFRFT. It was also established that RDFRFT has one more degree of freedom, while satisfying the same set of properties as satisfied by MPDFRFT with equal computational complexity. Although, Serbes and Durak-Ata [28] proposed a modified version of eigenvalue decomposition technique with slight betterment but this does not redeem the advantage offered by RDFRFT, that is, security of the data. This makes the RDFRFT as a better choice in security applications.


## 3.1. Sampling-Based Method

In this approach, by considering the samples of the kernel matrix under some predefined constraint, two different algorithms are available in the literature and they are named as sampling type DFRFT and closed-form type DFRFT. These two methods are described below:


### 3.1.1. Sampling Type DFRFT

The sampling theorem for the FRFT of band-limited and time-limited signals is followed from Shannon’s sampling theorem [4–9]. With this approach, the simplest way to derive the DFRFT is sampling the continuous FRFT and computing it directly from the samples, but by sampling the continuous FRFT directly, the resultant discrete transform obtained will lose many important properties. The most serious problem with this type DFRFT will be noncompliance of unitary and reversibility properties. In addition, the DFRFT obtained by direct sampling of the FRFT lacks closed-form expressions and is nonadditive, which affects its application domain. In order to maintain some of the FRFT properties, a type of DFRFT was derived as a special case of the continuous FRFT as reported in [10]. Specifically, the input function was assumed as a periodic, equally spaced impulse train. Since this type of DFRFT is a special case of continuous FRFT, many properties of the FRFT exist and have the fast algorithm. However, this type of DFRFT cannot be defined for all values of α due to various imposed constraints. At the same time, Ozaktas et al. [11] had proposed two innovative approaches for obtaining the DFRFT through the sampling of the FRFT. Both methods were based on the idea of manipulating the expression of the FRFT by sampling appropriately.In the first method, the DFRFT of a signalx(t) can be obtained by multiplying the signal with a chirp function as followes:
(4)x′(t)=e-jπt2tan(α/2)x(t).For this chirp multiplication, the bandwidth and time-bandwidth product ofx′(t) can be as large as twice of x(t). Thus, the sampling of x′(t) will be done at intervals of (1/2Δx). If the samples of x(t) are spaced at (1/Δx), then this needs to be interpolated before multiplying it by the samples of the chirp function to obtain the desired samples of x′(t).The next step was to convolve the signalx′(t) with another chirp function. To perform this convolution, the chirp signal was replaced by its band-limited version because x′(t) is assumed as a band-limited signal. 										(5)x′′(t)=Aα∫-∞∞x(τ)ejπ(t-τ)2cscαdτ.The band-limited version of chirp function can be expressed in the form of Fresnel integral. Finally, again performing a chirp multiplication, the DFRFT of signalx(t) can be obtained as:
(6)ℑ[x(t)]=Xα=e-jπt2tan(α/2)x′′(t).The convolution operation in the above equation can be achieved by samplingx(t) and performing the convolution using the fast Fourier transform (FFT). Overall, the procedure starts with N samples spaced at (1/Δx), which uniquely characterizes the function x(t) and returns the same for Xα. For x and Xα, denote column vectors with N elements containing the samples of x(t) and its DFRFT, respectively. The above procedure can be given in a matrix notation as:
(7)Xα=DΛHΛJx.Here,D and J are matrices representing the decimation and interpolation operations, Λ is a diagonal matrix that corresponds to chirp multiplication, and H corresponds to the convolution operation. The above method had allowed obtaining the samples of the ath transform in terms of the samples of the original function.First method includes the analysis of one step in terms of Fresnel integral, which requires large number of computation time particularly in inverse DFRFT. To remove this discrepancy, a second method was introduced by Ozaktas et al. [11]. In this method, one assumption was taken that the Wigner distribution of function x(t) could be zero outside a circle of diameter Δu centered at the origin and by limiting the order “a” to the interval 0.5≤|a|≤1.5. Simultaneously, the amount of vertical shear in Wigner space resulting from the chirp modulation is bounded by Δu/2. The FRFT can be written as:
(8)ℑ[x(t)]=Xα(u)=Aαejπu2cotα∫-∞∞e-jπtucscα{ejπt2cotαx(t)}dt.
Thus, ejπt2cotαx(t) can be written by using Shannon’s interpolation formula as:
(9)ejπt2cotαx(t)=∑n=-NNejπ(n/2Δu)2cotαx(n2Δu)sinc(2Δu(t-n2Δu)),
where N=(Δu)2; putting (9) in (8) and interchanging the order of integration and summation,
(10)Xα(u)=Aαejπu2cotα∑n=-NNejπ(n/2Δu)2cotαx(n2Δu)×∫-∞∞e-jπtucscαsinc(2Δu(t-n2Δu))dt.
And after some algebraic manipulation, (10) can be restructured as:
(11)Xα(m2Δu)=Aα2Δuejπ(cotα-cscα)(m/2Δu)2×∑n=-NNejπ(m-n/2Δu)2cscαejπ(cotα-cscα)(n/2Δu)2×x(n2Δu).This summation can be interpreted as the convolution ofejπ(n/2Δu)2cscα and the chirp modulated function ejπ(cotα-cscα)(n/2Δu)2x(n/2Δu). This convolution can be computed by fast Fourier transform (FFT). And the overall methodology can be written in matrix form as:
(12)Xα=DKαJx,
where the Kα matrix was given as:
(13)Kα(m,n)=Aα2Δuejπ{cotα(m/2Δu)2-2cscα(mn/(2Δu)2)+cotα(n/2Δu)2}.In the same article, Ozaktas et al. [11] gave alternate methods to calculate the DFRFT by introducing scaling process. All these methods resulted in the same time complexity. However, in practice, these alternate methods are not preferable because they require coordinate scaling. Scaling is not an advantageous process if the original data has already been sampled. This is because scaling requires additional interpolations and so forth, which will require additional computation.However, in both cases, the Wigner distribution of the signal was confined to a circle of diameterΔu around the origin. Thus, there might be several discrete fractional Fourier transform matrices, some of which were more elegantly expressible than the others, which gave the same result within the accuracy of this approximation. All these matrices would yield results that were in increasingly better agreement as the signal energy contained in the circle was increasingly closer to the total energy.


### 3.1.2. Closed-Form Type DFRFT

Pei and Ding [12] had derived a different type of DFRFT which was neat in concept and reversible but lacks the additivity property. Due to the orientation of practical usage, two types of DFRFT were derived which are different in parameterization. The parameters of the first type have directly linked to the continuous FRFT and suit the applications of computing the continuous FRFT. On the other hand, the second type has the simpler parameters set and allows more elegant expression for the operator kernels and more is suitable for other applications of DFRFT, such as the filter design, pattern recognition, and phase retrieval. Considering the expression of continuous FRFT:
(14)Fα(u)=1-jcotα2πej(u2/2)cotα∫-∞∞e-jtucscαej(t2/2)cotαf(t)dt.
And to derive the DFRFT, sampling of f(t) and Fα(u) with the interval Δt and Δu needs to be performed as:
(15)y(n)=f(n*Δt),Yα(m)=Fα(m*Δu),
where n=-N,…,N and m=-M,…,M. By using this sampling, FRFT can be converted [13] as:
(16)Yα(m)=1-jcotα2πΔtej(m2Δu2/2)cotα×∑n=-NNe-jnmΔtΔucscαej(n2Δt2/2)cotαy(n).The above can be written as:(17)Yα(m)=∑n=-NNFα(m,n)y(n),
where
(18)Fα(m,n)=1-jcotα2πΔtej(m2Δu2/2)cotα×e-jnmΔtΔucscαej(n2Δt2/2)cotα.
And Yα(m) can be reversible only if the inverse transform of Fα(m,n), for all M≥N, would be Hermitian (conjugate and transpose), given [12] as:
(19)y(n)=∑m=-MMFα*(m,n)Yα(m)forM≥N.From (17) and (19),
(20)y(n)=∑m=-MM∑k=-NNFα*(m,n)Fα(m,k)y(k).After some manipulation,(21)y(n)=Δt22π|sinα|∑m=-MM∑k=-NN×ej((k2-n2)Δt2/2)cotαejm(n-k)ΔtΔucscαy(k).The summation form in (21) should be equal to δ(n-k), then
(22)ΔtΔu=S2πsinα2m+1,
where |S| is some integer prime to (2m+1); in this case, (18) becomes
(23)Fα(m,n)=1-jcotα2πΔtej(m2Δu2/2)cotαe-j(S2πnm/(2m+1))×ej(n2Δt2/2)cotα,∑m=-MM∑k=-NNFα*(m,n)Fα(m,k)y(k)=2M+12π|sinα|Δt2y(n).
After normalization, the transform matrix [12] can be given as:
(24)Fα(m,n)=|sinα|-jsgn(sinα)cosα2M+1ej(m2Δu2/2)cotα×e-j(sgn(sinα)2πnm/(2M+1))ej(n2Δt2/2)cotα.
The definition of DFRFT is further divided into two parts [12]. The first part was defined for “sinα>0” as:
(25)Yα(m)=sinα-jcosα2M+1ej(m2Δu2/2)cotα×∑n=-NNe-j(2πnm/(2M+1))ej(n2Δt2/2)cotαy(n).
And the second part was defined for “sinα<0” as:
(26)Yα(m)=-sinα+jcosα2M+1ej(m2Δu2/2)cotα×∑n=-NNej(2πnm/(2M+1))ej(n2Δt2/2)cotαy(n).
This definition of DFRFT possesses reversibility property and the periodicity property:
(27)Yα(-m)=Yα+π(m),Yα(m)=Yα+2π(m).The DFRFT of type 1 has a very important advantage in terms of its efficiency in calculating and implementing the DFRFT [12]. Due to two chirp multiplications and one FFT required for the implementation of this type of DFRFT, the total number of the multiplication operations required was “2P+(P/2)log2P,” for length of the output as P=2M+1.Subsequently, another type of DFRFT, named as type 2, was introduced by modifying the definitions obtained for discrete affine Fourier transform (DAFT) and presented as:(28)Yp(M)=12M+1ej(m2p/2)×∑n=-NNe±j(2πnm/(2M+1))ej(n2p/2)y(n),∀M≥N.For the DFRFT of type 2, the parameter “p” is used to control the variation of the chirp in frequency domain [13]. The DFRFT of type 2 also needed “2P+(P/2)log2P” multiplication operations. It does not follow the additivity property but observes the convertibility operation.These definitions of DFRFT are better to previous ones, that is, sampling type, linear combination type, and eigen-function type, as it require less number of multiplications, having less complexity. Also it has a closed form expression being rich in properties.


## 3.1.1. Sampling Type DFRFT

The sampling theorem for the FRFT of band-limited and time-limited signals is followed from Shannon’s sampling theorem [4–9]. With this approach, the simplest way to derive the DFRFT is sampling the continuous FRFT and computing it directly from the samples, but by sampling the continuous FRFT directly, the resultant discrete transform obtained will lose many important properties. The most serious problem with this type DFRFT will be noncompliance of unitary and reversibility properties. In addition, the DFRFT obtained by direct sampling of the FRFT lacks closed-form expressions and is nonadditive, which affects its application domain. In order to maintain some of the FRFT properties, a type of DFRFT was derived as a special case of the continuous FRFT as reported in [10]. Specifically, the input function was assumed as a periodic, equally spaced impulse train. Since this type of DFRFT is a special case of continuous FRFT, many properties of the FRFT exist and have the fast algorithm. However, this type of DFRFT cannot be defined for all values of α due to various imposed constraints. At the same time, Ozaktas et al. [11] had proposed two innovative approaches for obtaining the DFRFT through the sampling of the FRFT. Both methods were based on the idea of manipulating the expression of the FRFT by sampling appropriately.In the first method, the DFRFT of a signalx(t) can be obtained by multiplying the signal with a chirp function as followes:
(4)x′(t)=e-jπt2tan(α/2)x(t).For this chirp multiplication, the bandwidth and time-bandwidth product ofx′(t) can be as large as twice of x(t). Thus, the sampling of x′(t) will be done at intervals of (1/2Δx). If the samples of x(t) are spaced at (1/Δx), then this needs to be interpolated before multiplying it by the samples of the chirp function to obtain the desired samples of x′(t).The next step was to convolve the signalx′(t) with another chirp function. To perform this convolution, the chirp signal was replaced by its band-limited version because x′(t) is assumed as a band-limited signal. 										(5)x′′(t)=Aα∫-∞∞x(τ)ejπ(t-τ)2cscαdτ.The band-limited version of chirp function can be expressed in the form of Fresnel integral. Finally, again performing a chirp multiplication, the DFRFT of signalx(t) can be obtained as:
(6)ℑ[x(t)]=Xα=e-jπt2tan(α/2)x′′(t).The convolution operation in the above equation can be achieved by samplingx(t) and performing the convolution using the fast Fourier transform (FFT). Overall, the procedure starts with N samples spaced at (1/Δx), which uniquely characterizes the function x(t) and returns the same for Xα. For x and Xα, denote column vectors with N elements containing the samples of x(t) and its DFRFT, respectively. The above procedure can be given in a matrix notation as:
(7)Xα=DΛHΛJx.Here,D and J are matrices representing the decimation and interpolation operations, Λ is a diagonal matrix that corresponds to chirp multiplication, and H corresponds to the convolution operation. The above method had allowed obtaining the samples of the ath transform in terms of the samples of the original function.First method includes the analysis of one step in terms of Fresnel integral, which requires large number of computation time particularly in inverse DFRFT. To remove this discrepancy, a second method was introduced by Ozaktas et al. [11]. In this method, one assumption was taken that the Wigner distribution of function x(t) could be zero outside a circle of diameter Δu centered at the origin and by limiting the order “a” to the interval 0.5≤|a|≤1.5. Simultaneously, the amount of vertical shear in Wigner space resulting from the chirp modulation is bounded by Δu/2. The FRFT can be written as:
(8)ℑ[x(t)]=Xα(u)=Aαejπu2cotα∫-∞∞e-jπtucscα{ejπt2cotαx(t)}dt.
Thus, ejπt2cotαx(t) can be written by using Shannon’s interpolation formula as:
(9)ejπt2cotαx(t)=∑n=-NNejπ(n/2Δu)2cotαx(n2Δu)sinc(2Δu(t-n2Δu)),
where N=(Δu)2; putting (9) in (8) and interchanging the order of integration and summation,
(10)Xα(u)=Aαejπu2cotα∑n=-NNejπ(n/2Δu)2cotαx(n2Δu)×∫-∞∞e-jπtucscαsinc(2Δu(t-n2Δu))dt.
And after some algebraic manipulation, (10) can be restructured as:
(11)Xα(m2Δu)=Aα2Δuejπ(cotα-cscα)(m/2Δu)2×∑n=-NNejπ(m-n/2Δu)2cscαejπ(cotα-cscα)(n/2Δu)2×x(n2Δu).This summation can be interpreted as the convolution ofejπ(n/2Δu)2cscα and the chirp modulated function ejπ(cotα-cscα)(n/2Δu)2x(n/2Δu). This convolution can be computed by fast Fourier transform (FFT). And the overall methodology can be written in matrix form as:
(12)Xα=DKαJx,
where the Kα matrix was given as:
(13)Kα(m,n)=Aα2Δuejπ{cotα(m/2Δu)2-2cscα(mn/(2Δu)2)+cotα(n/2Δu)2}.In the same article, Ozaktas et al. [11] gave alternate methods to calculate the DFRFT by introducing scaling process. All these methods resulted in the same time complexity. However, in practice, these alternate methods are not preferable because they require coordinate scaling. Scaling is not an advantageous process if the original data has already been sampled. This is because scaling requires additional interpolations and so forth, which will require additional computation.However, in both cases, the Wigner distribution of the signal was confined to a circle of diameterΔu around the origin. Thus, there might be several discrete fractional Fourier transform matrices, some of which were more elegantly expressible than the others, which gave the same result within the accuracy of this approximation. All these matrices would yield results that were in increasingly better agreement as the signal energy contained in the circle was increasingly closer to the total energy.


## 3.1.2. Closed-Form Type DFRFT

Pei and Ding [12] had derived a different type of DFRFT which was neat in concept and reversible but lacks the additivity property. Due to the orientation of practical usage, two types of DFRFT were derived which are different in parameterization. The parameters of the first type have directly linked to the continuous FRFT and suit the applications of computing the continuous FRFT. On the other hand, the second type has the simpler parameters set and allows more elegant expression for the operator kernels and more is suitable for other applications of DFRFT, such as the filter design, pattern recognition, and phase retrieval. Considering the expression of continuous FRFT:
(14)Fα(u)=1-jcotα2πej(u2/2)cotα∫-∞∞e-jtucscαej(t2/2)cotαf(t)dt.
And to derive the DFRFT, sampling of f(t) and Fα(u) with the interval Δt and Δu needs to be performed as:
(15)y(n)=f(n*Δt),Yα(m)=Fα(m*Δu),
where n=-N,…,N and m=-M,…,M. By using this sampling, FRFT can be converted [13] as:
(16)Yα(m)=1-jcotα2πΔtej(m2Δu2/2)cotα×∑n=-NNe-jnmΔtΔucscαej(n2Δt2/2)cotαy(n).The above can be written as:(17)Yα(m)=∑n=-NNFα(m,n)y(n),
where
(18)Fα(m,n)=1-jcotα2πΔtej(m2Δu2/2)cotα×e-jnmΔtΔucscαej(n2Δt2/2)cotα.
And Yα(m) can be reversible only if the inverse transform of Fα(m,n), for all M≥N, would be Hermitian (conjugate and transpose), given [12] as:
(19)y(n)=∑m=-MMFα*(m,n)Yα(m)forM≥N.From (17) and (19),
(20)y(n)=∑m=-MM∑k=-NNFα*(m,n)Fα(m,k)y(k).After some manipulation,(21)y(n)=Δt22π|sinα|∑m=-MM∑k=-NN×ej((k2-n2)Δt2/2)cotαejm(n-k)ΔtΔucscαy(k).The summation form in (21) should be equal to δ(n-k), then
(22)ΔtΔu=S2πsinα2m+1,
where |S| is some integer prime to (2m+1); in this case, (18) becomes
(23)Fα(m,n)=1-jcotα2πΔtej(m2Δu2/2)cotαe-j(S2πnm/(2m+1))×ej(n2Δt2/2)cotα,∑m=-MM∑k=-NNFα*(m,n)Fα(m,k)y(k)=2M+12π|sinα|Δt2y(n).
After normalization, the transform matrix [12] can be given as:
(24)Fα(m,n)=|sinα|-jsgn(sinα)cosα2M+1ej(m2Δu2/2)cotα×e-j(sgn(sinα)2πnm/(2M+1))ej(n2Δt2/2)cotα.
The definition of DFRFT is further divided into two parts [12]. The first part was defined for “sinα>0” as:
(25)Yα(m)=sinα-jcosα2M+1ej(m2Δu2/2)cotα×∑n=-NNe-j(2πnm/(2M+1))ej(n2Δt2/2)cotαy(n).
And the second part was defined for “sinα<0” as:
(26)Yα(m)=-sinα+jcosα2M+1ej(m2Δu2/2)cotα×∑n=-NNej(2πnm/(2M+1))ej(n2Δt2/2)cotαy(n).
This definition of DFRFT possesses reversibility property and the periodicity property:
(27)Yα(-m)=Yα+π(m),Yα(m)=Yα+2π(m).The DFRFT of type 1 has a very important advantage in terms of its efficiency in calculating and implementing the DFRFT [12]. Due to two chirp multiplications and one FFT required for the implementation of this type of DFRFT, the total number of the multiplication operations required was “2P+(P/2)log2P,” for length of the output as P=2M+1.Subsequently, another type of DFRFT, named as type 2, was introduced by modifying the definitions obtained for discrete affine Fourier transform (DAFT) and presented as:(28)Yp(M)=12M+1ej(m2p/2)×∑n=-NNe±j(2πnm/(2M+1))ej(n2p/2)y(n),∀M≥N.For the DFRFT of type 2, the parameter “p” is used to control the variation of the chirp in frequency domain [13]. The DFRFT of type 2 also needed “2P+(P/2)log2P” multiplication operations. It does not follow the additivity property but observes the convertibility operation.These definitions of DFRFT are better to previous ones, that is, sampling type, linear combination type, and eigen-function type, as it require less number of multiplications, having less complexity. Also it has a closed form expression being rich in properties.


## 3.2. Linear Combination Type DFRFT

In [14–17], the discrete fractional Fourier transform was derived by using the linear combination of identity operation (F0), discrete Fourier transform (Fπ/2), time inverse operation (Fπ), and inverse discrete Fourier transform (F3π/2). In [15], the concept of DFRFT was given by linear combination of Lagrange interpolation polynomial of degree 3. For any function x(v)(29)x(v)=∑k=03ϕk(v)g(λk),
where
(30)λk=ejkπ/2,for0≤k≤3.Andϕk(v) is the Lagrange interpolation polynomial of degree 3 taking value “one” at v=λk and “zero” otherwise. It follows that
(31)ϕk(v)=14∑i=03uie-jkiπ/2.
Finally, by taking the principle nth-root of eigenvalue λk, the fractional transform can be given as:
(32)Fn=∑i=03Fiαi(n),
where
(33)αi(n)=14∑k=14e-jk(n-i)π/2,for0≤i≤3.Subsequently, in [14, 15], it had also been shown that the operator defined by (17) is unitary; satisfies angle additivity and angle multiplicity properties, and the operator is periodic with a fundamental period of four. One analogous approach of this method was presented in [17], in which the DFRFT was calculated by first chirping the signal in time domain, then taking its Fourier transform, and finally chirping again in frequency domain. Depending on the character of the chirped signal, the embedded FT can be one of the two known classes in discrete case—discrete-time and periodic discrete-time.For discrete-time DFRFT, the expression for DFRFT of signalx(t) was given as:
(34)Xα(f)=∑k=-∞∞Tx(kT)Kα(f,kT).
And its inverse transform as:
(35)x(kT)=∫0FpXα(f)Kα*(f,kT)df,
where Fp=2/Tcscα and Kα(t,u) is the kernel of FRFT. Similarly, for periodic discrete-time signal x(t) with period Tp, the expression for DFRFT was given as:
(36)Xα(nf)=∑k=-∞∞Tx(kT)Kα(nF,kT).
And its inverse transform as:
(37)x(kT)=∑n=0N-1FXα(nF)Kα*(nF,kT),
where Tp=NT and Fp=NF=2/Tcscα.In this type of DFRFT, the transform matrix is orthogonal and satisfies the additivity property along with the reversibility property. However, the main problem with this method is that the transform result is not matched to the continuous FRFT. Besides, it performed very similarly to the original Fourier transform or the identity operation but lost the important characteristic of fractionalization.


## 3.3. Eigenvector Type DFRFT

In this classification, the basic approach of eigenvalues and eigenvectors of the kernel matrix of the FRFT is considered to define the discrete FRFT. Here, the eigenvalues of the discrete Fourier transform (DFT) are considered first and then by obtaining the powers of these eigenvectors of DFT matrix, the corresponding eigenvalues and eigenvectors of DFRFT matrix are derived. This classification includes three different approaches, named as, eigenvector decomposition method, random type method, andV¯ method DFRFT.


### 3.3.1. Eigenvector Decomposition Type DFRFT

The authors had derived another type of discrete fractional Fourier transform as reported in [13, 18–25] by searching the eigenvectors and eigenvalues of the DFT matrix followed by computing the fractional power of the DFT matrix based on Hermite function. This type of DFRFT worked very similarly to the continuous FRFT, and it fulfills the properties of orthogonality, additivity, and reversibility. The fractional power of matrix could be calculated from its eigendecomposition and the power of eigenvalues. Unfortunately, there exist two types of ambiguity in deciding the fractional power of the DFT kernel matrix.(i)
Ambiguity in Deciding the Fractional Powers of Eigenvalues.The square roots of unity are “1” and “–1” from elementary mathematics. This indicates that there exists root ambiguity in deciding the fractional power of eigenvalues.(ii)
Ambiguity in Deciding the Eigenvectors of the DFT Kernel Matrix.The DFT eigenvectors constitute four major eigensubspaces; therefore, the choices for the DFT eigenvectors to construct the DFRFT kernel are multiple and not unique.Because of the above mentioned ambiguities, several DFRFT kernel matrices were possible which could obey the rotational properties. The idea for developing the DFRFT is to find the discrete form of the following:(38)Kα(t,u)=∑n=0∞e-jnαHn(t)Hn(u),
where α indicates the rotation angle associated with FRFT of signal and Hn(t) represents the nth-order normalized Hermite function with unit variance. The nth-order normalized Hermite function with variance σ is defined as:
(39)Hn(t)=12nn!σπhn(tσ)e-(t2/2σ2),
where, hn(·) represents the nth-order Hermite polynomial. However, the normalized Hermite functions with unitary variance Hn(·) indicate the eigenfunction of the FRFT. With the help of Hermite functions, the FRFT of a signal can be computed [24] as:
(40)Xα(u)=∫-∞∞x(t)Kα(t,u)dt=∑n=0∞Hn(u)e-jnα{∫-∞∞x(t)Hn(t)dt}.The FRFT can be interpreted as a weighting summation of Hermite functions, as given in (40). The weighting coefficients are obtained by multiplying the phase term e-jnα and the inner product of the input signal with the corresponding Hermite functions [24]. The development of this DFRFT was based on the eigendecomposition of the DFT kernel. The eigenvalues of DFT kernel “F” can be given as {1,-j,-1,j}, and corresponding to each eigenvalue one eigenvector can be evaluated. The continuous FRFT has a Hermite function with unitary variance as its eigenfunction. The corresponding eigenfunction property for the DFT can be given as:
(41)F2α/π[u^n]=e-jnαu^n,
where u^n represents the eigenvector of DFT corresponding to the nth-order discrete Hermite functions. In order to retain the eigenfunction properties, the unit variance Hermite function is sampled with a period Ts=2π/N. In the case of continuous FRFT, the terms of the Hermite functions are summed up from order zero to infinity. However, for the discrete case, only eigenvectors for the DFT Hermite eigenvectors can be added. The selection of the DFT Hermite eigenvectors is usually made from low to high orders, due to small approximation error of the low DFT Hermite eigenvectors [24]. Finally, the transform kernel of DFRFT can be defined as:
(42)F2α/π=U^D2α/πU^T={∑k=0N-1e-jkαu^ku^kT,∀N=4m+1{∑k=0N-2e-jkαu^ku^kT}+e-jNαu^Nu^NT,u^5NT∀N=4m,
where U^=[u^0|u^1|⋯u^N-1|] for N as odd and U^=[u^0|u^1|⋯u^N-2|u^N] for N is even, u^k represents the normalized eigenvector corresponding to the kth-order discrete Hermite function, and D can be defined as follows.ForN odd,
(43)D2α/π=[e-j00e-jα⋱0e-jα(N-1)].
And for N even,
(44)D2α/π=[e-j00⋱e-jα(N-2)0e-jαN]The convergence for the eigenvectors obtained from these matrixes is not so fast for the high-order Hermite functions. To ensure orthogonality of the DFT Hermite eigenvectorsu^k, the Gram-Schmidt algorithm (GSA), or the orthogonal procrustes algorithm (OPA) can be used [24]. The GSA minimizes the errors between the samples of the Hermite functions and the orthogonal DFT Hermite eigenvectors. On the other hand, the OPA minimizes the total errors between those samples. The main difference between the approach proposed in [24] and similar approaches proposed in [21–23] is found in the eigenvectors obtained. The eigenvectors obtained by method reported in [24] were discrete Mathieu function, although the Mathieu functions can converge to Hermite functions as obtained by the methods reported in [21–23]. Given the number of points N and the rotation angle α, the DFT Hermite eigenvectors can be computed followed by determining the eigenvalues of DFRFT. The computation of the DFRFT can be implemented only by a transform kernel matrix multiplication. The complexity of computing the DFRFT is O(N2) same as in the DFT case. By adjusting the rotational angles, the method for implementing the DFRFT of a signal x can be given as:
(45)F2α/π=U^D2α/πU^Tx=∑k=0N-1ake-jkαu^k.The coefficientsak’s are the inner products of signal and eigenvectors, and they can be computed in advance. If rotation angle changed then only the diagonal matrix D2α/π needs to be recomputed. Additionally, the authors [24] investigated the relationship between the FRFT and the DFRFT and found that for a sampling period equal to 2π/N, the DFRFT performs a circular rotation of the signal in the time-frequency plane. However, the DFRFT becomes an elliptical rotation in the continuous time-frequency plane for sampling periods different from 2π/N [24]. Therefore, for these elliptical rotations an angle modification and a postphase compensation in the DFRFT have been required to obtain results similar to the continuous FRFT [24]. This approach has been extended to the so-called multiple-parameter discrete fractional Fourier transform (MPDFRFT) [20, 25]. In fact, the MPDFRFT maintains all desired properties and reduces to the DFRFT when all of its order parameters are the same.A similar approach to [24] has been proposed in [19]. However, authors in [19] believe that the discrete time counterparts of the continuous time Hermite-Gaussians maintained the same properties because these discrete time counterparts exhibited better approximations than the other proposed approaches [18]. In order to resolve this issue about the approximation of Hermite-Gaussian functions, a nearly tridiagonal commuting matrix of the DFT and a corresponding version of the DFRFT were proposed in [13]. Most of the eigenvectors of this proposed nearly tridiagonal matrix result in a good approximation of the continuous Hermite-Gaussian functions by providing a smaller approximation error in comparison to the previous approaches.


### 3.3.2. Random Type DFRFT

Pei and Hsue [26] had introduced a modified method to calculate the DFRFT based on the random DFT eigenvectors and eigenvalues. In this method, a new commuting matrix with random DFT eigenvectors was constructed. Then, a random discrete fractional Fourier transform (RDFRFT) kernel matrix with random DFT eigenvectors and eigenvalues was developed. The magnitude and phase of the transformed output are found to be random in RDFRFT. This can be considered as a special feature of RDFRFT. Previously, the DFT was randomized to define the discrete random Fourier transform (DRFT) in [27] by taking random powers of eigenvalues of the DFT matrix. However, the eigenvectors of the DRFT were not random and were computed using the same method as proposed in [23]. Although, the phase of the DRFT was random, its magnitude was not random, as shown in [27]. The eigen-decomposition of DFT matrix “F” was given as:
(46)F=∑k=0N-1λkekekT.
The DFRFT with one order parameter “a” was defined by [25]:
(47)Fa=∑k=0N-1λkaekekT.The definition of RDFRFT was based on the definition of multiple-parameter DFRFT (MPDFRFT) [25], given as:
(48)Fa^=∑k=0N-1λkakekekT,
where the 1×N parameter vector a^ was defined as:
(49)a^=[a0,a1,…,aN-1].This order parameter of the MPDFRFT can be taken as independent random numbers to define an MPDFRFT with the following random eigenvalues [26], λkak, k=0,1,…,N-1, and random eigenvalues were given as:
(50)λkak={(e-j2π)ak,forλk=1(e-jπ/2)ak,forλk=-j(e-jπ)ak,forλk=-1(e-j3π/2)ak,forλk=j.The generation of RDFRFT is given in the following section. In this generation of RDFRFT, the DFT-commuting matrix with orthonormal random DFT eigenvectors is used.First defining aN×N real random Matrix “A” being “K” symmetric and generated by A=KAK, where K represents circular reversal matrix, which is given by:
(51)K=F2=F-2=[100JN-1],
where, JN-1 is  the (N-1)×(N-1) reversal matrix, having nonzero entries “ones” placed on the antidiagonal locations. Taking the Ksymmetric part “E” of the random matrix “A” as B={(A+KAK)/2}, with the help of this matrix, a random matrix “C” was generated as C=(B+BT)/2, where “T” denotes the matrix transpose operation, and finally a DFT commuting matrix “D” will be formed [27] as D=C+FCF-1. Subsequently, with the random DFT-commuting matrix “D” and the MPDFRFT, the random DFRFT (RDFRFT) with the 1×N parameter vector a^=[a0,a1,…,aN-1] of a signal “x” can be determined [26] as:
(52)X(a^,D)=FDa^*x={∑k=0N-1λkakrkrkT}*x,
where rk represents the orthonormal random DFT eigenvectors computed from “D” and λk represents the DFT eigenvalue corresponding to rk, given by (50).The RDFRFT kernel matrixFDa^ satisfied various properties of the FRFT like—additivity, FT convertibility, and angle additivity. The computational complexity of the -point RDFRFT was O(N2) same as that of the MPDFRFT, but RDFRFT has an advantage over MPDFRFT in terms of one more free parameter associated with it as parameter vector a^. The magnitude and phase responses of the transformed output, with the RDFRFT, are random because the eigenvectors and eigenvalues of the RDFRFT matrix are random [26], as shown in (52). This aspect establishes its advantageous role in security areas like—cryptography, encryption, and watermarking.


### 3.3.3. TheV¯ Method DFRFT

In this work [28], authors focus on a different approach to derive eigenvectors of the DFT matrix in order to build the DFRFT matrix that mimics the properties of its continuous counterpart. A straightforward and refined derivation of the eigenvectors of the DFT matrix is proposed without using any commuting matrices, but using only the factored form of a basic property of the DFT matrix stating that four consecutive Fourier transforms is the identity transform
(53)WN4=IN,
where WN4 and IN are the centered DFT (CDFT) and the identity matrices of order N, respectively. We utilize factorization of the CDFT matrix to obtain new matrices, whose columns are eigenvectors of the CDFT matrix. Thereafter, GSA was used to find the orthonormal eigenvectors obtained by factorization of the CDFT matrix. By employing a DFT-shift matrix K, eigenvectors of the ordinary-DFT matrix are obtained by using the CDFT matrix. Later on, the S and T matrices together are also used to boost the performance.In [24], Pei et al. compute the samples of the Hermite-Gauss functions and employ Gram-Schmidt algorithm (GSA) to orthogonalize these samples by projecting on the S matrix, which is equivalent to orthogonalizing the samples of the Hermite-Gauss functions without using the S matrix. Hanna et al. [29] claim to find the eigenvectors of the DFT matrix using a similar method. However, they use the ordinary-DFT matrix as a basis to their work, but Serbes and Durak-Ata [28] employs the CDFT matrix. Additionally, the eigenvectors are nonorthogonal and dissimilar to the samples of Hermite-Gauss functions. Serbes and Durak-Ata [28] work is different from these methods in the sense that the non-orthogonal eigenvectors of the CDFT matrix is determined without utilizing any S matrix or samples of the Hermite-Gauss functions, and then orthogonalization of eigenvectors is obtained by employing the GSA.The DFT maps the signal from[0,N-1] discrete input space to [0,2π] discrete frequency space, whereas the CDFT maps [-(N-1)/2,(N-1)/2] to [-π,π]. This allows CDFT to define even and odd functions such as Hermite-Gauss functions, that is contrary to the ordinary DFT definition. In order to approximate the samples of the continuous FrFT and to imply the rotation property, Hermite-Gauss-like eigenvectors of the DFT matrix have to be obtained. Serbes and Durak-Ata [28] use the CDFT matrix instead of the ordinary DFT matrix, because the Hermite-Gauss-like vectors are eigenvectors of only the CDFT matrix, not the DFT matrix, that is, when a discrete Hermite-Gauss-like eigenvector is transformed using the ordinary DFT matrix, the output is the shifted version of the Hermite-Gauss vector multiplied by a complex sinusoidal.The FRFT operator of orderα can be defined as the αth-power of the ordinary DFT operator WN. Hence, the DFRFT matrix can be expressed by means of its eigenvector decomposition,
(54)WNα=UNΛNαUNT,
where ΛNα is given as—ΛNα=diag(e-j0,e-jπ/2a,…,e-jπ/2a(N-2),e-jπ/2a(N-1)) for the centered-DFRFT.The centered-DFRFT matrix can be obtained by calculating the Hermite-Gaussian-like eigenvectors of DFRFT as followes:(55)WNα=V1¯Λ¯k1αV¯1T+V2¯Λ¯k2αV¯2T+V3¯Λ¯k3αV¯3T+V4¯Λ¯k4αV¯4T,
where
(56)Λ¯k1α=diag(e-j2π(k1-1)a,…,e-j2πa,e-j0),Λ¯k2α=diag(e-j(π+(k2-1)2π)a,…,e-j3πa,e-jπa),Λ¯k3α=diag(e-j(3π/2+(k3-1)2π)a,…,e-j7πa/2,e-j3πa/2),Λ¯k4α=diag(e-j(π/2+(k4-1)2π)a,…,e-j5πa/2,e-jπa/2).The above mentioned method satisfies unitary and angle additivity properties. Also, it reduces to the ordinary DFT whena=1 and approximates to the samples of the continuous FRFT for fractional values of a.Performance of the Serbes and Durak-Ata [28] algorithm can be improved by using it in combination with S and T matrices. Pei et al. [13] uses the S method together with the Grunbaum’s T matrix [30] as a linear combination, in which it has been observed that Hermite-Gauss-like eigenvectors of (S+15T) are more accurate than both S and T alone. Serbes and Durak-Ata [28] proposed that eigenvectors of linear combinations of S, T, and VT matrices as k1S+k1T+k1VT produce more accurate eigenvectors. Pei and Hsue [25] showed that using the linear combination of S and T as (S+15T) produces better results than the S or T alone. Whereas, Serbes and Durak-Ata [28] showed that using (S+30T-7VT) produces a more accurate approximation to the continuous FRFT. Although, better commuting matrix generation method is available in the literature [31–33], but all these methods are explained in the context of DFT and some other transform only, excluding DFRFT. Therefore, the betterment offered by these methods can be utilized in devising the new algorithm of DFRFT based on these commuting matrixes.


## 3.3.1. Eigenvector Decomposition Type DFRFT

The authors had derived another type of discrete fractional Fourier transform as reported in [13, 18–25] by searching the eigenvectors and eigenvalues of the DFT matrix followed by computing the fractional power of the DFT matrix based on Hermite function. This type of DFRFT worked very similarly to the continuous FRFT, and it fulfills the properties of orthogonality, additivity, and reversibility. The fractional power of matrix could be calculated from its eigendecomposition and the power of eigenvalues. Unfortunately, there exist two types of ambiguity in deciding the fractional power of the DFT kernel matrix.(i)
Ambiguity in Deciding the Fractional Powers of Eigenvalues.The square roots of unity are “1” and “–1” from elementary mathematics. This indicates that there exists root ambiguity in deciding the fractional power of eigenvalues.(ii)
Ambiguity in Deciding the Eigenvectors of the DFT Kernel Matrix.The DFT eigenvectors constitute four major eigensubspaces; therefore, the choices for the DFT eigenvectors to construct the DFRFT kernel are multiple and not unique.Because of the above mentioned ambiguities, several DFRFT kernel matrices were possible which could obey the rotational properties. The idea for developing the DFRFT is to find the discrete form of the following:(38)Kα(t,u)=∑n=0∞e-jnαHn(t)Hn(u),
where α indicates the rotation angle associated with FRFT of signal and Hn(t) represents the nth-order normalized Hermite function with unit variance. The nth-order normalized Hermite function with variance σ is defined as:
(39)Hn(t)=12nn!σπhn(tσ)e-(t2/2σ2),
where, hn(·) represents the nth-order Hermite polynomial. However, the normalized Hermite functions with unitary variance Hn(·) indicate the eigenfunction of the FRFT. With the help of Hermite functions, the FRFT of a signal can be computed [24] as:
(40)Xα(u)=∫-∞∞x(t)Kα(t,u)dt=∑n=0∞Hn(u)e-jnα{∫-∞∞x(t)Hn(t)dt}.The FRFT can be interpreted as a weighting summation of Hermite functions, as given in (40). The weighting coefficients are obtained by multiplying the phase term e-jnα and the inner product of the input signal with the corresponding Hermite functions [24]. The development of this DFRFT was based on the eigendecomposition of the DFT kernel. The eigenvalues of DFT kernel “F” can be given as {1,-j,-1,j}, and corresponding to each eigenvalue one eigenvector can be evaluated. The continuous FRFT has a Hermite function with unitary variance as its eigenfunction. The corresponding eigenfunction property for the DFT can be given as:
(41)F2α/π[u^n]=e-jnαu^n,
where u^n represents the eigenvector of DFT corresponding to the nth-order discrete Hermite functions. In order to retain the eigenfunction properties, the unit variance Hermite function is sampled with a period Ts=2π/N. In the case of continuous FRFT, the terms of the Hermite functions are summed up from order zero to infinity. However, for the discrete case, only eigenvectors for the DFT Hermite eigenvectors can be added. The selection of the DFT Hermite eigenvectors is usually made from low to high orders, due to small approximation error of the low DFT Hermite eigenvectors [24]. Finally, the transform kernel of DFRFT can be defined as:
(42)F2α/π=U^D2α/πU^T={∑k=0N-1e-jkαu^ku^kT,∀N=4m+1{∑k=0N-2e-jkαu^ku^kT}+e-jNαu^Nu^NT,u^5NT∀N=4m,
where U^=[u^0|u^1|⋯u^N-1|] for N as odd and U^=[u^0|u^1|⋯u^N-2|u^N] for N is even, u^k represents the normalized eigenvector corresponding to the kth-order discrete Hermite function, and D can be defined as follows.ForN odd,
(43)D2α/π=[e-j00e-jα⋱0e-jα(N-1)].
And for N even,
(44)D2α/π=[e-j00⋱e-jα(N-2)0e-jαN]The convergence for the eigenvectors obtained from these matrixes is not so fast for the high-order Hermite functions. To ensure orthogonality of the DFT Hermite eigenvectorsu^k, the Gram-Schmidt algorithm (GSA), or the orthogonal procrustes algorithm (OPA) can be used [24]. The GSA minimizes the errors between the samples of the Hermite functions and the orthogonal DFT Hermite eigenvectors. On the other hand, the OPA minimizes the total errors between those samples. The main difference between the approach proposed in [24] and similar approaches proposed in [21–23] is found in the eigenvectors obtained. The eigenvectors obtained by method reported in [24] were discrete Mathieu function, although the Mathieu functions can converge to Hermite functions as obtained by the methods reported in [21–23]. Given the number of points N and the rotation angle α, the DFT Hermite eigenvectors can be computed followed by determining the eigenvalues of DFRFT. The computation of the DFRFT can be implemented only by a transform kernel matrix multiplication. The complexity of computing the DFRFT is O(N2) same as in the DFT case. By adjusting the rotational angles, the method for implementing the DFRFT of a signal x can be given as:
(45)F2α/π=U^D2α/πU^Tx=∑k=0N-1ake-jkαu^k.The coefficientsak’s are the inner products of signal and eigenvectors, and they can be computed in advance. If rotation angle changed then only the diagonal matrix D2α/π needs to be recomputed. Additionally, the authors [24] investigated the relationship between the FRFT and the DFRFT and found that for a sampling period equal to 2π/N, the DFRFT performs a circular rotation of the signal in the time-frequency plane. However, the DFRFT becomes an elliptical rotation in the continuous time-frequency plane for sampling periods different from 2π/N [24]. Therefore, for these elliptical rotations an angle modification and a postphase compensation in the DFRFT have been required to obtain results similar to the continuous FRFT [24]. This approach has been extended to the so-called multiple-parameter discrete fractional Fourier transform (MPDFRFT) [20, 25]. In fact, the MPDFRFT maintains all desired properties and reduces to the DFRFT when all of its order parameters are the same.A similar approach to [24] has been proposed in [19]. However, authors in [19] believe that the discrete time counterparts of the continuous time Hermite-Gaussians maintained the same properties because these discrete time counterparts exhibited better approximations than the other proposed approaches [18]. In order to resolve this issue about the approximation of Hermite-Gaussian functions, a nearly tridiagonal commuting matrix of the DFT and a corresponding version of the DFRFT were proposed in [13]. Most of the eigenvectors of this proposed nearly tridiagonal matrix result in a good approximation of the continuous Hermite-Gaussian functions by providing a smaller approximation error in comparison to the previous approaches.


## 3.3.2. Random Type DFRFT

Pei and Hsue [26] had introduced a modified method to calculate the DFRFT based on the random DFT eigenvectors and eigenvalues. In this method, a new commuting matrix with random DFT eigenvectors was constructed. Then, a random discrete fractional Fourier transform (RDFRFT) kernel matrix with random DFT eigenvectors and eigenvalues was developed. The magnitude and phase of the transformed output are found to be random in RDFRFT. This can be considered as a special feature of RDFRFT. Previously, the DFT was randomized to define the discrete random Fourier transform (DRFT) in [27] by taking random powers of eigenvalues of the DFT matrix. However, the eigenvectors of the DRFT were not random and were computed using the same method as proposed in [23]. Although, the phase of the DRFT was random, its magnitude was not random, as shown in [27]. The eigen-decomposition of DFT matrix “F” was given as:
(46)F=∑k=0N-1λkekekT.
The DFRFT with one order parameter “a” was defined by [25]:
(47)Fa=∑k=0N-1λkaekekT.The definition of RDFRFT was based on the definition of multiple-parameter DFRFT (MPDFRFT) [25], given as:
(48)Fa^=∑k=0N-1λkakekekT,
where the 1×N parameter vector a^ was defined as:
(49)a^=[a0,a1,…,aN-1].This order parameter of the MPDFRFT can be taken as independent random numbers to define an MPDFRFT with the following random eigenvalues [26], λkak, k=0,1,…,N-1, and random eigenvalues were given as:
(50)λkak={(e-j2π)ak,forλk=1(e-jπ/2)ak,forλk=-j(e-jπ)ak,forλk=-1(e-j3π/2)ak,forλk=j.The generation of RDFRFT is given in the following section. In this generation of RDFRFT, the DFT-commuting matrix with orthonormal random DFT eigenvectors is used.First defining aN×N real random Matrix “A” being “K” symmetric and generated by A=KAK, where K represents circular reversal matrix, which is given by:
(51)K=F2=F-2=[100JN-1],
where, JN-1 is  the (N-1)×(N-1) reversal matrix, having nonzero entries “ones” placed on the antidiagonal locations. Taking the Ksymmetric part “E” of the random matrix “A” as B={(A+KAK)/2}, with the help of this matrix, a random matrix “C” was generated as C=(B+BT)/2, where “T” denotes the matrix transpose operation, and finally a DFT commuting matrix “D” will be formed [27] as D=C+FCF-1. Subsequently, with the random DFT-commuting matrix “D” and the MPDFRFT, the random DFRFT (RDFRFT) with the 1×N parameter vector a^=[a0,a1,…,aN-1] of a signal “x” can be determined [26] as:
(52)X(a^,D)=FDa^*x={∑k=0N-1λkakrkrkT}*x,
where rk represents the orthonormal random DFT eigenvectors computed from “D” and λk represents the DFT eigenvalue corresponding to rk, given by (50).The RDFRFT kernel matrixFDa^ satisfied various properties of the FRFT like—additivity, FT convertibility, and angle additivity. The computational complexity of the -point RDFRFT was O(N2) same as that of the MPDFRFT, but RDFRFT has an advantage over MPDFRFT in terms of one more free parameter associated with it as parameter vector a^. The magnitude and phase responses of the transformed output, with the RDFRFT, are random because the eigenvectors and eigenvalues of the RDFRFT matrix are random [26], as shown in (52). This aspect establishes its advantageous role in security areas like—cryptography, encryption, and watermarking.


## 3.3.3. TheV¯ Method DFRFT

In this work [28], authors focus on a different approach to derive eigenvectors of the DFT matrix in order to build the DFRFT matrix that mimics the properties of its continuous counterpart. A straightforward and refined derivation of the eigenvectors of the DFT matrix is proposed without using any commuting matrices, but using only the factored form of a basic property of the DFT matrix stating that four consecutive Fourier transforms is the identity transform
(53)WN4=IN,
where WN4 and IN are the centered DFT (CDFT) and the identity matrices of order N, respectively. We utilize factorization of the CDFT matrix to obtain new matrices, whose columns are eigenvectors of the CDFT matrix. Thereafter, GSA was used to find the orthonormal eigenvectors obtained by factorization of the CDFT matrix. By employing a DFT-shift matrix K, eigenvectors of the ordinary-DFT matrix are obtained by using the CDFT matrix. Later on, the S and T matrices together are also used to boost the performance.In [24], Pei et al. compute the samples of the Hermite-Gauss functions and employ Gram-Schmidt algorithm (GSA) to orthogonalize these samples by projecting on the S matrix, which is equivalent to orthogonalizing the samples of the Hermite-Gauss functions without using the S matrix. Hanna et al. [29] claim to find the eigenvectors of the DFT matrix using a similar method. However, they use the ordinary-DFT matrix as a basis to their work, but Serbes and Durak-Ata [28] employs the CDFT matrix. Additionally, the eigenvectors are nonorthogonal and dissimilar to the samples of Hermite-Gauss functions. Serbes and Durak-Ata [28] work is different from these methods in the sense that the non-orthogonal eigenvectors of the CDFT matrix is determined without utilizing any S matrix or samples of the Hermite-Gauss functions, and then orthogonalization of eigenvectors is obtained by employing the GSA.The DFT maps the signal from[0,N-1] discrete input space to [0,2π] discrete frequency space, whereas the CDFT maps [-(N-1)/2,(N-1)/2] to [-π,π]. This allows CDFT to define even and odd functions such as Hermite-Gauss functions, that is contrary to the ordinary DFT definition. In order to approximate the samples of the continuous FrFT and to imply the rotation property, Hermite-Gauss-like eigenvectors of the DFT matrix have to be obtained. Serbes and Durak-Ata [28] use the CDFT matrix instead of the ordinary DFT matrix, because the Hermite-Gauss-like vectors are eigenvectors of only the CDFT matrix, not the DFT matrix, that is, when a discrete Hermite-Gauss-like eigenvector is transformed using the ordinary DFT matrix, the output is the shifted version of the Hermite-Gauss vector multiplied by a complex sinusoidal.The FRFT operator of orderα can be defined as the αth-power of the ordinary DFT operator WN. Hence, the DFRFT matrix can be expressed by means of its eigenvector decomposition,
(54)WNα=UNΛNαUNT,
where ΛNα is given as—ΛNα=diag(e-j0,e-jπ/2a,…,e-jπ/2a(N-2),e-jπ/2a(N-1)) for the centered-DFRFT.The centered-DFRFT matrix can be obtained by calculating the Hermite-Gaussian-like eigenvectors of DFRFT as followes:(55)WNα=V1¯Λ¯k1αV¯1T+V2¯Λ¯k2αV¯2T+V3¯Λ¯k3αV¯3T+V4¯Λ¯k4αV¯4T,
where
(56)Λ¯k1α=diag(e-j2π(k1-1)a,…,e-j2πa,e-j0),Λ¯k2α=diag(e-j(π+(k2-1)2π)a,…,e-j3πa,e-jπa),Λ¯k3α=diag(e-j(3π/2+(k3-1)2π)a,…,e-j7πa/2,e-j3πa/2),Λ¯k4α=diag(e-j(π/2+(k4-1)2π)a,…,e-j5πa/2,e-jπa/2).The above mentioned method satisfies unitary and angle additivity properties. Also, it reduces to the ordinary DFT whena=1 and approximates to the samples of the continuous FRFT for fractional values of a.Performance of the Serbes and Durak-Ata [28] algorithm can be improved by using it in combination with S and T matrices. Pei et al. [13] uses the S method together with the Grunbaum’s T matrix [30] as a linear combination, in which it has been observed that Hermite-Gauss-like eigenvectors of (S+15T) are more accurate than both S and T alone. Serbes and Durak-Ata [28] proposed that eigenvectors of linear combinations of S, T, and VT matrices as k1S+k1T+k1VT produce more accurate eigenvectors. Pei and Hsue [25] showed that using the linear combination of S and T as (S+15T) produces better results than the S or T alone. Whereas, Serbes and Durak-Ata [28] showed that using (S+30T-7VT) produces a more accurate approximation to the continuous FRFT. Although, better commuting matrix generation method is available in the literature [31–33], but all these methods are explained in the context of DFT and some other transform only, excluding DFRFT. Therefore, the betterment offered by these methods can be utilized in devising the new algorithm of DFRFT based on these commuting matrixes.


## 3.4. Weighted Summation Type DFRFT

Yeh and Pei [34] had developed a new method for DFRFT computation. The idea of the developed method was to compute the DFRFT at any angle by a weighted summation of the DFRFT’s with the special angles. In this method DFRFT computation with odd point of length can be realized by the weighted summation of the DFRFTs in special angles. The special angles were multiples of 2π/N for the odd case and multiples of 2π/(N+1) for even length. Regardless of even- or odd-length cases, the weighting coefficients were obtained from an IDFT operation. The N-point IDFT was needed for odd length, and the (N+1) - point was computed for even length.For discrete signal “x” having odd length “N”, the DFRFT of “x” for rotation angle “α” can be computed as:
(57)Xα=∑n=0N-1Bn,αXn,β∀β=2πN,
where the weighting coefficients Bn,α was given as:
(58)Bn,α=IDFT[e-jkα]∀k=0,1,…,N-1=1N∑n=0N-1e-jkαej(2π/N)nk.
The weighted summation Bn,α has a closed-form solution:
(59)Bn,α={1N1-e-j(N-1)(α-nβ)1-e-j(α-nβ),∀α≠kβδ(n-k),∀α=kβ.Similarly, for discrete signal “x” having even length N, the DFRFT of “x” for rotation angle “α” can be computed by the following expression:
(60)Xα=∑n=0NBn,αXn,β∀β=2πN+1.
The weighting coefficients, Bn,α can be computed as:
(61)Bn,α=IDFT[e-jkα]∀k=0,1,…,N=1N+1∑n=0Ne-jkαej(2π/(N+1))nk.
Similar to the odd length, the weighted summation Bn,α for even case has also a closed-form solution:
(62)Bn,α={1N+11-e-jN(α-nβ)1-e-j(α-nβ)∀α≠kβδ(n-k)∀α=kβThe eigenvector decomposition methods requireN2 storages to store the transform kernel. If the computing angle is changed, the transform kernel is also changed and needs to be recomputed. In this algorithm, the DFRFT of the specified angle was stored. The DFRFT at any angle can be obtained by weighted summation of these specified DFRFTs. The weighted computation in this algorithm still takes O(N2) multiplications; therefore, its computation load is still O(N2). In [34], two implementation methods for the DFRFT computation algorithm were introduced. One was called the parallel method and the other cascade method. Similar to the parallel method, the special angle and numbers of terms of cascade forms are also different for the even and odd cases. Both of the implementation methods can have the advantages over each other for different application. The parallel method was suitable for the signal, whose DFRFTs with special angles were already known. Hence, the computation of the DFRFT will become only a linear combination of the DFRFT’s in special angles. And, the cascade method means that the computation of the DFRFT can be realized from the DFRFT with only one specified angle. If the DFRFT with the specified angle can be computed efficiently, the computation of the DFRFT will become efficient. Therefore, it has been found suitable for VLSI implementation [34].Although, the existence of one well-accepted closed-form definition of DFRFT is not available; however many definitions with their pros and cons are accessible. A comparison of all the proposed classes of DFRFT based on some properties and complexity is shown in Table1. In this Table, yes is included against the method which possesses given property/parameter, otherwise NO is mentioned.Table 1
Comparison of Different Types of DFRFT Algorithms.
Properties/parameters
Sampling-based method
Linear combination-based method
Eigenvector-based method
Weighted summation-based method
Sampling type
Closed-form type
Linearcombination type
Eigenvector decomposition type
Random type
V
¯ Method type
Weighted summation Type
Reversibility
NO
YES
YES
YES
YES
YES
YES
Additivity
NO
NO
YES
YES
YES
YES
YES
Periodicity
YES
YES
YES
YES
YES
YES
YES
Symmetry
YES
YES
NO
YES
YES
YES
YES
Angle additivity
NO
YES
NO
YES
YES
YES
YES
Closed-form expression
YES
YES
YES
YES
YES
YES
YES
Similarity with FRFT
YES
YES
NO
YES
YES
YES
NO
Complexity
O(N log2N)
O(N log2N)
O(N log2N)
O(N2)
O
(
N
2
)
O(N log2N)
O(N2)
Constraints
Many
Fewer
Many
Many
Very less
Many
FewerFrom the above mentioned comparative analysis, it can be recapitulated that RDFRFT provides a better solution for computation of the DFRFT. The computation of RDFRFT requires more time than sampling type, linear combination type, and closed form type; however, it satisfies more properties than any other method with very less number of imposed constraints. Previously, IEEE has conferred to use the Eigenvectors decomposition type DFRFT as the most appropriate definition of DFRFT. In the year 2006, Pei et al. [25] had established the superiority of their method (MPDFRFT) over the earlier reported methods. Subsequently, in the year 2009, the MPDFRFT was modified by Pei et al. [26] and given as RDFRFT. It was also established that RDFRFT has one more degree of freedom, while satisfying the same set of properties as satisfied by MPDFRFT with equal computational complexity. Although, Serbes and Durak-Ata [28] proposed a modified version of eigenvalue decomposition technique with slight betterment but this does not redeem the advantage offered by RDFRFT, that is, security of the data. This makes the RDFRFT as a better choice in security applications.


## 4. Application of DFRFT

In this section, the application of DFRFT has been investigated in two areas. First, the detection of Doppler shifted DSB-SC signal will be performed than DFRFT has been applied in the image processing application, particularly in steganography.


### 4.1. Removal of Doppler Shift in DSB-SC AM Signal

The weighted autocorrelation theorem defined for FRFT [2] is utilized to remove the unwanted frequency components introduced in double sideband suppressed-carrier amplitude modulated (DSB-SC-AM) signal. It is assumed that a Doppler component is introduced in the carrier frequency due to motion of transmitter and/or receiver. The transmitted DSB-SC amplitude-modulated signal represented by “e(t)” with modulating signal frequency as fm and carrier frequency as fc is given as:
(63)e(t)=cos{2πfmt}cos{2πfct}.Similarly, the received signal, which is corrupted by noise (Doppler frequency component—fd), represented as “r(t)” is given by:
(64)r(t)=cos{2πfmt}cos{2π(fc+fd)t}.The values of “fm, fc, and fd” are taken as 1 Hz, 10 Hz and as 10% of carrier frequency, respectively. The simulation is performed on the platform of MATLAB software (version 7.1) on a system having configuration: processor as Pentium-4, Intel (R) CPU-1.8 GHz processor having 1 GB RAM. The effect of Doppler shift in the carrier is illustrated in Figure 2(a), where Fourier transform of original signal and Doppler corrupted signal is plotted.Removal of Doppler component introduced in DSB-SC signal. (a) FRFT of weighted autocorrelation of DSB-SC signal and corrupted DSB-SC signal at the FRFT order ofa=1 (i.e., Fourier transform domain) and (b) FRFT of weighted auto-correlation of corrupted DSB-SC signal at different FRFT order.
(a)
(b)It is evident from Figure2(a) that for a received Doppler corrupted DSB-SC signal, the magnitude spectrum of the signal in Fourier domain is not giving the exact information about original frequency component present in the transmitted signal due to only one peak present in the spectrum in Fourier domain. Whereas, when the same corrupted DSB-SC signal is analyzed in fractional Fourier domain, the magnitude spectrum of FRFT of weighted auto-correlation of corrupted DSB-SC signal gives two clear peak at an angle corresponding to the value of a=0.16, as shown in Figure 2(b). Subsequently, the actual frequency components present in transmitted DSB-SC can be estimated by observing the location of two peaks in this FRFT domain. From the observed value of FRFT domain variable “upeak” at which the peak is obtained, the frequency component present in the signal can be estimated by:
(65)Estimatedvalue=upeakcos{(1-a)(π/2)}.In Figure2(b), two simulated values of upeak is obtained as upeak,low and upeak,high for a=0.16. Then these values of upeak is used in (65) to obtain the desired frequency components. The assumed frequency component value and the estimated ones are shown in Table 2.Table 2
Frequency estimation in DSB-SC signal.
Assume value
Calculated value
Simulated value
Estimated value
f
m
f
c
f
c
+
f
m
f
c
-
f
m
u
peak
,
high
u
peak
,
low
f
c
+
f
m
f
c
-
f
m
1
10
11
9
2.769
2.22
11.13
8.927From Table2, it can be clearly observed that the calculated value of higher and lower frequency components present in DSB-SC AM signal is matched nicely with the estimated value of higher and lower frequency components even when the signal is corrupted with the presence of a Doppler frequency component in the signal itself.


### 4.2. Data Hiding by Steganography in Fractional Fourier Domain

Steganography is the technique of hiding confidential information within any other media. The objective of steganography is to hide a secret message within a cover-media in such a way that others cannot determine the presence of the hidden message. Hiding information into a media requires the following elements:(i)
the cover media (C) that will hold the hidden data,(ii)
the secret message (M), may be text, image, or any type of data,(iii)
the stego function and its inverse,(iv)
stego-key (K) may be used to hide and unhide the message.In this exercise, the cover and secret messages are assumed as Lena image and Mandrill image, respectively. The stego function is defined as(66)StegoImage=ℱa-1[ℱa{CoverImage}+b*ℱa{SecretImage}].Here, first, the DFRFT (ℱa) of cover image and secret image is obtained, then the steganographed image can be achieved by using the stego function defined in (66), where ℱa-1 is representing the inverse DFRFT and b is used as stego weight. Similarly, the inverse stego function for obtaining the secret information can be written as:
(67)Secretimage=ℱa1-1[ℱa1{StegoImage}-ℱa1{SecretImage}b].In this process, the values of stego weight “b” and DFRFT order “a” have been assumed as 0.05 and 0.2, respectively. For retrieval of secret image, the DFRFT order “a1” in (67) should be equal to the negative of the DFRFT order “a” in (66). The results obtained are shown in Figure 3. In this simulation exercise, the detection of hidden image is also performed for two other DFRFT angle parameters, both are in other side of actual, that is, at a1=-0.4 and a1=-0.22 in Figures 3(e) and 3(f), respectively. It can be easily verified from these results that for obtaining the hidden image correctly, the correct value of angular parameter of DFRFT should be known; hence, this correct value can also be used as a key other than stego weight. This extra key is only possible with the DFRFT and not with DFT, as there no such variation of angle is possible with DFT. This again confirms the superiority of DFRFT over DFT in the security applications.Steganography using DFRFT.
(a)
Cover image
(b)
Secret image
(c)
Steganography Image
(d)
Retrieved secret image ata1 = −0.2
(e)
Retrieved secret image ata1 = −0.4
(f)
Retrieved secret image ata1 = −0.22


## 4.1. Removal of Doppler Shift in DSB-SC AM Signal

The weighted autocorrelation theorem defined for FRFT [2] is utilized to remove the unwanted frequency components introduced in double sideband suppressed-carrier amplitude modulated (DSB-SC-AM) signal. It is assumed that a Doppler component is introduced in the carrier frequency due to motion of transmitter and/or receiver. The transmitted DSB-SC amplitude-modulated signal represented by “e(t)” with modulating signal frequency as fm and carrier frequency as fc is given as:
(63)e(t)=cos{2πfmt}cos{2πfct}.Similarly, the received signal, which is corrupted by noise (Doppler frequency component—fd), represented as “r(t)” is given by:
(64)r(t)=cos{2πfmt}cos{2π(fc+fd)t}.The values of “fm, fc, and fd” are taken as 1 Hz, 10 Hz and as 10% of carrier frequency, respectively. The simulation is performed on the platform of MATLAB software (version 7.1) on a system having configuration: processor as Pentium-4, Intel (R) CPU-1.8 GHz processor having 1 GB RAM. The effect of Doppler shift in the carrier is illustrated in Figure 2(a), where Fourier transform of original signal and Doppler corrupted signal is plotted.Removal of Doppler component introduced in DSB-SC signal. (a) FRFT of weighted autocorrelation of DSB-SC signal and corrupted DSB-SC signal at the FRFT order ofa=1 (i.e., Fourier transform domain) and (b) FRFT of weighted auto-correlation of corrupted DSB-SC signal at different FRFT order.
(a)
(b)It is evident from Figure2(a) that for a received Doppler corrupted DSB-SC signal, the magnitude spectrum of the signal in Fourier domain is not giving the exact information about original frequency component present in the transmitted signal due to only one peak present in the spectrum in Fourier domain. Whereas, when the same corrupted DSB-SC signal is analyzed in fractional Fourier domain, the magnitude spectrum of FRFT of weighted auto-correlation of corrupted DSB-SC signal gives two clear peak at an angle corresponding to the value of a=0.16, as shown in Figure 2(b). Subsequently, the actual frequency components present in transmitted DSB-SC can be estimated by observing the location of two peaks in this FRFT domain. From the observed value of FRFT domain variable “upeak” at which the peak is obtained, the frequency component present in the signal can be estimated by:
(65)Estimatedvalue=upeakcos{(1-a)(π/2)}.In Figure2(b), two simulated values of upeak is obtained as upeak,low and upeak,high for a=0.16. Then these values of upeak is used in (65) to obtain the desired frequency components. The assumed frequency component value and the estimated ones are shown in Table 2.Table 2
Frequency estimation in DSB-SC signal.
Assume value
Calculated value
Simulated value
Estimated value
f
m
f
c
f
c
+
f
m
f
c
-
f
m
u
peak
,
high
u
peak
,
low
f
c
+
f
m
f
c
-
f
m
1
10
11
9
2.769
2.22
11.13
8.927From Table2, it can be clearly observed that the calculated value of higher and lower frequency components present in DSB-SC AM signal is matched nicely with the estimated value of higher and lower frequency components even when the signal is corrupted with the presence of a Doppler frequency component in the signal itself.


## 4.2. Data Hiding by Steganography in Fractional Fourier Domain

Steganography is the technique of hiding confidential information within any other media. The objective of steganography is to hide a secret message within a cover-media in such a way that others cannot determine the presence of the hidden message. Hiding information into a media requires the following elements:(i)
the cover media (C) that will hold the hidden data,(ii)
the secret message (M), may be text, image, or any type of data,(iii)
the stego function and its inverse,(iv)
stego-key (K) may be used to hide and unhide the message.In this exercise, the cover and secret messages are assumed as Lena image and Mandrill image, respectively. The stego function is defined as(66)StegoImage=ℱa-1[ℱa{CoverImage}+b*ℱa{SecretImage}].Here, first, the DFRFT (ℱa) of cover image and secret image is obtained, then the steganographed image can be achieved by using the stego function defined in (66), where ℱa-1 is representing the inverse DFRFT and b is used as stego weight. Similarly, the inverse stego function for obtaining the secret information can be written as:
(67)Secretimage=ℱa1-1[ℱa1{StegoImage}-ℱa1{SecretImage}b].In this process, the values of stego weight “b” and DFRFT order “a” have been assumed as 0.05 and 0.2, respectively. For retrieval of secret image, the DFRFT order “a1” in (67) should be equal to the negative of the DFRFT order “a” in (66). The results obtained are shown in Figure 3. In this simulation exercise, the detection of hidden image is also performed for two other DFRFT angle parameters, both are in other side of actual, that is, at a1=-0.4 and a1=-0.22 in Figures 3(e) and 3(f), respectively. It can be easily verified from these results that for obtaining the hidden image correctly, the correct value of angular parameter of DFRFT should be known; hence, this correct value can also be used as a key other than stego weight. This extra key is only possible with the DFRFT and not with DFT, as there no such variation of angle is possible with DFT. This again confirms the superiority of DFRFT over DFT in the security applications.Steganography using DFRFT.
(a)
Cover image
(b)
Secret image
(c)
Steganography Image
(d)
Retrieved secret image ata1 = −0.2
(e)
Retrieved secret image ata1 = −0.4
(f)
Retrieved secret image ata1 = −0.22


## 5. Conclusion

In this paper, the different approaches of obtaining discrete fractional Fourier transform (DFRFT) have been analyzed. These different approaches are classified in four categories based on their methodology of evaluation. Further, with the help of a comparative analysis, it is established that the random type DFRFT emerged out as a better choice for calculating DFRFT. Thereafter, the random type DFRFT along with weighted auto-correlation theorem defined for the FRFT has been employed for the successful removal of the Doppler frequency component introduced in the DSB-SC AM signal due to time-variant motion of transmitter or receiver. Subsequently, the role of DFRFT in steganography has also been successfully performed with the availability of one extra key, in the form of DFRFT angular parameter. The results obtained are in the conformity of the superiority of the fractional Fourier transform over Fourier transform.



---

*Source: 214650-2013-04-29.xml*