Discrete octonion Fourier transform and the analysis of discrete 3-D data

The purpose of this article is to further develop the theory of octonion Fourier transformations (OFT), but from a different perspective than before. It follows the earlier work by Błaszczyk and Snopek, where they proved a few essential properties of the OFT of real-valued functions of three continuous variables. The research described in this article applies to discrete transformations, i.e. discrete-space octonion Fourier transform (DSOFT) and discrete octonion Fourier transform (DOFT). The described results combine the theory of Fourier transform with the analysis of solutions for difference equations, using for this purpose previous research on algebra of quadruple-complex numbers. This hypercomplex generalization of the discrete Fourier transformation provides an excellent tool for the analysis of 3-D discrete linear time-invariant (LTI) systems and 3-D discrete data.


Introduction
Research in the field of signal and system processing is not only an analysis of the functions of a continuous variable (interpreted as time or space) and their representations in the frequency domain, but also in the field of discrete variables, which we naturally obtain as a result of the process of digitization (i.e. in connection with sampling of signals). As in the case of classical signal processing, so the discrete counterpart of this theory has so far mainly focused on signals with real and complex values, as well as their complex spectra. In recent years, Communicated by Apala Majumdar.  however, more and more works have started to appear, which authors use in their research hypercomplex algebras, among others quaternions and octonions (Brackx et al. 2013;Hahn and Snopek 2016;Lian 2019;Snopek 2015;Wang et al. 2017). The area of applications is focused so far on the study of neural networks (Popa 2016(Popa , 2018, analysis of color and multispectral images (Ell et a. 2014;Gao and Lam 2014;Gomes et al. 2017;Grigoryan and Agaian 2018;Lazendić et al. 2018a, b;Sheng et al. 2018), as well as the biomedical signals processing (Delsuc 1988;Klco et al. 2017).
Octonions form a non-associative and a non-commutative algebra (Baez 2002).
In previous studies, we focused on discussing the properties of the octonion Fourier transform (OFT) of real-valued functions of three variables (and in the case of some properties, we extended this to octonion-valued functions) (Błaszczyk , 2019(Błaszczyk , 2020Błaszczyk and Snopek 2017). OFT is given by a formula x 3 ) · e −2π e 1 f 1 x 1 · e −2π e 2 f 2 x 2 · e −2π e 4 f 3 x 3 dx 1 dx 2 dx 3 , where e 1 , e 2 and e 4 are three of the seven octonion imaginary units and the octonion exponential function is defined by the power series (Błaszczyk 2020), similarly as for the complex numbers and quaternions. In previous investigations (Błaszczyk 2020;Błaszczyk and Snopek 2017) it has been shown that (1) is well defined (i.e. the inverse transform theorem is valid) and some properties of the OFT, analogous to the properties of the complex and quaternion Fourier transform (e.g. Hermitian symmetry, shift theorem, Plancherel-Parseval theorems and derivative theorem) were also derived.

Fig. 2 Multiplication rules in F
A major difficulty in these works was the lack of significant properties of octonion multiplication, i.e. commutativity and associativity. Non-commutativity is also encountered with Fourier transforms based on quaternion algebra and (in general) Clifford algebras (Brackx et al. 2013). In order to deal with these problems, we have defined the algebra of quadruplecomplex numbers (based on the double-complex numbers introduced by Ell (1993); Kurman (1958)). Like octonions, we define the algebra of order 8 over the field of real numbers and each element of this algebra will be identified with the 8-tuple of real numbers (Błaszczyk 2019(Błaszczyk , 2020, i.e. Addition in F is defined in a classical way-element-wise, and multiplication (denoted as ) rules are shown in Fig. 2. We can see that imaginary units in F follow the analogous (although significantly different) rules as the octonion multiplication, i.e.
It is worth pointing out that, just as in the case of octonions, algebra F is generated by three imaginary units, i.e. E 1 , E 2 and E 4 . Each other imaginary unit can be obtained by multiplying these three units, i.e. E 3 = E 1 E 2 , E 5 = E 1 E 4 , E 6 = E 2 E 4 and E 7 = E 1 E 2 E 4 . However, unlike octonions, the three basic imaginary units of algebra F are commutative. The multiplication in F is commutative and associative, however not every non-zero element of F has its inverse, which is a property common to Clifford algebras. This algebra has already appeared in the literature, in particular in Felsberg et al. (2001) (where it was denoted as hypercomplex algebra H 4 ). By using operations in the quadruple algebra the notation of formulas for determining many OFT properties could be reduced to simple equations known from classical theory. We discussed that in detail in Błaszczyk (2019Błaszczyk ( , 2020. From the point of view of numerical calculations and practical applications, it is also important to develop discrete equivalents of the transformations under consideration and their properties. In the case of the classical Fourier transformation, this subject is well known, as in the case of quaternions (Bahri and Surahman 2013;Ell et a. 2014;Sangwine and Bihan 2005;Sangwine 1997). In the case of octonions, references to the analysis of 1-and 2-dimensional signals appear in the literature (Grigoryan and Agaian 2018), but there is still no definition of discrete octonion Fourier transform of 3-dimensional signals.
The paper is organized as follows. In Sect. 2 we recall previous results concerning the octonion Fourier transform and introduce some important properties. In Sect. 3 we focus on the discrete-time LTI systems, which leads to the definition of the discrete-space octonion Fourier transform. Its well-posedness and some of its properties are the main part of Sect. 4. Sections 5 and 6 are devoted to the discrete octonion Fourier transform-its definition, properties and implementation aspects. The implementation of the DOFT algorithm and the results of the first tests are also presented there. The paper is concluded in Sect. 7 with a short discussion of obtained results.

Octonion Fourier transform and some properties
As mentioned earlier, the OFT of the real-valued function u : R 3 → R is given by the formula (1). In order for the integral (1) to exist, it is necessary for the function u to be at least integrable. However, in general, conditions of existence of the OFT are the same as for the classical (complex) Fourier transform (Błaszczyk 2020). Choice and order of imaginary units in the exponents is not accidental (see Błaszczyk and Snopek 2017). In our previous studies (see Błaszczyk and Snopek 2017), we proved the inverse formula (for continuous functions u : R 3 → R such that both u and its OFT are integrable): and multiplication is performed from left to right. We have also proved some important features, among which there are those that will be useful in the analysis of discrete-variable signals (Błaszczyk 2020;Błaszczyk and Snopek 2017). Below, we quote their wording, assuming in each of these statements that the OFTs are well-defined and invertible. Let U be the OFT of the real-valued function u and let u x i denote the partial derivative of u with respect to x i .
Theorem 2 Let U ∂ x 1 , U ∂ x 2 and U ∂ x 3 denote the OFTs of u x 1 , u x 2 and u x 3 , respectively. Then It seems that in the above form these statements are of little use. It is worth noting, however, that on the one hand we have a theorem on Hermitian symmetry (Błaszczyk and Snopek 2017), and on the other hand we can reformulate the given formulas using the operation . Then these theorems take the form known from classical theory.
Corollary 2 Let U ∂ x 1 , U ∂ x 2 and U ∂ x 3 denote the OFTs of u x 1 , u x 2 and u x 3 , respectively. Then

Remark 1
The claims of the Corollaries 1 and 2 should be understood in a specific way. We will show it on the example of the last formula in Corollary 2. The OFT of the real-valued function u can be expressed as U = U 0 + U 1 e 1 + U 2 e 2 + U 3 e 3 + U 4 e 4 + U 5 e 5 + U 6 e 6 + U 7 e 7 , where U 0 , . . . , U 7 : R 3 → R are real-valued functions (we omit the argument of the function so as not to lose readability of the formula). Then U · e 4 = −U 4 − U 5 e 1 − U 6 e 2 − U 7 e 3 + U 1 e 4 + U 2 e 5 + U 3 e 6 + U 4 e 7 .
On the other hand we can treat the octonion-valued function U as the 8-tuple of real-valued functions, which allows us to write this function as a F-valued function: which, again treating this expression as 8-tuple of real-valued functions, gives the statement of the Corollary. Other expressions should be understood analogously.
Proofs of these claims are based on direct calculations and we omit the details here. It should be remembered that the notation related to the operation serves only to increase the readability of performed operations and facilitate their interpretation. More importantly, it allows the direct use of the OFT for the analysis of LTI systems, which are described both by partial differential equations and difference equations (of three variables).

Discrete-time LTI systems
Consider linear time-invariant stationary system of three variables. We know from the classical signal and system theory that such systems are described by their impulse responses h : R 3 → R (also sometimes called Green functions) and then the input-output relation of signals u : R 3 → R and v : R 3 → R, respectively, is given by the formula: The function at the output of the system is therefore a convolution of the function on the input with the impulse response, which is schematically presented in the block form as in Fig. 3.
From the convolution-multiplication duality theorem for the classic Fourier transformation, it immediately follows that the following equality occurs:  and function H = F CFT {h} is called the frequency response of the system. In Błaszczyk (2019Błaszczyk ( , 2020 we showed that the similar relation is valid also in case of octonion Fourier transform: where H OFT is the octonion Fourier transform of the impuls response h (and therefore called the octonion frequency response) and we use the multiplication in the algebra of quadruplecomplex numbers (in the sense described in Remark 1).
Thanks to these relations (and also Corollary 2), we can see that every linear partial differential equation with constant coefficients can be reduced to the algebraic equation (in the sense of multiplication in F). In Błaszczyk (2019Błaszczyk ( , 2020 we presented a few examples showing the possibilities of using this notation.
Analogous reasoning can be performed for discrete-time systems. They are described mostly by difference equations, i.e. equations of the form: Using the octonion Fourier transform and Theorem 1, we can write this equality in a different form. However, only by composing Corollary 1 (and treating octonions as 8-tuples of real numbers, and therefore as elements of F) we can reduce this relation to that which we know from classical theory: It can be shown with direct calculations that where the multiplication in the octonion algebra is performed from left to right. The above equality should be understood in the sense described in Remark 1, i.e. as the equality of 8-tuples of real numbers. Therefore we have where expressions in brackets are some octonion counterparts to discrete Fourier transforms of vectors a = (a i ), b = (b j ) indexed by 3-D discrete variables.

Discrete-space octonion Fourier transform
The above considerations lead to the definition of the octonion counterpart of the discrete Fourier transform of real-valued sequences.
Definition 1 Let a : N 3 → R be a sequence indexed by a 3-D discrete variable and let Octonion Fourier transform of the sequence a is defined by the formula where multiplication is performed from left to right.
The above definition is a three-dimensional equivalent of the Fourier transformation of discrete time (discrete-time Fourier transform-DTFT), in relation to which the given formula can be abbreviatied as DSOFT (discrete-space octonion Fourier transform). Like its classic equivalent, DSOFT is a periodic function with relation to each variable, with a period of 1. Using the methods presented in the proof of the inverse theorem (Błaszczyk and Snopek 2017), the following formula can be derived for the inverse transform.

Theorem 3 Let a : N 3 → R be a sequence indexed by a 3-D discrete variable and let
where multiplication is performed from left to right.
This theorem can be generalized to octonion-valued sequences, proving it with the same methods as in Błaszczyk ( , 2020. From simple calculations the proof of the following theorem on the transformation of the rescaled function follows.
Proof From straightforward calculations we have that which concludes the proof.
It is worth noting that proofs of other properties of the octonion Fourier transformation (of functions of the continuous variable) proceeded in the same way as proofs of equivalents of these properties for the classic Fourier transform. Differences in statements of those theorems resulted only from the properties of multiplication of octonions, and, as a consequence, of operations on exponential functions in the kernel of transformation. Therefore, the equivalents of these properties for DSOFT will look the same, and their proofs will be very similar. We will therefore omit the details and quote only the statements.
Theorem 5 Let a = (a i ) : N 3 → R and let A denote the DSOFT of a. Moreover, let f 0 ∈ R and denote a cos,k = (a i ), a cos,k i = a i · cos(2π f 0 i k ), and A cos,k as the DSOFT of a cos,k , k = 1, 2, 3. Then Theorem 6 Let a = (a i ) : N 3 → R and let A denote the DSOFT of a. Moreover, let f 0 ∈ R and denote a sin,k = (a i ), a sin,k i = a i · sin(2π f 0 i k ), and A cos,k as the DSOFT of a sin,k , k = 1, 2, 3. Then From the reasoning in the previous section, the shift theorem also immediately follows.
Theorem 7 Let a = (a i ) : N 3 → R and let A denote the DSOFT of a. Moreover, let α, β, γ ∈ Z and denote As we mentioned earlier, deriving the DSOFT definition, it can be used to analyze discrete systems described by difference equations. It is also easy to see that these methods will be used in the analysis of finite difference schemes for partial differential equations. The theory considered so far mainly used classic Fourier transforms (discrete), which were applied to a variable interpreted as space (one-or two-dimensional), but by defining an octonion transformation we can try to transform the whole scheme, both for time and space.

Discrete octonion Fourier transform
In the previous section we considered signals (sequences) of infinite length. However, in practice, finite-length signals are usually found, which, as in the classical case, leads to the definition of discrete octonion Fourier transform. Similarly to the complex and quaternion case, the following definition has its basis in the discretization of the frequency domain in discrete-space octonion Fourier transform.
, be a finite-length sequence indexed by a 3-D discrete variable and let a = (a n ), n = (n 1 , n 2 , n 3 ). Discrete octonion Fourier transform (DOFT) A OFT = (A OFT,k ) of a is defined by the formula a (n 1 ,n 2 ,n 3 ) · e −e 1 2π k 1 n 1 /N 1 · e −e 2 2π k 2 n 2 /N 2 · e −e 4 2π k 3 n 3 /N 3 , (2) and multiplication is performed from left to right. This is the direct octonion equivalent of the 3-D discrete Fourier transformation. The formula for inverse transformation is proved analogously to the corresponding formula for inverse OFT (Błaszczyk 2020). In particular, the following equality is used: valid for every octonion imaginary unit e i , i = 1, . . . , 7. We omit the proof and only quote the statement of the inverse-transform theorem.
, be a finite-length sequence indexed by a 3-D discrete variable and let a = (a n ), n = (n 1 , n 2 , n 3 ). If A OFT = (A OFT,k ) is the DOFT of a, then a n = 1 N 1 N 2 N 3 A OFT,(k 1 ,k 2 ,k 3 ) · e e 4 2π k 3 n 3 /N 3 · e e 2 2π k 2 n 2 /N 2 · e e 1 2π k 1 n 1 /N 1 , where multiplication is performed from left to right.
From a computational point of view, consideration should be given to the possibility of implementing a fast DOFT calculation version. In the quaternion case, various solutions to this problem have been proposed-direct implementation of the quaternion version of the Fast Fourier Transform algorithm (FFT) and the use of the complex (original) version of this algorithm. The latter option turned out to be better and less computationally demanding (Ell et a. 2014).
In the case of octonions, one can do the same. In Błaszczyk (2020) we proved that the octonion Fourier transformation of octonion-valued functions can be calculated using the classical Fourier transformation. This was formulated in the following statement (for proof see Błaszczyk 2020).
It is important to notice that the field of complex numbers C used in the abovementioned theorem is the specific subfield of the octonion algebra, i.e. C = {x 0 +x 1 e 1 ∈ O : x 0 , x 1 ∈ R}. Consequently, only the imaginary units e 2 and e 4 appear in the thesis. Equivalent versions of the above theorem can be derived which take into account other pairs of imaginary units.
Analogous theorems for calculating the inverse transformation are also known (Błaszczyk 2020, Theorem 5 and Corollary 2). Due to the long formulation of the theorem, we omit it here. These claims remain true also for discrete transformations. The operation of changing the sign of a variable should be understood in the sense of modulo operations, as in the classic case-DOFT, just like DFT, can be treated as a periodic function.
Thanks to this, it is possible to use all the advantages of the FFT algorithm, with a small additional calculation effort-the octonion FFT algorithm for functions with octonion values requires the calculation of four transforms of different functions with complex values.
Tools for operations on hypercomplex numbers are available in many programming environments, including in MATLAB®. MATLAB® environment is one of the more popular tools supporting the work of engineers and packages expanding programming capabilities in this environment appear quite often. The qtfm package developed by the team of S. Sangwine and N. Le Bihan Sangwine and Bihan (2005) focuses on numerical calculations in quaternion algebra (not only basic arithmetic operations are available, but also highly developed tools for calculating quaternion Fourier transforms), but on the other hand, it also allows simple operations in octonion algebra. However, it lacks more advanced features that would give users the opportunity to calculate octonion Fourier transforms.
The qtfm package has been extended by, among others, five additional functions -doft3 and idoft3 (forward and inverse DOFT using the direct formula), -offt3 and iofft3 (forward and inverse DOFT using the FFT algorithm), fftreflection (symmetrical reflection of functions relative to the indicated axes).
The functions have been implemented using the qtfm package syntax. Octonions are entered in it, among others as octonion(r0, r1, r2, r3, r4, r5, r6, r7), where r0, …, r7 are numbers or matrixes of numbers. It is also possible to automatically project real numbers, complex numbers and quaternions to the appropriate octonions, using the fact that they are created in the Cayley-Dickson construction (this is done using the available function cd).
The implementation of the forward and inverse discrete octonion Fourier transformation was performed in two ways-using the direct formula (2) and using Theorem 9, where first the classical (complex) Fourier transforms are calculated, and then they are combined in an appropriate manner. In the case of the direct algorithm, in order to avoid problems with multiplication of octonions, the matrix representation of multiplication presented in Tian (2000) was used. One of the important elements of the octonion FFT algorithm is the ability to symmetrically reflect functions with relation to the respective axes. This has been implemented in the fftreflection(X, A) function, where X is a mirrored function (matrix), and A is a vector that indicates which variables to reflect (e.g. [1,3] means that X is mirrored with respect to 1st and 3rd variable). It uses the convention that MATLAB® adopts with the FFT algorithm -zero frequencies always appear on the first coordinates of the matrix (in each dimension), and the discrete Fourier transform is a periodic function.
Implemented functions can be found in the GitHub repository: https://github.com/blaszczykl/matlab-octonions and the code is also included in Appendix 1 (Listings 1-5). This is the original code developed by the author, however, it requires the qtfm package to work correctly (it is based on the syntax used in this package).
The correctness of implemented functions can be tested in two ways. First, it was checked whether applying the reverse transformation to the result of the forward transformation returned the original matrix. This check was done for both versions of the algorithm by generating random N × N × N octonion matrices, where N = 4 and each coordinate was generated from a uniform distribution over the interval [0, 1] (for details see Listing 7). In both cases, the maximum relative error did not exceed 2.2 · 10 −11 %, which can be treated as a numerical approximation error.
It is worth noting the (expected) problem of calculation time for various quantities of N . In the case of the direct algorithm, the execution time is proportional to N 6 , while for the FFT algorithm the time is proportional to N 3 log N . This is not a surprising feature, it also applies to the classical DFT, however, due to operations on octonions, the time needed to perform calculations may be larger than in the case of the classical one. It can be seen that this issue is important.
To check whether an error was made when implementing the algorithm using Theorem 9 and whether it actually returns a discrete OFT, it was assumed that the direct algorithm returns the correct result and then it was compared with the result returned by the octonion FFT. In this case, the error did not exceed 1.5 · 10 −12 %, which suggests that the implementation is correct.

Symmetry properties of the DOFT
As in the case of continuous OFT, most discrete classical equivalents of Fourier transform properties can be proved for discrete OFT. One of the first results for continuous OFT was to show the equivalent of Hermitian symmetry (Błaszczyk and Snopek 2017, Theorem 4.6). In the discrete case, the following theorem can be proved by repeating the reasoning presented in Błaszczyk and Snopek (2017).
The function α i 1 ,...,i 4 has also been implemented on the basis of the qtfm package as a function alpha(o, i1, i2, i3, i4) (see Listing 6), where o is the transformed octonion (or matrix of octonions) and i1, …, i4 are four (different) indices of imaginary parts of o whose sign should be changed. The tests that could be performed thanks to this (presented in Listing 7) illustrate Theorem 10 and show that indeed DOFT of real-value matrices has the mentioned symmetry properties.

Discussion and conclusions
The results presented show that discrete Fourier transforms can be generalized to the case of higher order algebras (e.g. octonions). What's more, using the properties of algebra of quadruple-complex numbers, this generalization can lead to a very similar form.
One could ask the question why it is worth working in two different algebras in parallel instead of just working in the algebra of quadruple-complex numbers F. While going through calculations in F facilitates the interpretation of results, there are more and more papers in the literature devoted to the use of octonion algebra (Gao and Lam 2014;Grigoryan and Agaian 2018;Hahn and Snopek 2016;Klco et al. 2017;Lazendić et al. 2018a, b;Lian 2019;Popa 2016Popa , 2018Sheng et al. 2018;Snopek 2015;Wang et al. 2017). Therefore, it seems important to develop tools enabling work in this algebra as well.
The properties of discrete octonion Fourier transforms show that they can be used without difficulty for the analysis of difference equations, as well as for the analysis of finite difference schemes for differential equations. Detailed work in this area remains in the plans for further actions, as well as the further development of this theory.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.