One-Dimensional Quaternion Discrete Fourier Transform and an Approach to Its Fast Computation

Abstract

This paper proposes a new method for calculating the quaternion discrete Fourier transform for one-dimensional data. Although the computational complexity of the proposed method still belongs to the $O\left(N{log}_{2}N\right)$ class, it allows us to reduce the total number of arithmetic operations required to perform it compared to other known methods for computing this transform. Moreover, compared to the method using symplectic decomposition, the presented method does not require changing the basis in the subspace of pure quaternions and, consequently, calculating the new basis vectors and change-of-basis matrix.

1. Introduction

The discrete Fourier transform (DFT) is a fundamental tool for data processing and analysis in many fields of science and technology [1]. The DFT is especially used in digital signal and image processing. If a one-dimensional discrete Fourier transform is used in the processing and analysis of one-dimensional signals, such as speech and other audio signals [2,3], then in the case of image processing and recognition, we are dealing with a two-dimensional DFT [4,5,6,7,8,9].

The rapid development of signal and image processing technologies has led to advanced data processing techniques using data representation and processing in the hypercomplex domain. This led to the proposal to represent each pixel of a color image as a quaternion [10,11,12,13,14].

The old approach to using color channels to process color images was to use separate channels for the primary colors red, green, and blue. This approach divided the color image into three separate channels. In fact, three 2D matrices were formed, as a result of the presence of three channels for color images, and the computations with each matrix were carried out independently in grayscale. With this approach, due to channel separation, the cross-correlation between channels was lost at the start of the process because the color image was not considered as a whole [14,15].

In the quaternion-based method, a color image was considered as one matrix consisting of pure quaternion numbers. Such a matrix is processed as a whole. In order to convert an RGB image into a quaternion matrix, the RGB channels were simply inserted into the vector part of the quaternion matrix. In addition, such a representation of the pixel matrix made it possible to significantly reduce the computational complexity of color image processing methods compared to three-channel processing.

A natural development of research within the framework of the quaternion-valued representation of pixels in color images was the definition of a two-dimensional quaternion discrete Fourier transform (2-D QDFT) and the development of fast algorithms for 2-D QDFT implementation [12,16,17,18,19].

As for the one-dimensional QDFT, this issue has not been studied in detail, since many computer scientists simply did not see the use of such a transform. That is why there are relatively few research studies on this topic [16,19,20,21,22,23,24]. In this article, we would like to fill this gap and give a detailed definition of 1-D QDFT, as well as propose a fast method for its implementation.

2. Quaternions

Quaternions were introduced by Hamilton in 1843 [25]. They fulfill the usual rules of algebra, but their multiplication is not commutative. Quaternions have four components, one real and three imaginary. In usual notation, which comes from complex numbers, the quaternion q is written in the form

$$q=a+bi+cj+dk$$

(1)

where a, b, c, and d are real numbers, and i, j, k are imaginary units that obey the following rules of multiplication:

$$i^2=j^2=k^2=-1$$

(2)

$$ij=k,ji=-k,jk=i,kj=-i,ki=j,ik=-j$$

(3)

Quaternions are also called hypercomplex numbers (together with coquaternions, biquaternions, tessarines, octonions, sedenions, etc.). A quaternion has a real part—a in (1)—and an imaginary part, which has three components and is usually denoted by $V\left(q\right)=bi+cj+dk$. The imaginary part $V\left(q\right)$ of a quaternion may be associated with a 3-space vector ${[b,c,d]}^T$, and it is often called the vector part of the quaternion. For this reason, the real part is often called the scalar part of the quaternion q, and denoted by $S\left(q\right)$. The whole quaternion may be represented as $q=S\left(q\right)+V\left(q\right)$, i.e., the sum of its scalar and vector parts. A quaternion in which the real (or scalar) part is equal to zero is called a pure quaternion.

The conjugate of a quaternion q is designed by $\overline{q}$ and obtained, as in complex numbers, by negating the imaginary part, so

$$\overline{q}=a-bi-cj-dk$$

(4)

The modulus of a quaternion q follows the definition for complex numbers, so

$$\left|q\right|=\sqrt{q\overline{q}}=\sqrt{a^2+b^2+c^2+d^2}$$

(5)

A quaternion with a unit modulus is called a unit quaternion. The inverse of a quaternion q, denoted by $q^{-1}$, such that $q{q}^{-1}=q^{-1}q=1$ is given by

$$q^{-1}=\frac{\overline{q}}{{\left|q\right|}^2}$$

(6)

Euler’s formula for the complex exponential is also true for quaternions, so

$$e^{\mu \varphi }=cos\varphi +\mu sin\varphi$$

(7)

where $\mu$ is a unit pure quaternion and $\varphi$ is a real number–angle. Any quaternion q may be represented in the polar form as

$$q=\left|q\right|e^{\mu \varphi }$$

(8)

where $\mu$ is called the eigenaxis, and $\varphi$ is called the eigenangle of the quaternion q.

A quaternion q may be also represented as a generalized complex number [26]

$$q=w+rj$$

(9)

where its real and imaginary parts are ordinary complex numbers

$$w=a+bir=c+di$$

(10)

The imaginary units i and j should be different and orthogonal (where the orthogonality of pure quaternions is understood as for vectors in ${\mathbb{R}}^3$ vector space). Definition (9) of a quaternion is known as the Cayley–Dickson form. Using this form, the quaternion q may be written as

$$q=(a+bi)+(c+di)j$$

(11)

where it can be seen by multiplying out, using the rules (3), that the result is the same as in the definition (1). The two parts $w=a+bi$ and $r=c+dj$ are called the simplex and perplex parts of the quaternion q, respectively, and the decomposition of the quaternion into these parts is known as symplectic decomposition [26]. Such a decomposition of a quaternion can be used to calculate the QDFT, as will be shown in Section 4.

The product of two quaternions $q=a+bi+cj+dk=S\left(q\right)+V\left(q\right)$ and $\widehat{q}=\widehat{a}+\widehat{b}i+\widehat{c}j+\widehat{d}k=S\left(\widehat{q}\right)+V\left(\widehat{q}\right)$ can be represented by the products of their scalar and vector parts

$$q\widehat{q}=a\widehat{a}-(b\widehat{b}+c\widehat{c}+d\widehat{d})+a(\widehat{b}i+\widehat{c}j+\widehat{d}k)+$$

$$\widehat{a}(bi+cj+dk)+(c\widehat{d}-d\widehat{c})i+(d\widehat{b}-b\widehat{d})j+(b\widehat{c}-c\widehat{b})k=$$

$$S\left(q\right)S\left(\widehat{q}\right)-V\left(q\right)·V\left(\widehat{q}\right)+S\left(q\right)V\left(\widehat{q}\right)+S\left(\widehat{q}\right)V\left(q\right)+V\left(q\right)\times V\left(\widehat{q}\right)$$

(12)

where “·” denotes the scalar product and “×” denotes the cross product of vectors. Generally, quaternion multiplication is neither commutative nor anti-commutative because the dot product is commutative and the cross product is anti-commutative. When q and $\widehat{q}$ are pure quaternions, then $S\left(q\right)=S\left(\widehat{q}\right)=0$, so $q\widehat{q}=V\left(q\right)·V\left(\widehat{q}\right)+V\left(q\right)\times V\left(\widehat{q}\right)$. If, moreover, $V\left(q\right)·V\left(\widehat{q}\right)=0$, we say that pure quaternions q and $\widehat{q}$ are perpendicular ($q\perp \widehat{q}$). Then $q\widehat{q}=V\left(q\right)\times V\left(\widehat{q}\right)$ and $q\widehat{q}=-\widehat{q}q$. In general, two quaternions whose vector parts are parallel, in the usual vector sense, are called parallel quaternions. Likewise, if their vector parts are perpendicular, they are called perpendicular quaternions [11].

3. Two Types of Quaternion Discrete Fourier Transform

Although the quaternion Fourier transform is mainly used for two-dimensional data, i.e., images, it has been specified for one-dimensional data as well [19,21]. We will only deal with the discrete version of the one-dimensional quaternion Fourier transform.

Let $x^{\left(r\right)}\left(n\right)$, $x^{\left(i\right)}\left(n\right)$, $x^{\left(j\right)}\left(n\right)$, and $x^{\left(k\right)}\left(n\right)$ be the components of a discrete quaternion signal $x\left(n\right)$ with the number of samples N. Then

$$x\left(n\right)=x^{\left(r\right)}\left(n\right)+x^{\left(i\right)}\left(n\right)i+x^{\left(j\right)}\left(n\right)j+x^{\left(k\right)}\left(n\right)k$$

(13)

for $n=0,\dots ,N-1$, where $x^{\left(r\right)}\left(n\right)$, $x^{\left(i\right)}\left(n\right)$, $x^{\left(j\right)}\left(n\right)$, and $x^{\left(k\right)}\left(n\right)$ are real numbers.

Let $\mu$ be a unit pure quaternion and ${\mu }_i$, ${\mu }_j$, ${\mu }_k$ be its coefficients associated with the imaginary units i, j, k, respectively:

$$\mu ={\mu }_{i}i+{\mu }_{j}j+{\mu }_{k}k$$

(14)

where

$${\left|\mu \right|}^2=\mu \overline{\mu }={\left({\mu }_i\right)}^2+{\left({\mu }_j\right)}^2+{\left({\mu }_k\right)}^2=1$$

(15)

The choice of $\mu$ is arbitrary, but it matters. In applications, the case ${\mu }_i={\mu }_j={\mu }_k$ is often used; thus, taking into account (15), we obtain $\mu =(i+j+k)/\sqrt{3}$.

Since quaternion multiplication is not commutative, two different 1-D quaternion discrete Fourier transforms were defined, which are referred to as right-side and left-side QDFTs.

The right-side 1-D quaternion discrete Fourier transform is defined as follows:

$$y\left(m\right)=\sum _{n=0}^{N-1}x\left(n\right)e^{-\frac{\mu 2\pi mn}{N}}$$

(16)

for $m=0,\dots ,N-1$, where the quaternion exponents are to the right of $x\left(n\right)$, and $\mu$ is any unit pure quaternion.

It should be noted that for any choice of $\mu$, ${\mu }^2=-1$, so this transform can be regarded as a generalization of the standard complex discrete Fourier transform in which the imaginary unit i has been generalized to a vector imaginary unit $\mu$. This means that if $\mu =i$ and the input signal is complex-valued, that is all samples have values $(x^{\left(r\right)}\left(n\right)+x^{\left(i\right)}i)$, then this transform is reduced to the usual complex-valued discrete Fourier transform.

According to (7), the quaternion exponent can be written in the form

$$e^{-\frac{\mu 2\pi mn}{N}}=cos\frac{2\pi mn}{N}-\mu sin\frac{2\pi mn}{N}$$

(17)

or, if we introduce the denotations

$$a_{m,n}=cos\frac{2\pi mn}{N},b_{m,n}=sin\frac{2\pi mn}{N}$$

(18)

and take into account (14), the quaternion exponent can by written as

$$e^{-\frac{\mu 2\pi mn}{N}}=a_{m,n}-({\mu }_{i}i+{\mu }_{j}j+{\mu }_{k}k)b_{m,n}$$

(19)

Putting expressions (13) and (19) into (16), the right-side quaternion discrete Fourier transform can be written in the form

$$y\left(m\right)=\sum _{n=0}^{N-1}[x^{\left(r\right)}\left(n\right)+x^{\left(i\right)}\left(n\right)i+x^{\left(j\right)}\left(n\right)j+x^{\left(k\right)}\left(n\right)][a_{m,n}-({\mu }_{i}i+{\mu }_{j}j+{\mu }_{k}k)b_{m,n}]$$

(20)

After multiplying the expressions in square brackets, according to the rules (2) and (3) of multiplication of imaginary units, we can write the result as

$$y\left(m\right)=y^{\left(r\right)}\left(m\right)+y^{\left(i\right)}\left(m\right)i+y^{\left(j\right)}\left(m\right)j+y^{\left(k\right)}\left(m\right)k$$

(21)

where $y^{\left(r\right)}\left(m\right)$, $y^{\left(i\right)}\left(m\right)$, $y^{\left(j\right)}\left(m\right)$, and $y^{\left(k\right)}\left(m\right)$ are real numbers. It is easy to check that they take the following forms:

$$y^{\left(r\right)}\left(m\right)=\sum _{n=0}^{N-1}\left\{x^{\left(r\right)}\left(n\right)a_{m,n}+{\mu }_i{x}^{\left(i\right)}\left(n\right)b_{m,n}+{\mu }_j{x}^{\left(j\right)}\left(n\right)b_{m,n}+{\mu }_k{x}^{\left(k\right)}\left(n\right)b_{m,n}\right\}$$

(22)

$$y^{\left(i\right)}\left(m\right)=\sum _{n=0}^{N-1}\left\{x^{\left(i\right)}\left(n\right)a_{m,n}-{\mu }_i{x}^{\left(r\right)}\left(n\right)b_{m,n}+{\mu }_j{x}^{\left(k\right)}\left(n\right)b_{m,n}-{\mu }_k{x}^{\left(j\right)}\left(n\right)b_{m,n}\right\}$$

(23)

$$y^{\left(j\right)}\left(m\right)=\sum _{n=0}^{N-1}\left\{x^{\left(j\right)}\left(n\right)a_{m,n}-{\mu }_i{x}^{\left(k\right)}\left(n\right)b_{m,n}-{\mu }_j{x}^{\left(r\right)}\left(n\right)b_{m,n}+{\mu }_k{x}^{\left(i\right)}\left(n\right)b_{m,n}\right\}$$

(24)

$$y^{\left(k\right)}\left(m\right)=\sum _{n=0}^{N-1}\left\{x^{\left(k\right)}\left(n\right)a_{m,n}+{\mu }_i{x}^{\left(j\right)}\left(n\right)b_{m,n}-{\mu }_j{x}^{\left(i\right)}\left(n\right)b_{m,n}-{\mu }_k{x}^{\left(r\right)}\left(n\right)b_{m,n}\right\}$$

(25)

The left-side 1-D quaternion discrete Fourier transform is defined as

$$z\left(m\right)=\sum _{n=0}^{N-1}{e}^{-\frac{\mu 2\pi mn}{N}}x\left(n\right)$$

(26)

for $m=0,\dots ,N-1$, where the quaternion exponents are on the left side in multiplications. Taking into account (13) and (19), the right-side quaternion discrete Fourier transform (26) can be written in the form

$$z\left(m\right)=\sum _{n=0}^{N-1}[a_{m,n}-({\mu }_{i}i+{\mu }_{j}j+{\mu }_{k}k)b_{m,n}][x^{\left(r\right)}\left(n\right)+x^{\left(i\right)}\left(n\right)i+x^{\left(j\right)}\left(n\right)j+x^{\left(k\right)}\left(n\right)]$$

(27)

After multiplying the expressions in square brackets, according to the rules (2) and (3) of multiplication of imaginary units, we can write the result as

$$z\left(m\right)=z^{\left(r\right)}\left(m\right)+z^{\left(i\right)}\left(m\right)i+z^{\left(j\right)}\left(m\right)j+z^{\left(k\right)}\left(m\right)k$$

(28)

where $z^{\left(r\right)}\left(m\right)$, $z^{\left(i\right)}\left(m\right)$, $z^{\left(j\right)}\left(m\right)$, and $z^{\left(k\right)}\left(m\right)$ are real numbers. It is easy to check that they take the following forms:

$$z^{\left(r\right)}\left(m\right)=\sum _{n=0}^{N-1}\left\{x^{\left(r\right)}\left(n\right)a_{m,n}+{\mu }_i{x}^{\left(i\right)}\left(n\right)b_{m,n}+{\mu }_j{x}^{\left(j\right)}\left(n\right)b_{m,n}+{\mu }_k{x}^{\left(k\right)}\left(n\right)b_{m,n}\right\}$$

(29)

$$z^{\left(i\right)}\left(m\right)=\sum _{n=0}^{N-1}\left\{x^{\left(i\right)}\left(n\right)a_{m,n}-{\mu }_i{x}^{\left(r\right)}\left(n\right)b_{m,n}-{\mu }_j{x}^{\left(k\right)}\left(n\right)b_{m,n}+{\mu }_k{x}^{\left(j\right)}\left(n\right)b_{m,n}\right\}$$

(30)

$$z^{\left(j\right)}\left(m\right)=\sum _{n=0}^{N-1}\left\{x^{\left(j\right)}\left(n\right)a_{m,n}+{\mu }_i{x}^{\left(k\right)}\left(n\right)b_{m,n}-{\mu }_j{x}^{\left(r\right)}\left(n\right)b_{m,n}-{\mu }_k{x}^{\left(i\right)}\left(n\right)b_{m,n}\right\}$$

(31)

$$z^{\left(k\right)}\left(m\right)=\sum _{n=0}^{N-1}\left\{x^{\left(k\right)}\left(n\right)a_{m,n}-{\mu }_i{x}^{\left(j\right)}\left(n\right)b_{m,n}+{\mu }_j{x}^{\left(i\right)}\left(n\right)b_{m,n}-{\mu }_k{x}^{\left(r\right)}\left(n\right)b_{m,n}\right\}$$

(32)

Unlike the complex-valued DFT, the left-side and right-side QDFTs are slightly different transforms. From (22) and (29), we can see that the real parts of both transforms are the same but the imaginary parts differ from each other in the signs of some components of the sums (compare (23) with (30), (24) with (31), and (25) with (32)). Obviously, this follows from the non-commutative nature of quaternions multiplication. However, if we choose $\mu =i$ (i.e., ${\mu }_i=1$ and ${\mu }_j={\mu }_k=0$) and the complex-valued input signal $x\left(n\right)$ (i.e., $x^{\left(j\right)}\left(n\right)=x^{\left(k\right)}\left(n\right)=0$), then both of these transforms will be reduced to the usual complex-valued DFT.

Moreover, by analyzing (22)–(25) and (29)–(32), it is easy to see that all components of both QDFTs are linear combinations of the real and imaginary parts of four ordinary Fourier transforms for the real signals $x^{\left(r\right)}\left(n\right)$, $x^{\left(i\right)}\left(n\right)$, $x^{\left(j\right)}\left(n\right)$, and $x^{\left(k\right)}\left(n\right)$—components of the quaternion signal $x\left(n\right)$. This will be further explained in Section 5. Therefore, for determining the right-side or left-side QDFT, the calculation of four discrete Fourier transforms for the four real-valued components of the quaternion signal is needed.

Figure 1 shows graphical representations of four components of the right-side and the left-side QDFTs for the quaternion signal $x\left(n\right)=x^{\left(r\right)}\left(n\right)+x^{\left(i\right)}\left(n\right)i+x^{\left(j\right)}\left(n\right)j+x^{\left(k\right)}\left(n\right)k$ with $N=16$ samples, where

$$x^{\left(r\right)}\left(n\right)=sin\frac{n\pi }{N}-n$$

(33)

$$x^{\left(i\right)}\left(n\right)=-4+n$$

(34)

$$x^{\left(j\right)}\left(n\right)=10+sin\frac{2n\pi }{N}$$

(35)

$$x^{\left(k\right)}\left(n\right)=5+2cos\frac{n\pi }{N}$$

(36)

4. Calculation of QDFTs by Symplectic Decomposition

The use of the symplectic decomposition method allows both QDFTs to be computed by performing only two complex-valued discrete Fourier transforms instead of four. This approach was used in [26] to calculate the two-dimensional left-side QDFT. In this work, we will show how to also use the symplectic decomposition method to reduce the computational complexity of calculating one-dimensional right-side and left-side QDFTs.

If we take arbitrary two unit pure quaternions ${\mu }_1$ and ${\mu }_2$, such that ${\mu }_1\perp {\mu }_2$ (pure quaternions are perpendicular when their dot product is equal to zero), and the third unit pure quaternion ${\mu }_3={\mu }_1{\mu }_2$, we can represent the quaternion signal

$$x\left(n\right)=x^{\left(r\right)}\left(n\right)+x^{\left(i\right)}\left(n\right)i+x^{\left(j\right)}\left(n\right)j+x^{\left(k\right)}\left(n\right)k$$

(37)

as

$$x\left(n\right)=x^{\left(r\right)}\left(n\right)+x^{\left({\mu }_1\right)}\left(n\right){\mu }_1+x^{\left({\mu }_2\right)}\left(n\right){\mu }_2+x^{\left({\mu }_3\right)}\left(n\right){\mu }_3$$

(38)

The change of the imaginary units from i, j, k to ${\mu }_1$, ${\mu }_2$, ${\mu }_3$ corresponds to the change of the basis in ${\mathbb{R}}^3$ vector space from ${[1,0,0]}^T$, ${[0,1,0]}^T$, ${[0,0,1]}^T$ to another basis of three mutually orthogonal unit vectors; so for each sample, we obtain the new coordinates by multiplying the matrix inverse to the change-of-basis matrix by the old coordinates vector

$$\left[\begin{array}{c}{x}^{\left({\mu }_1\right)}\left(n\right)\\ x^{\left({\mu }_2\right)}\left(n\right)\\ x^{\left({\mu }_3\right)}\left(n\right)\end{array}\right]={\mathbf{S}}_3^{-1}\left[\begin{array}{c}{x}^{\left(i\right)}\left(n\right)\\ x^{\left(j\right)}\left(n\right)\\ x^{\left(k\right)}\left(n\right)\end{array}\right]$$

(39)

The change-of-basis matrix ${\mathbf{S}}_3$ depends on the choice of new imaginary units. For example, if we take the imaginary units ${\mu }_1$, ${\mu }_2$, ${\mu }_3$, as in [26]

$${\mu }_1=\frac{1}{\sqrt{3}}i+\frac{1}{\sqrt{3}}j+\frac{1}{\sqrt{3}}k$$

(40)

$${\mu }_2=\frac{1}{\sqrt{2}}j-\frac{1}{\sqrt{2}}k$$

(41)

$${\mu }_3=-\frac{2}{\sqrt{6}}i+\frac{1}{\sqrt{6}}j+\frac{1}{\sqrt{6}}k$$

(42)

then the change-of-basis matrix has the form

$${\mathbf{S}}_3=\left[\begin{array}{ccc}\frac{1}{\sqrt{3}}& 0& -\frac{2}{\sqrt{6}}\\ \frac{1}{\sqrt{3}}& \frac{1}{\sqrt{2}}& \frac{1}{\sqrt{6}}\\ \frac{1}{\sqrt{3}}& -\frac{1}{\sqrt{2}}& \frac{2}{\sqrt{6}}\end{array}\right]$$

(43)

In Cayley–Dickson form (38) has the form

$$x\left(n\right)=[x^{\left(r\right)}\left(n\right)+x^{\left({\mu }_1\right)}\left(n\right){\mu }_1]+[x^{\left({\mu }_2\right)}\left(n\right)+x^{\left({\mu }_3\right)}\left(n\right){\mu }_1]{\mu }_2$$

(44)

This is a symplectic decomposition of each sample $x\left(n\right)$ into a simplex part

$$x^{\left(s\right)}\left(n\right)=x^{\left(r\right)}\left(n\right)+x^{\left({\mu }_1\right)}\left(n\right){\mu }_1$$

(45)

and a perplex part

$$x^{\left(p\right)}\left(n\right)=x^{\left({\mu }_2\right)}\left(n\right)+x^{\left({\mu }_3\right)}\left(n\right){\mu }_1$$

(46)

both isomorphic to the standard complex number. The result is a pair of complex signals. They may be stored using the same space as the quaternion signal $x\left(n\right)$. The symplectic decomposition splits a quaternion into two perpendicular planes. These planes intersect only at the origin (in 4-space). Each of these planes may be regarded as an Argand plane. The first is the Argand plane of the simplex part; its real axis is identical to the scalar axis of quaternion space and its imaginary axis is ${\mu }_1$. The second is the Argand plane of the perplex part; this plane is perpendicular to both axes of the Argand plane of the simplex part. The real axis of this plane is identical to ${\mu }_2$ and its imaginary axis is ${\mu }_3$.

Let the axis of transform kernel of the right-side QDFT in (16) be chosen as ${\mu }_1$, i.e., $\mu ={\mu }_1$. Then (16), can be rewritten as

$$y\left(m\right)=\sum _{n=0}^{N-1}[x^{\left(s\right)}\left(n\right)+x^{\left(p\right)}\left(n\right){\mu }_2]e^{-\frac{{\mu }_{1}2\pi mn}{N}}$$

(47)

The sum in (47) can be written as the sum of two terms. The first term, which is the simplex part of $y\left(m\right)$, has the form

$$y^{\left(s\right)}\left(m\right)=\sum _{n=0}^{N-1}{x}^{\left(s\right)}\left(n\right)e^{-\frac{{\mu }_{1}2\pi mn}{N}}=\sum _{n=0}^{N-1}[x^{\left(r\right)}\left(n\right)+x^{\left({\mu }_1\right)}\left(n\right){\mu }_1]e^{-\frac{{\mu }_{1}2\pi mn}{N}}$$

(48)

This is isomorphic to the discrete Fourier transform for a complex input. The only difference is that the imaginary unit is called ${\mu }_1$ instead of i. Hence, calculation of $y^{\left(s\right)}\left(m\right)=y^{\left(r\right)}\left(n\right)+y^{\left({\mu }_1\right)}\left(n\right){\mu }_1$ comes down to the computation of a complex DFT for the simplex part $x^{\left(s\right)}\left(n\right)$ of $x\left(n\right)$. This can be performed using the FFT algorithm. The second term

$$y^{\left(p\right)}\left(m\right){\mu }_2=\sum _{n=0}^{N-1}{x}^{\left(p\right)}\left(n\right){\mu }_2{e}^{-\frac{{\mu }_{1}2\pi mn}{N}}=\sum _{n=0}^{N-1}{x}^{\left(p\right)}\left(n\right){\mu }_2\left(cos\frac{2\pi mn}{N}-{\mu }_{1}sin\frac{2\pi mn}{N}\right)=$$

$$\sum _{n=0}^{N-1}{x}^{\left(p\right)}\left(n\right)\left(cos\frac{2\pi mn}{N}+{\mu }_{1}sin\frac{2\pi mn}{N}\right){\mu }_2=\left(\sum _{n=0}^{N-1}[x^{\left({\mu }_2\right)}\left(n\right)+x^{\left({\mu }_3\right)}\left(n\right){\mu }_1]e^{\frac{{\mu }_{1}2\pi mn}{N}}\right){\mu }_2$$

(49)

where, during the transformation, we used the equality ${\mu }_2{\mu }_1=-{\mu }_1{\mu }_2$, which is true for pure, perpendicular quaternions. It is easy to see that $y^{\left(p\right)}\left(n\right)=y^{\left({\mu }_2\right)}\left(n\right)+y^{\left({\mu }_3\right)}\left(n\right){\mu }_1$ is the perplex part of $y\left(n\right)$. So, the calculation of $y^{\left(p\right)}\left(m\right)$ comes down to the computation of a complex inverse discrete Fourier transform (IDFT) for the perplex part $x^{\left(p\right)}\left(n\right)$ of $x\left(n\right)$, but without multiplying the result by $1/N$. This can be performed using the IFFT algorithm.

After calculating the simplex and perplex parts of the output signal $y\left(m\right)=y^{\left(r\right)}\left(n\right)+y^{\left({\mu }_1\right)}\left(n\right){\mu }_1+y^{\left({\mu }_2\right)}\left(n\right){\mu }_3+y^{\left({\mu }_3\right)}\left(n\right){\mu }_3$, we want to write $y\left(m\right)$ using the imaginary units i, j, and k instead of ${\mu }_1$, ${\mu }_2$, and ${\mu }_3$, i.e.,

$$y\left(m\right)=y^{\left(r\right)}\left(m\right)+y^{\left(i\right)}\left(m\right)i+y^{\left(j\right)}\left(m\right)j+y^{\left(k\right)}\left(m\right)k$$

(50)

In the associated ${\mathbb{R}}^3$ vector space, this corresponds to returning to the old base ${[1,0,0]}^T$, ${[0,1,0]}^T$, ${[0,0,1]}^T$, so it can be done by the operation inverse to (39), i.e.,

$$\left[\begin{array}{c}{y}^{\left(i\right)}\left(m\right)\\ y^{\left(j\right)}\left(m\right)\\ y^{\left(k\right)}\left(m\right)\end{array}\right]={\mathbf{S}}_3\left[\begin{array}{c}{y}^{\left({\mu }_1\right)}\left(m\right)\\ y^{\left({\mu }_2\right)}\left(m\right)\\ y^{\left({\mu }_3\right)}\left(m\right)\end{array}\right]$$

(51)

Summarizing the above considerations, in order to calculate the right-side QDFT using symplectic decomposition, the input signal should be represented by new imaginary units such that the first becomes the axis of the transform kernel. This corresponds to a change of the basis in the associated ${\mathbb{R}}^3$ vector space (for each sample of the quaternion-valued input signal). Then, a single complex DFT and a single complex IDFT should be calculated for the simplex and perplex parts of the quaternion-valued input signal, respectively. Finally, we have to return to the standard imaginary units i, j, and k, which corresponds to a return to the old basis in the associated ${\mathbb{R}}^3$ vector space.

A similar, and even simpler, consideration may be given to the left-side QDFT. Assuming that the axis of the transform kernel of the left-side QDFT in (26) is chosen as ${\mu }_1$, and using the symplectic decomposition of the quaternion input signal $x\left(n\right)$, (26), can be rewritten as

$$z\left(m\right)=\sum _{n=0}^{N-1}{e}^{-\frac{{\mu }_{1}2\pi mn}{N}}[x^{\left(s\right)}\left(n\right)+x^{\left(p\right)}\left(n\right){\mu }_2]$$

(52)

The sum in (52) can be written as the sum of two terms. The first term, which is the simplex part of $z\left(m\right)$ has the form

$$z^{\left(s\right)}\left(m\right)=\sum _{n=0}^{N-1}{e}^{-\frac{{\mu }_{1}2\pi mn}{N}}{x}^{\left(s\right)}\left(n\right)=\sum _{n=0}^{N-1}{e}^{-\frac{{\mu }_{1}2\pi mn}{N}}[x^{\left(r\right)}\left(n\right)+x^{\left({\mu }_1\right)}\left(n\right){\mu }_1]=$$

$$\sum _{n=0}^{N-1}[x^{\left(r\right)}\left(n\right)+x^{\left({\mu }_1\right)}\left(n\right){\mu }_1]e^{-\frac{{\mu }_{1}2\pi mn}{N}}$$

(53)

Calculation of $z^{\left(s\right)}\left(m\right)$ comes down to the computation of complex-valued DFTs for the simplex part $x^{\left(s\right)}\left(n\right)$ of $x\left(n\right)$. This can be performed using the FFT algorithm. In addition, $z^{\left(s\right)}\left(n\right)=y^{\left(s\right)}\left(n\right)$.

The second term has the form

$$z^{\left(p\right)}\left(m\right){\mu }_2=\sum _{n=0}^{N-1}{e}^{-\frac{{\mu }_{1}2\pi mn}{N}}{x}^{\left(p\right)}\left(n\right){\mu }_2=\left(\sum _{n=0}^{N-1}[x^{\left({\mu }_2\right)}\left(n\right)+x^{\left({\mu }_3\right)}\left(n\right){\mu }_1]e^{-\frac{{\mu }_{1}2\pi mn}{N}}\right){\mu }_2$$

(54)

It is easy to see that $z^{\left(p\right)}\left(n\right)$ is the perplex part of $z\left(n\right)$. Calculation of $z^{\left(p\right)}\left(m\right)$ also comes down to the computation of complex-valued DFTs for the perplex part $x^{\left(p\right)}\left(n\right)$ of $x\left(n\right)$. This can be performed using the FFT algorithm. Summarizing the above considerations, it can be argued that in order to calculate the left-side QDFT using a symplectic decomposition, the input signal should be represented by new imaginary units, the first of which will be equal to the axis of the transform kernel. This step is the same as in calculating the right-side QDFT. Then two complex DFTs should be calculated—for the simplex and perplex parts of the quaternion-valued input signal. In the end, as with the right-side QDFT calculation, the return to the standard imaginary units i, j, and k is needed.

5. The Simplest Method of QDFT Calculation

We will use the matrix notation to describe the QDFTs. Let ${\mathbf{x}}_N={[x_0,x_1,\dots ,x_{N-1}]}^T$ be a vector of samples of the quaternion input signal, where $x_n=x\left(n\right)$ for $n=0,1,\dots ,N-1$, so it can be written in the form

$${\mathbf{x}}_N=\left[\begin{array}{c}{x}_0^{\left(r\right)}\\ x_1^{\left(r\right)}\\ ⋮\\ x_{N-1}^{\left(r\right)}\end{array}\right]+\left[\begin{array}{c}{x}_0^{\left(i\right)}\\ x_1^{\left(i\right)}\\ ⋮\\ x_{N-1}^{\left(i\right)}\end{array}\right]i+\left[\begin{array}{c}{x}_0^{\left(j\right)}\\ x_1^{\left(j\right)}\\ ⋮\\ x_{N-1}^{\left(j\right)}\end{array}\right]j+\left[\begin{array}{c}{x}_0^{\left(k\right)}\\ x_1^{\left(k\right)}\\ ⋮\\ x_{N-1}^{\left(k\right)}\end{array}\right]k$$

(55)

Similarly ${\mathbf{y}}_N={[y_0,y_1,\dots ,y_{N-1}]}^T$ and ${\mathbf{z}}_N={[z_0,z_1,\dots ,z_{N-1}]}^T$ are the quaternion output vectors of the right-side and left-side QDFTs, respectively, where $y_n=y\left(n\right)$ and $z_n=z\left(n\right)$, so they can be written as

$${\mathbf{y}}_N=\left[\begin{array}{c}{y}_0^{\left(r\right)}\\ y_1^{\left(r\right)}\\ ⋮\\ y_{N-1}^{\left(r\right)}\end{array}\right]+\left[\begin{array}{c}{y}_0^{\left(i\right)}\\ y_1^{\left(i\right)}\\ ⋮\\ y_{N-1}^{\left(i\right)}\end{array}\right]i+\left[\begin{array}{c}{y}_0^{\left(j\right)}\\ y_1^{\left(j\right)}\\ ⋮\\ y_{N-1}^{\left(j\right)}\end{array}\right]j+\left[\begin{array}{c}{y}_0^{\left(k\right)}\\ y_1^{\left(k\right)}\\ ⋮\\ y_{N-1}^{\left(k\right)}\end{array}\right]k$$

(56)

$${\mathbf{z}}_N=\left[\begin{array}{c}{z}_0^{\left(r\right)}\\ z_1^{\left(r\right)}\\ ⋮\\ z_{N-1}^{\left(r\right)}\end{array}\right]+\left[\begin{array}{c}{z}_0^{\left(i\right)}\\ z_1^{\left(i\right)}\\ ⋮\\ z_{N-1}^{\left(i\right)}\end{array}\right]i+\left[\begin{array}{c}{z}_0^{\left(j\right)}\\ z_1^{\left(j\right)}\\ ⋮\\ z_{N-1}^{\left(j\right)}\end{array}\right]j+\left[\begin{array}{c}{z}_0^{\left(k\right)}\\ z_1^{\left(k\right)}\\ ⋮\\ z_{N-1}^{\left(k\right)}\end{array}\right]k$$

(57)

Let square matrices of cosine and sine coefficients be denoted by ${\mathbf{A}}_N$ and ${\mathbf{B}}_N$, respectively, so

$${\mathbf{A}}_N=\left[\begin{array}{cccc}{a}_{0,0}& a_{0,1}& \dots & a_{0,N-1}\\ a_{1,0}& a_{1,1}& \dots & a_{1,N-1}\\ ⋮& ⋮& \ddots & ⋮\\ \\ a_{N-1,0}& a_{N-1,1}& \dots & a_{N-1,N-1}\end{array}\right]$$

(58)

$${\mathbf{B}}_N=\left[\begin{array}{cccc}{b}_{0,0}& b_{0,1}& \dots & b_{0,N-1}\\ b_{1,0}& a_{1,1}& \dots & b_{1,N-1}\\ ⋮& ⋮& \ddots & ⋮\\ \\ b_{N-1,0}& b_{N-1,1}& \dots & b_{N-1,N-1}\end{array}\right]$$

(59)

where the entries $a_{m,n}$ and $b_{m,n}$ of the matrices ${\mathbf{A}}_N$ and ${\mathbf{B}}_N$, respectively, are defined in (18). Note that the matrix $({\mathbf{A}}_N-{\mathbf{B}}_{N}i)={\mathbf{F}}_N$, where ${\mathbf{F}}_N$ is the DFT matrix.

In matrix notation, we can write (22)–(25), which describe the right-side QDFT, in the forms

$${\mathbf{y}}_N^{\left(r\right)}={\mathbf{A}}_N{\mathbf{x}}_N^{\left(r\right)}+{\mu }_i{\mathbf{B}}_N{\mathbf{x}}_N^{\left(i\right)}+{\mu }_j{\mathbf{B}}_N{\mathbf{x}}_N^{\left(j\right)}+{\mu }_k{\mathbf{B}}_N{\mathbf{x}}_N^{\left(k\right)}$$

(60)

$${\mathbf{y}}_N^{\left(i\right)}={\mathbf{A}}_N{\mathbf{x}}_N^{\left(i\right)}-{\mu }_i{\mathbf{B}}_N{\mathbf{x}}_N^{\left(r\right)}+{\mu }_j{\mathbf{B}}_N{\mathbf{x}}_N^{\left(k\right)}-{\mu }_k{\mathbf{B}}_N{\mathbf{x}}_N^{\left(j\right)}$$

(61)

$${\mathbf{y}}_N^{\left(j\right)}={\mathbf{A}}_N{\mathbf{x}}_N^{\left(j\right)}-{\mu }_i{\mathbf{B}}_N{\mathbf{x}}_N^{\left(k\right)}-{\mu }_j{\mathbf{B}}_N{\mathbf{x}}_N^{\left(r\right)}+{\mu }_k{\mathbf{B}}_N{\mathbf{x}}_N^{\left(i\right)}$$

(62)

$${\mathbf{y}}_N^{\left(k\right)}={\mathbf{A}}_N{\mathbf{x}}_N^{\left(k\right)}+{\mu }_i{\mathbf{B}}_N{\mathbf{x}}_N^{\left(j\right)}-{\mu }_j{\mathbf{B}}_N{\mathbf{x}}_N^{\left(i\right)}-{\mu }_k{\mathbf{B}}_N{\mathbf{x}}_N^{\left(r\right)}$$

(63)

Figure 2 shows a data flow diagram for calculating the right-side QDFT. The straight lines in the figures indicate the data transmission buses, so that each bus simultaneously transmits all elements of a real-valued vector. By default, we assume that data flows from left to right. Therefore, we do not use arrows, so as not to clutter up the drawings. The points at which the lines converge denote the element-wise summation of the input vectors. The rectangles show the operation of multiplying a vector by a matrix, the designation of which is inscribed inside the rectangle, and the circles show the operation of multiplying all elements of the input vector by a constant inscribed in the circle.

Now, we will develop a similar algorithm for the left-side QDFT. In matrix notation, the equivalents of (29)–(32), which describe the left-side QDFT, are as follows:

$${\mathbf{z}}_N^{\left(r\right)}={\mathbf{A}}_N{\mathbf{x}}_N^{\left(r\right)}+{\mu }_i{\mathbf{B}}_N{\mathbf{x}}_N^{\left(i\right)}+{\mu }_j{\mathbf{B}}_N{\mathbf{x}}_N^{\left(j\right)}+{\mu }_k{\mathbf{B}}_N{\mathbf{x}}_N^{\left(k\right)}$$

(64)

$${\mathbf{z}}_N^{\left(i\right)}={\mathbf{A}}_N{\mathbf{x}}_N^{\left(i\right)}-{\mu }_i{\mathbf{B}}_N{\mathbf{x}}_N^{\left(r\right)}-{\mu }_j{\mathbf{B}}_N{\mathbf{x}}_N^{\left(k\right)}+{\mu }_k{\mathbf{B}}_N{\mathbf{x}}_N^{\left(j\right)}$$

(65)

$${\mathbf{z}}_N^{\left(j\right)}={\mathbf{A}}_N{\mathbf{x}}_N^{\left(j\right)}+{\mu }_i{\mathbf{B}}_N{\mathbf{x}}_N^{\left(k\right)}-{\mu }_j{\mathbf{B}}_N{\mathbf{x}}_N^{\left(r\right)}-{\mu }_k{\mathbf{B}}_N{\mathbf{x}}_N^{\left(i\right)}$$

(66)

$${\mathbf{z}}_N^{\left(k\right)}={\mathbf{A}}_N{\mathbf{x}}_N^{\left(k\right)}-{\mu }_i{\mathbf{B}}_N{\mathbf{x}}_N^{\left(j\right)}+{\mu }_j{\mathbf{B}}_N{\mathbf{x}}_N^{\left(i\right)}-{\mu }_k{\mathbf{B}}_N{\mathbf{x}}_N^{\left(r\right)}$$

(67)

Figure 3 shows the data flow diagram for calculation of the left-side QDFT.

To calculate both QDFTs—the right-side and the left-side—you need to multiply the matrices ${\mathbf{A}}_N$ and ${\mathbf{B}}_N$ by each of the four components of the quaternion-valued input vector ${\mathbf{x}}_N$. Since the components of the input vectors are real values, these operations can be performed by computing the DFT for each component of the input vector and dividing the real and imaginary parts of the result, i.e.,

$${\mathbf{A}}_N{\mathbf{x}}_N^{\left(r\right)}=Re\left({\mathbf{F}}_N{\mathbf{x}}_N^{\left(r\right)}\right),{\mathbf{B}}_N{\mathbf{x}}_N^{\left(r\right)}=-Im\left({\mathbf{F}}_N{\mathbf{x}}_N^{\left(r\right)}\right)$$

(68)

$${\mathbf{A}}_N{\mathbf{x}}_N^{\left(i\right)}=Re\left({\mathbf{F}}_N{\mathbf{x}}_N^{\left(i\right)}\right),{\mathbf{B}}_N{\mathbf{x}}_N^{\left(i\right)}=-Im\left({\mathbf{F}}_N{\mathbf{x}}_N^{\left(i\right)}\right)$$

(69)

$${\mathbf{A}}_N{\mathbf{x}}_N^{\left(j\right)}=Re\left({\mathbf{F}}_N{\mathbf{x}}_N^{\left(j\right)}\right),{\mathbf{B}}_N{\mathbf{x}}_N^{\left(j\right)}=-Im\left({\mathbf{F}}_N{\mathbf{x}}_N^{\left(j\right)}\right)$$

(70)

$${\mathbf{A}}_N{\mathbf{x}}_N^{\left(k\right)}=Re\left({\mathbf{F}}_N{\mathbf{x}}_N^{\left(k\right)}\right),{\mathbf{B}}_N{\mathbf{x}}_N^{\left(k\right)}=-Im\left({\mathbf{F}}_N{\mathbf{x}}_N^{\left(k\right)}\right)$$

(71)

6. 1-D QDFT Proposal

By computing both QDFTs, the main computational cost is related to multiplying the four components of the quaternion-valued input vector by the matrices ${\mathbf{A}}_N$ and ${\mathbf{B}}_N$. This is equivalent to calculating the DFTs for each of the four components and separating the real and imaginary parts, as was described in Section 5. It seems that in order to calculate any QDFT, it is necessary to compute four DFTs. However, it is known that the same operation can be done by performing only two Fourier transforms for two complex-valued input vectors, since each of the two DFTs for real-valued input vectors can be computed with one DFT for the complex-valued input vector. To explain this, let us introduce two real-valued vectors ${\mathbf{x}}_N^{\left(1\right)}$ and ${\mathbf{x}}_N^{\left(2\right)}$ formed from the complex-valued vector ${\mathbf{x}}_N={\mathbf{x}}_N^{\left(1\right)}+i{\mathbf{x}}_N^{\left(2\right)}$. Let ${\mathbf{c}}_N^{\left(1\right)}$ be the DFT output for input ${\mathbf{x}}_N^{\left(1\right)}$, ${\mathbf{c}}_N^{\left(2\right)}$—the DFT output for input ${\mathbf{x}}_N^{\left(2\right)}$, and ${\mathbf{c}}_N$—the DFT output for complex-valued input ${\mathbf{x}}_N$. Although the DFTs outputs for the real-valued inputs are complex-valued, their real and imaginary parts have special symmetries, namely

$$Re\left(c_n^{\left(1\right)}\right)=Re\left(c_{N-n}^{\left(1\right)}\right),Im\left(c_n^{\left(1\right)}\right)=-Im\left(c_{N-n}^{\left(1\right)}\right)$$

(72)

$$Re\left(c_n^{\left(2\right)}\right)=Re\left(c_{N-n}^{\left(2\right)}\right),Im\left(c_n^{\left(2\right)}\right)=-Im\left(c_{N-n}^{\left(2\right)}\right)$$

(73)

for $n=0,1,\dots ,N-1$. It is said that for a real-valued input signal, the real part of its DFT is symmetric and the imaginary part is antisymmetric.

For each coefficient $c_n$ of DFT output for complex-valued signal ${\mathbf{x}}_N={\mathbf{x}}_N^{\left(1\right)}+i{\mathbf{x}}_N^{\left(2\right)}$, we can write

$$c_n=Re\left(c_n\right)+Im\left(c_n\right)i$$

(74)

However, the DFT is a linear transform, so it can be also written as

$$c_n=c_n^{\left(1\right)}+i{c}_n^{\left(2\right)}=Re\left(c_n^{\left(1\right)}\right)+iIm\left(c_n^{\left(1\right)}\right)+i[Re\left(c_n^{\left(2\right)}\right)+iIm\left(c_n^{\left(2\right)}\right)]$$

(75)

Thus,

$$Re\left(c_n\right)=Re\left(c_n^{\left(1\right)}\right)-Im\left(c_n^{\left(2\right)}\right)$$

(76)

and

$$Im\left(c_n\right)=Im\left(c_n^{\left(1\right)}\right)+Re\left(c_n^{\left(2\right)}\right)$$

(77)

Analogously,

$$Re\left(c_{N-n}\right)=Re\left(c_{N-n}^{\left(1\right)}\right)-Im\left(c_{N-n}^{\left(2\right)}\right)$$

(78)

and

$$Im\left(c_{N-n}\right)=Im\left(c_{N-n}^{\left(1\right)}\right)+Re\left(c_{N-n}^{\left(2\right)}\right)$$

(79)

Taking into account properties (72) and (73), (78) can be written in the form

$$Re\left(c_{N-n}\right)=Re\left(c_n^{\left(1\right)}\right)+Im\left(c_n^{\left(2\right)}\right)$$

(80)

and (79) as

$$Im\left(c_{N-n}\right)=-Im\left(c_n^{\left(1\right)}\right)+Re\left(c_n^{\left(2\right)}\right)$$

(81)

By adding and subtracting the sides of (80) and (76), we obtain, respectively

$$Re\left(c_n^{\left(1\right)}\right)=\frac{Re\left(c_n\right)+Re\left(c_{N-n}\right)}{2}$$

(82)

and

$$Im\left(c_n^{\left(2\right)}\right)=\frac{Re\left(c_{N-n}\right)-Re\left(c_n\right)}{2}$$

(83)

Similarly, by adding and subtracting the sides of (77) and (81) we obtain, respectively

$$Re\left(c_n^{\left(2\right)}\right)=\frac{Im\left(c_n\right)+Im\left(c_{N-n}\right)}{2}$$

(84)

and

$$Im\left(c_n^{\left(1\right)}\right)=\frac{Im\left(c_n\right)-Im\left(c_{N-n}\right)}{2}$$

(85)

We introduce the denotation ${}^{\uparrow }$, which means the signal reversed in time. So, for ${\mathbf{c}}_N={[c_0,c_1,\dots ,c_{N-1}]}^T$, we have

$${\mathbf{c}}_N^{\uparrow }=\left[\begin{array}{c}{c}_0\\ c_{N-1}\\ ⋮\\ c_1\end{array}\right]={\mathbf{P}}_N{\mathbf{c}}_N$$

(86)

where the matrix ${\mathbf{P}}_N$, which is responsible for reversing signal in time, has the form

$${\mathbf{P}}_N=\left[\begin{array}{ccccc}1& 0& \dots & 0& 0\\ 0& 0& \dots & 0& 1\\ 0& 0& \dots & 1& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 1& \dots & 0& 0\end{array}\right]$$

(87)

Then, (82) written in vector notation has the form

$$Re\left({\mathbf{c}}_N^{\left(1\right)}\right)=\frac{1}{2}[Re\left({\mathbf{c}}_N\right)+Re\left({\mathbf{c}}_N^{\uparrow }\right)]$$

(88)

Analogously, (85) written in vector notation has the form

$$Im\left({\mathbf{c}}_N^{\left(1\right)}\right)=\frac{1}{2}[Im\left({\mathbf{c}}_N\right)-Im\left({\mathbf{c}}_N^{\uparrow }\right)]$$

(89)

Similarly, (84) written in vector notation has the form

$$Re\left({\mathbf{c}}_N^{\left(2\right)}\right)=\frac{1}{2}[Im\left({\mathbf{c}}_N\right)+Im\left({\mathbf{c}}_N^{\uparrow }\right)]$$

(90)

In the end, (83) written in vector notation has the form

$$Im\left({\mathbf{c}}_N^{\left(2\right)}\right)=\frac{1}{2}[Re\left({\mathbf{c}}_N^{\uparrow }\right)-Re\left({\mathbf{c}}_N\right)]$$

(91)

Going back to QDFT transforms, the quaternion-valud input vector ${\mathbf{x}}_N$ has four real-valued component ${\mathbf{x}}_N^{\left(r\right)}$, ${\mathbf{x}}_N^{\left(i\right)}$, ${\mathbf{x}}_N^{\left(j\right)}$, ${\mathbf{x}}_N^{\left(k\right)}$. Let us denote their complex Fourier transforms by ${\mathbf{c}}_N^{\left(r\right)}$, ${\mathbf{c}}_N^{\left(i\right)}$, ${\mathbf{c}}_N^{\left(j\right)}$, ${\mathbf{c}}_N^{\left(k\right)}$. They are all complex-valued, and

$${\mathbf{A}}_N{\mathbf{x}}_N^{\left(r\right)}=Re\left({\mathbf{c}}_N^{\left(r\right)}\right),{\mathbf{B}}_N{\mathbf{x}}_N^{\left(r\right)}=-Im\left({\mathbf{c}}_N^{\left(r\right)}\right)$$

(92)

$${\mathbf{A}}_N{\mathbf{x}}_N^{\left(i\right)}=Re\left({\mathbf{c}}_N^{\left(i\right)}\right),{\mathbf{B}}_N{\mathbf{x}}_N^{\left(i\right)}=-Im\left({\mathbf{c}}_N^{\left(i\right)}\right)$$

(93)

$${\mathbf{A}}_N{\mathbf{x}}_N^{\left(j\right)}=Re\left({\mathbf{c}}_N^{\left(j\right)}\right),{\mathbf{B}}_N{\mathbf{x}}_N^{\left(j\right)}=-Im\left({\mathbf{c}}_N^{\left(j\right)}\right)$$

(94)

$${\mathbf{A}}_N{\mathbf{x}}_N^{\left(k\right)}=Re\left({\mathbf{c}}_N^{\left(k\right)}\right),{\mathbf{B}}_N{\mathbf{x}}_N^{\left(k\right)}=-Im\left({\mathbf{c}}_N^{\left(k\right)}\right)$$

(95)

From the components of the quaternion input vector ${\mathbf{x}}_N$, we can create two complex inputs vectors ${\tilde{\mathbf{x}}}_N={\mathbf{x}}_N^{\left(r\right)}+{\mathbf{x}}_N^{\left(i\right)}i$ and ${\widehat{\mathbf{x}}}_N={\mathbf{x}}_N^{\left(j\right)}+{\mathbf{x}}_N^{\left(k\right)}i$. Then, we calculate the complex discrete Fourier transforms for these, i.e., ${\tilde{\mathbf{c}}}_N$ and ${\widehat{\mathbf{c}}}_N$, using the fast algorithm, and separate their real and imaginary parts. Then, according to (88)–(91), for the first complex vector ${\tilde{\mathbf{x}}}_N$, we obtain

$$Re\left({\mathbf{c}}_N^{\left(r\right)}\right)=\frac{1}{2}[Re\left({\tilde{\mathbf{c}}}_N\right)+Re\left({\tilde{\mathbf{c}}}_N^{\uparrow }\right)]$$

(96)

$$Im\left({\mathbf{c}}_N^{\left(r\right)}\right)=\frac{1}{2}[Im\left({\tilde{\mathbf{c}}}_N\right)-Im\left({\tilde{\mathbf{c}}}_N^{\uparrow }\right)]$$

(97)

$$Re\left({\mathbf{c}}_N^{\left(i\right)}\right)=\frac{1}{2}[Im\left({\tilde{\mathbf{c}}}_N\right)+Im\left({\tilde{\mathbf{c}}}_N^{\uparrow }\right)]$$

(98)

$$Im\left({\mathbf{c}}_N^{\left(i\right)}\right)=\frac{1}{2}[Re\left({\tilde{\mathbf{c}}}_N^{\uparrow }\right)-Re\left({\tilde{\mathbf{c}}}_N\right)]$$

(99)

and, similarly, for the second complex-valued vector ${\widehat{\mathbf{x}}}_N$, we obtain

$$Re\left({\mathbf{c}}_N^{\left(j\right)}\right)=\frac{1}{2}[Re\left({\widehat{\mathbf{c}}}_N\right)+Re\left({\widehat{\mathbf{c}}}_N^{\uparrow }\right)]$$

(100)

$$Im\left({\mathbf{c}}_N^{\left(j\right)}\right)=\frac{1}{2}[Im\left({\widehat{\mathbf{c}}}_N\right)-Im\left({\widehat{\mathbf{c}}}_N^{\uparrow }\right)]$$

(101)

$$Re\left({\mathbf{c}}_N^{\left(k\right)}\right)=\frac{1}{2}[Im\left({\widehat{\mathbf{c}}}_N\right)+Im\left({\widehat{\mathbf{c}}}_N^{\uparrow }\right)]$$

(102)

$$Im\left({\mathbf{c}}_N^{\left(k\right)}\right)=\frac{1}{2}[Re\left({\widehat{\mathbf{c}}}_N^{\uparrow }\right)-Re\left({\widehat{\mathbf{c}}}_N\right)]$$

(103)

The data flow diagrams for our method of QDFT calculation are shown in Figure 4 for the right-side QDFT and in Figure 5 for the left-side QDFT.

It should be noted that the proposed method, unlike the method using symplectic decomposition, does not require changing the basis in the pure quaternions subspace and, consequently, calculating the new basis vectors and change-of-basis matrix.

7. Computational Complexity Discussion

We will now compare the computational complexity of the different ways to calculate a QDFT, taking into account the number of multiplications and additions of real numbers needed to compute this transform (these numbers are the same for the right-side and left-side transforms). Since, in each case, the ordinary FFT algorithm is used, we denote by $m_{FFT}\left(N\right)$ the number of real multiplications necessary to compute the DFT for a complex/real signal with N samples. Similarly, let $a_{FFT}\left(N\right)$ denote the number of additions of real numbers needed to determine the FFT of such a signal. In this comparison, we will omit the method of calculating QDFTs resulting directly from the definitions (16) or (26) because its computational complexity is much higher.

Let us start with the most basic calculation method resulting from (60)–(63) or (64)–(67), respectively. Since, according to (68), calculating the products of ${\mathbf{A}}_N{\mathbf{x}}_N^{\left(r\right)}$ and ${\mathbf{B}}_N{\mathbf{x}}_N^{\left(r\right)}$ is equivalent to finding the FFT for the ${\mathbf{x}}_N^{\left(r\right)}$ part of the quaternion vector ${\mathbf{x}}_N$, and similarly, for the other ${\mathbf{x}}_N^{\left(i\right)}$, ${\mathbf{x}}_N^{\left(j\right)}$, ${\mathbf{x}}_N^{\left(k\right)}$ parts of the quaternion input vector, we need to calculate four FFTs, which requires $4m_{FFT}\left(N\right)$ multiplications and $4a_{FFT}\left(N\right)$ additions of real numbers. Then, each of the products ${\mathbf{B}}_N{\mathbf{x}}_N^{\left(r\right)}$, ${\mathbf{B}}_N{\mathbf{x}}_N^{\left(i\right)}$, ${\mathbf{B}}_N{\mathbf{x}}_N^{\left(j\right)}$, ${\mathbf{B}}_N{\mathbf{x}}_N^{\left(k\right)}$ has to be multiplied by three real coefficients ${\mu }_i$, ${\mu }_j$, and ${\mu }_k$, which gives $3·4·N=12N$ multiplications of real numbers. The last step is to add the four components of the sum. Since each sum consists of four vectors, it takes $3·N$ additions to calculate the sum. There are four such sums, so we have $4·3·N=12N$ additions of real numbers. Summarizing the above considerations, in order to calculate the DQFT according to this method, $m_1\left(N\right)$ multiplications and $a_1\left(N\right)$ additions of real numbers should be performed, where

$$m_1\left(N\right)=4m_{FFT}\left(N\right)+12N$$

(104)

and

$$a_1\left(N\right)=4a_{FFT}\left(N\right)+12N$$

(105)

The next method in the comparison is the method based on symplectic decomposition. First, all samples of the quaternion signal should be expressed using the new imaginary units ${\mu }_1$, ${\mu }_2$, ${\mu }_3$, which corresponds to the multiplication of the matrix inverse to the change-of-basis matrix by the vector containing the imaginary components of the sample according to (39). Such an operation for one sample requires nine multiplications and six additions of real numbers. For N samples of the quaternion signal, we get $9N$ multiplications and $6N$ additions of real numbers. Then, a single ordinary complex FFT and a single IFFT (without dividing the result by N), or two complex FFTs, should be calculated, according to (48), (49) or (53), (54), which require $2m_{FFT}\left(N\right)$ multiplications and $2a_{FFT}\left(N\right)$ additions of real numbers. The resulting quaternion output signal is expressed using the new imaginary units ${\mu }_1$, ${\mu }_2$, ${\mu }_3$, and we want it to be represented by standard imaginary units i, j, k. This corresponds to the multiplication of the change-of-basis matrix by the vector containing the imaginary components of the resulting vector according to (51). To do this for N samples of the quaternion output signal, $9N$ multiplications and $6N$ additions of real numbers is needed. To recap the method based on symplectic decomposition, calculating the QDFT according to this method requires $m_2\left(N\right)$ multiplications and $a_2\left(N\right)$ additions of real numbers, where

$$m_2\left(N\right)=2m_{FFT}\left(N\right)+18N$$

(106)

and

$$a_2\left(N\right)=2a_{FFT}\left(N\right)+12N$$

(107)

Finally, we evaluate the numbers of arithmetic operations needed for our method of QDFT calculation. Our method is also based on (60)–(63) for the right-side QDFT or (64)–(67) for the left-side QDFT, but taking into account relations (92)–(95), the products ${\mathbf{A}}_N{\mathbf{x}}_N^{\left(r\right)}$, ${\mathbf{B}}_N{\mathbf{x}}_N^{\left(r\right)}$, ${\mathbf{A}}_N{\mathbf{x}}_N^{\left(i\right)}$, ${\mathbf{B}}_N{\mathbf{x}}_N^{\left(i\right)}$, ${\mathbf{A}}_N{\mathbf{x}}_N^{\left(j\right)}$, ${\mathbf{B}}_N{\mathbf{x}}_N^{\left(j\right)}$ and ${\mathbf{A}}_N{\mathbf{x}}_N^{\left(k\right)}$, ${\mathbf{B}}_N{\mathbf{x}}_N^{\left(k\right)}$ are calculated from (96)–(103). This method requires the calculation of two FFTs, hence $2m_{FFT}\left(N\right)$ multiplications and $2a_{FFT}\left(N\right)$ additions of real numbers is needed. Then, the real/imaginary parts of the FFT results are added/subtracted (some after reordering their elements, which corresponds to multiplying the matrix ${\mathbf{P}}_N$ by these elements, as was shown in Figure 3 and Figure 4—this does not require any arithmetic operations) and divided by 2. This may seem to require $8N$ additions and $8N$ multiplications of real numbers. In fact, ${\mathbf{c}}_N^{\left(r\right)}$, ${\mathbf{c}}_N^{\left(i\right)}$, ${\mathbf{c}}_N^{\left(j\right)}$, ${\mathbf{c}}_N^{\left(k\right)}$ are Fourier transforms of real-valued vectors, so they have symmetry, as in (72)–(73), and only half of their entries have to be calculated. Thus, the number of real additions reduces to $4N$. Multiplying by the factor $1/2$ in (96)–(103) does not require any multiplication operation. It is just a simple bit shift to the right. Then, as in the first method, each of the products ${\mathbf{B}}_N{\mathbf{x}}_N^{\left(r\right)}$, ${\mathbf{B}}_N{\mathbf{x}}_N^{\left(i\right)}$, ${\mathbf{B}}_N{\mathbf{x}}_N^{\left(j\right)}$, ${\mathbf{B}}_N{\mathbf{x}}_N^{\left(k\right)}$ has to be multiplied by three real coefficients ${\mu }_i$, ${\mu }_j$, and ${\mu }_k$, which gives $12N$ multiplications of real numbers. The last step is to add the four components of each sum in (60)–(63) for the right-side QDFT or (64)–(67) for the left-side QDFT, as in the first method. This requires $12N$ additions of real numbers. Summarizing the above considerations, in order to calculate any of the QDFT according to our method $m_3\left(N\right)$ multiplications and $a_3\left(N\right)$ additions of real numbers have to be performed, where

$$m_3\left(N\right)=2m_{FFT}\left(N\right)+12N$$

(108)

and

$$a_3\left(N\right)=2a_{FFT}\left(N\right)+16N$$

(109)

Table 1. Number of multiplications and additions of real numbers needed to calculate any 1-D QDFT for a quaternion signal of N samples.

Type of Operation	Simplest Method	Symplectic Method	Proposed Method
×	$4m_{FFT}\left(N\right)+12N$	$2m_{FFT}\left(N\right)+18N$	$2m_{FFT}\left(N\right)+12N$
+	$4a_{FFT}\left(N\right)+12N$	$2a_{FFT}\left(N\right)+12N$	$2a_{FFT}\left(N\right)+16N$

presents estimates for the numbers of multiplications and additions of real numbers necessary to determine any of the QDFTs in the three presented methods of QDFT calculation. It is easy to see that although our method has more additions than method 2, it requires fewer multiplications. Furthermore, the total number of arithmetic operations in our method is the smallest among the three compared methods for calculating the QDFT.

In terms of the number of arithmetic operations, our method is very similar to the symplectic decomposition method. To compare the 1-D QDFT computation times of these two methods, we implemented them in Matlab R2022b on a computer with an Intel(R) Core(TM) i5-7400 CPU and 8 GB RAM. For quaternion signals of different lengths N, which are powers of 2, and 100,000 repetitions of transform calculation, we obtained the times, in seconds, presented in

Table 2. Calculation times of 1-D QDFTs for signals of length N for the method using symplectic decomposition and the proposed method, for 100,000 repetitions.

N	Time for Simplectic Method [s]	Time for Proposed Method [s]
128	1.43	1.19
256	2.32	1.87
512	3.33	2.60
1024	5.40	4.49
2048	9.14	7.64
4096	19.36	14.90
8192	35.63	28.38

.

The dependencies in Table 2 are not an exact reflection of the dependencies in Table 1, because the calculation time also includes the time needed to transfer data from and to memory, which Table 1 does not take into account. Despite this, it can be observed from Table 2 that the times needed to calculate the 1-D QDFT with the proposed method are shorter than for the method using symplectic decomposition. However, it should be remembered that the time of computation depends very much on the way the code was written and on the implementation platform.

8. Conclusions

The article presents a new method of 1-D QDFT calculation. The proposed method allows a reduction in the number of multiplications in the calculation of the QDFT. Reducing the number of multipliers is especially important when designing dedicated on-board VLSI processors, since minimizing the number of multipliers required also reduces power dissipation and lowers the cost of implementing the entire system.

Furthermore, although the number of additions of real numbers in our method increases, the total number of arithmetic operations is still lower than in other ways of calculating this transform. This leads to a reduction in computational costs, which makes the proposed solution also preferable for its implementation on a general-purpose computer.

It should be added that the newly proposed method can be used to compute both the left-side and right-side QDFTs.