Discrete Pseudo-Fractional Fourier Transform and Its Fast Algorithm

In this article, we introduce a new discrete fractional transform for data sequences whose size is a composite number. The main kernels of the introduced transform are small-size discrete fractional Fourier transforms. Since the introduced transformation is not, in the generally known sense, a classical discrete fractional transform, we call it discrete pseudo-fractional Fourier transform. We also provide a generalization of this new transform, which depends on many fractional parameters. A fast algorithm for computing the introduced transform is developed and described.


Introduction
Discrete fractional transforms are the generalizations of the conventional discrete transforms with an additional fractional parameter. They have been introduced in the 1990s. Noteworthy discrete fractional transforms, such as the discrete fractional Fourier transform (DFrFT) [1][2][3], discrete fractional Hartley transform [4], discrete fractional cosine and sine transforms [5], discrete fractional Hadamard transform (DFrHT) [6] have been defined. They are widely used in various fields of science and technology, including image representation and compression [7], image encryption [8], digital watermarking [9], adaptive filtering [10,11] and others. Among these transforms, the DFrFT has been applied in most practical contexts. To date, several effective algorithms have been developed to implement this transform [5,12,13]. In [14], the DFrFT was generalized, so that it depends on many fractional parameters. It was named the multiple-parameter discrete fractional Fourier transform (MPDFrFT). In our previous work [12], we have proposed a fast algorithm for the DFrFT and MPDFrFT. However, even this algorithm was quite computationally intensive. Therefore, our goal was to further reduce the complexity of calculations while preserving the properties of fractionality. In this article, we present a new discrete pseudo-fractional Fourier transform and its multiple-parameter version as well as a fast algorithm to calculate them. A similar approach was proposed in our previous work [15] but concerned with the use of DFrHT kernels. As a result of this approach, a new orthogonal base was obtained and a fast algorithm for discrete transform in this base was synthesized. This transform can be implemented using small-size discrete fractional Fourier transforms, for which algorithms were described in [13]. The algorithm of the proposed transform has a parallel, simple and modular structure, which makes it convenient for hardware implementation.

Mathematical Background
The normalized discrete Fourier transform (DFT) matrix of order N is defined as follows: where w N = e −j 2π N , j is the imaginary unit, and 1/ √ N is the normalizing factor. Since F N is a unitary and symmetric matrix, it can be diagonalized, i.e., represented as a product [3] where Λ N is an N × N diagonal matrix, whose diagonal entries are the eigenvalues of F N . Z N is the matrix whose columns are normalized mutually orthogonal eigenvectors N corresponds to the eigenvalue λ k . The matrix F N has only four different eigenvalues: 1, −1, j, −j and, for N ≥ 4, the eigenvalues are degenerated and the set of eigenvectors is not unique [3]. For this reason, a particular eigenvector set should be specified.
The discrete fractional Fourier transform (DFrFT) matrix is a power of the DFT matrix with one real parameter a in the exponent. This parameter is usually a fraction, hence the term fractional. The DFrFT matrix could be calculated from the eigenvalue decomposition of DFT matrix [3] F For a = 0, the DFrFT matrix is equal to the identity matrix, and, for a = 1, it becomes the ordinary DFT matrix. The definition (3) of DFrFT was first introduced by Pei and Yeh [1,2]. They defined the DFrFT in terms of a particular set of eigenvectors, namely the discrete counterpart of the set of Hermite-Gaussian functions. In our article, we assume that the set of eigenvectors of the DFT matrix and their ordering were determined according to the procedure described in [3]. These eigenvectors are not determined in an analytic form. They can be calculated using numerical methods.

Definition of a New Transform
To determine a new transform matrix, we firstly represent the order of the transform matrix N, which is assumed to be a composite number, as a product of positive integers N = N 1 , N 2 , . . . , N k , where k ≥ 2. The numbers N 1 , N 2 , . . . , N K may but need not to be prime numbers.
We define the fractional-like Fourier transform matrix in the following way: where F a N k , k = 1, 2, . . . , K, are the discrete fractional Fourier transform matrices with fractional parameter a, which are determined according to Equation (3), and ⊗ means the Kronecker product of matrices operation [16]. The definition (4) of the matrix F (a,N 1 ,N 2 ,...,N K ) depends on the way how N is decomposed into the product of integers. Now, we will show that the matrix defined in this way satisfies most of the requirements imposed on the fractional transform matrix, i.e., the unitarity and index additivity. However, it is not reducible to the usual Fourier transform matrix when the fractional parameter a is equal to one. At the beginning, we prove unitarity of the matrix F where * means conjugate transpose of a matrix.
For the Kronecker product of any invertible matrices A and B, the following equality is true [16]: In case of proving unitarity of the matrix F (a,N 1 ,N 2 ,...,N K ) N , the following property of the Kronecker product of any matrices A and B should be used [16]: Using this property and Equations (5) and (7), we may write is unitary. The next requirement for a fractional transform is index additivity, which means that for any real fractional parameters a and b, the result of two successive transforms with parameters a and b is the same as the result of one transform with parameter a+b. To demonstrate that the defined transform matrix (4) satisfies this property, we can use the appropriate property of the matrix F a N k again. This matrix, being a matrix of discrete fractional transform, fulfills We will use this fact and the next property of the Kronecker product, which is true for matrices A, B, C and D of such sizes that the matrix products AC and BD exist [16] We obtain F Moreover, the defined matrix (4) is also symmetric, i.e., Its symmetry is caused by symmetry of the matrix F a N k . Taking into account the above considerations, we can conclude that the defined matrix has many properties in common with the DFrFT matrix. However, since this matrix is not equal to the DFrFT matrix and does not fulfill all the requirements for fractional transforms, we call it the discrete pseudo-fractional Fourier transform matrix and the transform determined by it-the discrete pseudo-fractional Fourier transform (DPFrFT).

Visualization of the DPFrFT
The DFrFT was specified for the input data vector, however, it is often used in image processing. We chose this type of input data because it is suitable for transform visualization. The two-dimensional transform of an image consists in determining one-dimensional transforms of all columns and all rows of the matrix representing this image. In matrix notation, it can be written as follows: where X N×M is the matrix representing an input image, F N and F M -matrices of onedimensional transforms, which operate on columns and rows of the input matrix, respectively, and Y N×M is the output matrix.
Since we chose the DFrFT as the starting point, we will compare the results of DFrFT and DPFrFT transforms obtained for some images. The calculations and visualizations have been made in the Matlab R2020b environment. Two well-known testing images-"Lena" and "Boat", of size 512 × 512 pixels, were chosen to compare the DFrFT and DPFrFT transforms. In the case of DPFrFT, the number N = 512 was decomposed into the product N 1 N 2 N 3 = 16 × 8 × 4, so according to Equation (4) Before calculating each transform, the grayscale image, represented by the matrix X 512 , is converted from the type uint8 into the type double, with values in the range [0, 1], using the Matlab function im2double. Then, the transforms were performed on the rows and columns of the input matrix, according to Equations (13) and (14), respectively. Since the results are complex matrices, the modulus, the real part and the imaginary part of each output will be presented. In order to visualize the outputs (the modulus, the real part and the imaginary part), they should also be scaled to the range [0, 1]. Logarithmic scaling was performed as follows: where Y can be the modulus, the real part, or the imaginary part of the output matrix of each transform. Y scaled is the matrix Y scaled to the range [0, 1]. | · | denotes the absolute value and max is the largest value in the matrix. In the end, each scaled matrix is converted from the type double into the type uint8, using the Matlab function im2uint8, and visualized using the Matlab function imshow. Figures 1-8 show the obtained images for the four values of the fractional parameter a. In each figure, at the top, there is the original image.
In the middle, from left to right, there are the images of the modulus, the real part and the imaginary part of its DFrFT, respectively. Analogously, in the bottom, there are the modulus, the real part and the imaginary part of DPFrFT. It is easy to notice that an increase in the fractional parameter a, both in the case of DFrFT and DPFrFT, makes the obtained images less and less similar to the original image. In other words, they are becoming more and more difficult to recognize. For the same value of the fractional parameter a in DPFrFT images, the original image is more difficult to recognize than in DFrFT images, which may indicate that DPFrFT is more useful for hiding/encrypting images than DFrFT. Furthermore, it can be seen that in the DFrFT images, the large-scale elements of the original image are more visible, while in the case of DPFrFT, the details are more visible.

Multiple-Parameter Discrete Pseudo-Fractional Fourier Transform
There are applications of DFrFT, such as image encryption, in which it is more preferable when the transform has many fractional parameters [14]. Such a transform gives more degrees of freedom when the signal is represented in the MPDFrFT domain, which makes its correct decryption more difficult. It is also possible to generalize in this way the DPFrFT matrix (4). Let us define the matrix F (a 1 ,a 2 ,...,a K ,N 1 ,N 2 where a 1 , a 2 , . . . , a K are some real fractional parameters. The transform to which the matrix F The matrix F (a 1 ,a 2 ,...,a K ,N 1 ,N 2 ,...,N K ) N is invertible and its inverse matrix is equal to (F (a 1 ,a 2 ,...,a K ,N 1 ,N 2 ,...,N K ) N As in the case of the DPFrFT matrix, it can be easily shown that the MPDPFrFT matrix is unitary and satisfies index additivity. When all fractional parameters are equal to each other then the MPDPFrFT reduces to ordinary DPFrFT, that is, for a 1 = a 2 = . . . = a K = a.

Fast Algorithm for the Multiple-Parameter Discrete Pseudo-Fractional Fourier Transform
Only the algorithm for MPDPFrFT will be presented because an algorithm for the ordinary DPFrFT is its special case in which all fractional parameters are equal to each other.
We use the Fundamental Tensor-Product Factorization Theorem, which was proven in [17]. It states that the tensor product of any square matrices A n 1 , A n 2 , . . . , A n t , where A n i is an n i × n i matrix, can be factorized as follows: Moreover, the factorization is correct for any permutation of the factors ( ). In this equation, I n is the identity matrix of order n, N(i) = n 1 n 2 . . . n i , N(0) = 1 and N = N(t), where t denotes the number of factors.
Using the mentioned theorem to the tensor product in Equation (16) and noticing that in this case t = K and N = N(K), we obtain If we change the order of the factors from the last to the first, which we can do according to the mentioned theorem, and introduce a new index j = K − i + 1, then we have Since the left side of the above equality is equal to the matrix of MPDPFrFT, so when we enter the designation for k = 1, 2, . . . , K, we may write F (a 1 ,a 2 ,...,a K ,N 1 ,N 2 ,..., so, we have the following matrix-vector procedure to calculate the output of MPDPFrFT y (a 1 ,a 2 ,...,a K ,N 1 ,N 2 ,..., For example, we consider the case of N = 12 = 2 × 2 × 3. In this case K = 3, N 1 = 2, N 2 = 2 and N 3 = 3. In this case, the algorithm for computing the output vector is represented by the following matrix-vector procedure: y (a 1 ,a 2 ,a 3 ,2,2,3) 12 (30) Figure 9 shows a data flow diagram of a fast algorithm for a 12-point MPDPFrFT. The data flow diagram is oriented from left to right and straight lines represent data transfer operations. In turn, the rectangles denote small-sized DFrFT kernels that implement the multiplications of the matrices F a k N k , k = 1, 2, 3, by the corresponding input vectors.

Computational Complexity
We will now estimate the total number of multiplications and additions of complex numbers needed to compute the output of MPDPFrFT (or DPFrFT if a 1 = a 2 = . . . = a K = a), for the complex-valued input vector x N , according to Equation (25). We assume that the matrices R (a k ,N k ,k) N , k = 1, 2, . . . , K are factorized as in Equation (23) and the matrices F a k N k have been prepared in advance for the selected fractional parameters a k , k = 1, 2, . . . , K. The numbers of arithmetic operations depend on the method of factorization of the number N into the product N 1 N 2 . . . N K and on the algorithm used for the calculation of the DFrFT for individual k. When performing the matrix-vector multiplication according to Equation (25), from the right side to the left, we have K steps (corresponding to the multiplication of the matrices R (a k ,N k ,k) N , k = 1, 2, . . . , K by a vector). At the step k, there are N/N k multiplications of the matrix F a k N k by the corresponding subvectors needed. If we carry out the multiplication of the latter matrix by the corresponding subvector in the classical way, we have N 2 k multiplications and N k (N k − 1) additions of complex numbers. Hence, the total numbers of multiplications L × and additions L + of complex numbers are equal to For the example shown in Figure 1, i.e., N = 12 = 2 × 2 × 3, the numbers of additions and multiplications of complex numbers are equal to 84 and 48, respectively.
It is necessary to emphasize that the matrices F a k N k , k = 1, 2, . . . , K have a special structure; therefore, the numbers of arithmetic operations when calculating the matrixvector product of each of them by the vector can be further reduced by about a half [12,13].

Conclusions
The paper describes a new type of discrete fractional-like transform called the discrete pseudo-fractional Fourier transform. We also present a generalization of this transform, called multiple-parameter discrete pseudo-fractional Fourier transform. Finally, we propose a fast algorithm for calculating this transform. A distinctive feature of this transform and its algorithm is that while maintaining all the advantages of fractionality and multiparametricity, a significant reduction in computational complexity is achieved compared to the true DFrFT [12]. In addition, the algorithm has a regular modular structure. Typical cores of the algorithm are standard small-size DFrFTs, which allows them to be used as unified building blocks in software or hardware implementation [13]. This will simplify and unify the process of developing programs and digital signal processors that implement the considered algorithm. Matrices of the proposed transform are defined by the smaller size DFrFT matrices. Taking into account the specific properties of the introduced transform, one can predict its interesting applications. However, these questions are beyond the scope of our article. Funding: This research received no external funding.

Data Availability Statement:
All data included in this study are available upon request by contacting with the corresponding author.