Optimal Restoration of Multiple Signals in Quaternion Algebra Artyom - - PowerPoint PPT Presentation

Optimal Restoration of Multiple Signals in Quaternion Algebra

Artyom M. Grigoryana and Sos S. Agaianb

aDepartment of Electrical and Computer Engineering

The University of Texas at San Antonio, San Antonio, Texas, USA, and

bComputer Science Department, College of Staten Island and the Graduate

Center, Staten Island, NY, USA amgrigoryan@utsa.edu, sos.agaian@csi.cuny.edu

7/8/2020 1

OUTLINE

• Introduction
• Problem Formulation
• Quaternion Numbers
• 1-D Convolution in Quaternion Algebra (QA)
• Inverse Problem in QA and Solution
• Quaternion Convolution plus Noise: Solution
• Example of Restored Quaternion Signal
• Summary
• References

7/8/2020 2

Abstract

• This paper offers a new multiple signal restoration tool to solve the inverse

problem, when signals are convoluted with a multiple impulse response and then degraded by an additive noise signal with multiple components.

Inverse problems arise practically in all areas of science and engineering and refers to as methods of estimating data/parameters, in our case of multiple signals that cannot directly be observed.

• The presented tool is based on the mapping multiple signals into the

quaternion domain, and then solving the inverse problem.

• Due to the non-commutativity of quaternion arithmetic, it is difficult to find

the optimal filter in the frequency domain for degraded quaternion signals.

7/8/2020 3

Presented Work

• As an alternative, we introduce an optimal filter by using special 4×4

matrices on the discrete Fourier transforms of signal components, at each frequency-point. The optimality of the solution is with respect to the mean-square-root error, as in the classical theory of the signal restoration by the Wiener filter.

• The Illustrative example of optimal filtration of multiple degraded

signals in the quaternion domain is given. The computer simulations validate the effectiveness of the proposed method.

7/8/2020 4

7/8/2020 5

• 1. PROBLEM OF MULTIPLE SIGNAL RESTORATION

In the space of quaternion signals, in the model described the signal 𝑟(𝑢) convoluted with the function ℎ(𝑢) plus a noise 𝑜(𝑢) 𝑗(𝑢) = 𝑟(𝑢) ∗ ℎ(𝑢) + 𝑜(𝑢), (1) the signal 𝑟(𝑢) is restoring from the degraded signal 𝑗(𝑢). The classic case: the inverse problem is solving by the optimal filter 𝑍 𝜕 = 𝐼(𝜕) 𝐼 𝜕 + 𝜚𝑂/𝑅(𝜕) . (2) Here, 𝐼(𝜕) is the Fourier transform of ℎ(𝑢), and 𝜚𝑂/𝑅(𝜕) is the noise-signal ratio, and 𝜚𝑂 𝜕 =< 𝑂 𝜕

2 > and 𝜚𝑅 𝜕 =< 𝑅 𝜕 2 > are spatial spectral

densities of the signal 𝑟(𝑢) and noise 𝑜(𝑢), respectively.

7/8/2020 6

Inroduction to Quanterions

The quaternion number is composed by one real part and three-component imaginary part, 𝑟 = 𝑏 + (𝑐𝑗 + 𝑑𝑘 + 𝑒𝑙) = 𝑏 + 𝑐𝑗 + 𝑑𝑘 + 𝑒𝑙, where 𝑏, 𝑐, 𝑑, and 𝑒 are real numbers. Together with unit 1, three imaginary units 𝑗, 𝑘, and 𝑙 are used with the multiplication laws, which are following: 𝑗2 = 𝑘2 = 𝑙2 = −1, 𝑙𝑗 = −𝑗𝑙 = 𝑘, 𝑗𝑘 = −𝑘𝑗 = 𝑙, 𝑘𝑙 = −𝑙𝑘 = 𝑗. The quaternion conjugate and modulus of 𝑟 are defined as 𝑟 = 𝑏 − (𝑐𝑗 + 𝑑𝑘 + 𝑒𝑙) and 𝑟 = 𝑏2 + 𝑐2 + 𝑑2 + 𝑒2. The multiplication of quaternions is not a commutative operation, i.e., 𝑟1𝑟2 ≠ 𝑟2𝑟1 for many quaternions 𝑟2 ≠ 𝑟1.

7/8/2020 7

In the definition of the 𝑂-point quaternion DFT (QDFT), the exponential kernel is used the exponential kernel is used 𝑋

𝜈 = exp(−𝜈2𝜌/𝑂) = cos(2𝜌/𝑂) − 𝜈 sin(2𝜌/𝑂),

where 𝜈 is a pure unit quaternion, 𝜈 = 𝑛1𝑗 + 𝑛2𝑘 + 𝑛3𝑙. For a such number, |𝜈| = 1 and 𝜈2 = −1. The 𝑂-point right-side QDFT 𝑅𝑞 = ෍

𝑜=0 𝑂−1

𝑟𝑜𝑋

𝜈 𝑜𝑞 ,

𝑞 = 0: 𝑂 − 1 . The fast algorithms for the 𝑂-points QDFT exist [1]. Because of not commutativity of multiplication in quaternion arithmetic, the main operation of the cyclic convolution is not reduced to the multiplication of the QDFTs, as for the traditional 𝑂-point DFT.

7/8/2020 8

1-D CONVOLUTION IN QUATERNIOIN ALGEBRA

In this section, we show a new technique for calculating the quaternion convolution by the Fourier transforms, which was developed by Grigoryan in

• 2015. The convolution will be transformed into the frequency domain and a new
• peration of multiplication of transforms by 4×4 matrices will be performed,

instead of point-wise multiplication of the DFTs in the traditional method. The quaternion signal 𝑔

𝑜, 𝑜 = 0: (𝑂 − 1), of length 𝑂

𝑔

𝑜 = (𝑔 1)𝑜, (𝑔 2)𝑜, (𝑔 3)𝑜, (𝑔 4)𝑜 = (𝑔 1)𝑜+𝑗(𝑔 2)𝑜+𝑘(𝑔 3)𝑜+𝑙(𝑔 4)𝑜

can be considered in the following matrix representation.

7/8/2020 9

We introduce the following 4 matrices with the multiplications shown in the table

𝐹 = 1 1 1 1 , 𝐽 = −1 1 1 −1 , 𝐾 = −1 −1 1 1 , 𝐹 = −1 1 −1 1 , (4)

. The quaternion signal at the point 𝑜 can be written as a matrix, (we use the same notation 𝑔

𝑜),

𝑔

𝑜 = 𝑔 1 𝑜𝐹 + 𝑔 2 𝑜𝐽 + 𝑔 3 𝑜𝐾 + 𝑔 4 𝑜𝐿 =

𝑔

1 𝑜

− 𝑔

2 𝑜

− 𝑔

3 𝑜

− 𝑔

4 𝑜

𝑔

2 𝑜

𝑔

1 𝑜

𝑔

4 𝑜

− 𝑔

3 𝑜

𝑔

3 𝑜

− 𝑔

4 𝑜

𝑔

1 𝑜

𝑔

2 𝑜

𝑔

4 𝑜

𝑔

3 𝑜

− 𝑔

2 𝑜

𝑔

1 𝑜

, (5)

7/8/2020 10

Another quaternion sequence, which we call the impulse response characteristic

• f a linear system, is denoted by ℎ𝑜 = ((ℎ1)𝑜, (ℎ2)𝑜, (ℎ3)𝑜, (ℎ4)𝑜) and can be

presented in matrix form as ℎ𝑜 = ℎ1 𝑜𝐹 + ℎ2 𝑜𝐽 + ℎ3 𝑜𝐾 + ℎ4 𝑜𝐿. (6) We define the circular linear convolution as 𝑕𝑜 = 𝑔 ∗ ℎ 𝑜 = ෍

𝑛=0 𝑂−1

𝑔

𝑜−𝑛ℎ𝑛 ≜ ෍ 𝑛=0 𝑂−1

𝑔

(𝑜−𝑛)mod𝑂ℎ𝑛 .

(7) The following two sums are different: 𝑔 ∗ ℎ 𝑜 ≠ ℎ ∗ 𝑔 𝑜.

7/8/2020 11

To simplify calculations, we separate all 𝐹, 𝐽, 𝐾, 𝐿-components of the signals

𝑔

𝐹 =

𝑔

1 0, 𝑔 1 1, 𝑔 1 2, … , 𝑔 1 𝑂−1 ,

𝑔

𝐽 = (𝑔 2)0, (𝑔 2)1, (𝑔 2)2, … , (𝑔 2)𝑂−1 ,

𝑔

𝐾 = (𝑔 3)0, (𝑔 3)1, (𝑔 3)2, … , (𝑔 3)𝑂−1 ,

𝑔

𝐿 = (𝑔 4)0, (𝑔 4)1, (𝑔 4)2, … , (𝑔 4)𝑂−1 ,

ℎ𝐹 = ℎ1 0, ℎ1 1, ℎ1 2, … , ℎ1 𝑂−1 , ℎ𝐽 = (ℎ2)0, (ℎ2)1, (ℎ2)2, … , (ℎ2)𝑂−1 , ℎ𝐾 = (ℎ3)0,(ℎ3)1, (ℎ3)2, … , (ℎ3)𝑂−1 , ℎ𝐿 = (ℎ4)0, (ℎ4)1, (ℎ4)2, … , (ℎ4)𝑂−1 .

The 𝐹, 𝐽, 𝐾-, and 𝐿-components of the quaternion convolution 𝑕𝑜 are denoted by 𝑕𝐹, 𝑕𝐽, 𝑕𝐾, and 𝑕𝐿, respectively. By using the table of multiplications for the basic four quaternion matrices, we can open the convolution

𝑕𝑜 = 𝑔

1 𝑜𝐹 + 𝑔 2 𝑜𝐽 + 𝑔 3 𝑜𝐾 + 𝑔 4 𝑜𝐿 ∗

ℎ1 𝑜𝐹 + ℎ2 𝑜𝐽 + ℎ3 𝑜𝐾 + ℎ4 𝑜𝐿

and consider the following system of Eqs for the 𝐹, 𝐽, 𝐾,𝐿-components of 𝑕𝑜:

ቋ 𝑕𝐹 = 𝑔

𝐹 ∗ ℎ𝐹 − 𝑔 𝐽 ∗ ℎ𝐽 − 𝑔 𝐾 ∗ ℎ𝐾 − 𝑔 𝐿 ∗ ℎ𝐿,

𝑕𝐽 = 𝑔

𝐹 ∗ ℎ𝐽 + 𝑔 𝐽 ∗ ℎ𝐹 − 𝑔 𝐾 ∗ ℎ𝐿 + 𝑔 𝐿 ∗ ℎ𝐾,

𝑕𝐾 = 𝑔

𝐹 ∗ ℎ𝐾 + 𝑔 𝐽 ∗ ℎ𝐿 + 𝑔 𝐾 ∗ ℎ𝐹 − 𝑔 𝐿 ∗ ℎ𝐽,

𝑕𝐿 = 𝑔

𝐹 ∗ ℎ𝐿 − 𝑔 𝐽 ∗ ℎ𝐾 + 𝑔 𝐾 ∗ ℎ𝐽 + 𝑔 𝐿 ∗ ℎ𝐹.

7/8/2020 12

For each real component of the signal 𝑔

𝑜 and the impulse response sequence ℎ𝑜,

we will use the corresponding capital letters for the corresponding DFTs. In the frequency domain, the system of Eq. 8 for frequency-point 𝑞 can be written as

𝐻𝐹 𝐻𝐽 𝐻𝐾 𝐻𝐿 = 𝐼𝐹 −𝐼𝐽 −𝐼𝐾 −𝐼𝐿 𝐼𝐽 𝐼𝐹 −𝐼𝐿 𝐼𝐾 𝐼𝐾 𝐼𝐿 𝐼𝐹 −𝐼𝐽 𝐼𝐿 −𝐼𝐾 𝐼𝐽 𝐼𝐹 𝐺

𝐹

𝐺

𝐽

𝐺

𝐾

𝐺

𝐿

.

This is a compact form of the equation 𝐻 = 𝑰𝐺 at one frequency-point 𝑞, which is

• mitted from the notation. For this point, the matrix 𝑰 is the 4×4 matrix

𝑰 = 𝐼𝐹 −𝐼𝐽 −𝐼𝐾 −𝐼𝐿 𝐼𝐽 𝐼𝐹 −𝐼𝐿 𝐼𝐾 𝐼𝐾 𝐼𝐿 𝐼𝐹 −𝐼𝐽 𝐼𝐿 −𝐼𝐾 𝐼𝐽 𝐼𝐹 .

This matrix is orthogonal.

7/8/2020 13

The convoluted quaternion signal

𝑕𝑜 = 𝑕1 𝑜𝐹 + 𝑕2 𝑜𝐽 + 𝑕3 𝑜𝐾 + 𝑕4 𝑜𝐿

can be calculated as

𝑕1 = 𝑮−1𝐻𝐹, 𝑕2= 𝑮−1𝐻𝐽, 𝑕3= 𝑮−1𝐻𝐾, 𝑕4= 𝑮−1𝐻𝐿.

Here, 𝑮−1 is the matrix of the inverse 𝑂-point DFT. Together with the unite matrix 𝐹, we consider the following three new quaternion matrices:

𝐽∗ = −1 1 −1 1 , 𝐾∗ = −1 1 1 −1 , 𝐿∗= −1 −1 1 1 .

7/8/2020 14

Statement 1. Given a frequency-point 𝑞, the inverse matrix 𝑰−1 at this point is 𝑰−1 = 1 det 𝑰 𝑰𝑈 = 1 det 𝑰 𝐹𝐼𝐹 − 𝐽∗𝐼𝐽 − 𝐾∗𝐼𝐾 − 𝐿∗𝐼𝐿 , (13) where the square root of the determinant of the matrix is calculated by det 𝑰 = 𝐼𝐹

2 + 𝐼𝐽 2 +

𝐼𝐾

2 + 𝐼𝐿

• 2. The operation 𝑰𝑈 denotes the transposition of the complex matrix 𝑰.

The reconstruction of the frequency-vector 𝐺 at the frequency-point 𝑞 is the inverse transform 𝐺

𝐹

𝐺

𝐽

𝐺

𝐾

𝐺

𝐿

= 1 det 𝑰 𝐼𝐹 𝐼𝐽 𝐼𝐾 𝐼𝐿 −𝐼𝐽 𝐼𝐹 𝐼𝐿 −𝐼𝐾 −𝐼𝐾 −𝐼𝐿 𝐼𝐹 𝐼𝐽 −𝐼𝐿 𝐼𝐾 −𝐼𝐽 𝐼𝐹 𝐻𝐹 𝐻𝐽 𝐻𝐾 𝐻𝐿 . (14) The original quaternion signal 𝑔

𝑜 = 𝑔 1 𝑜𝐹 + 𝑔 2 𝑜𝐽 + 𝑔 3 𝑜𝐾 + 𝑔 4 𝑜𝐿 is reconstructed from the

components of the inverse transform by 𝑔

1 = 𝑮−1𝐺 𝐹,

𝑔

2 = 𝑮−1𝐺 𝐽, 𝑔 3= 𝑮−1𝐺 𝐾, 𝑔 4= 𝑮−1𝐺 𝐿.

(15)

7/8/2020 15

Algorithm 1. Inverse problem: 1-D Quaternion convolution

• 1. Input signal 𝑔

𝑜 = ((𝑔 1)𝑜,(𝑔 2)𝑜, (𝑔 3)𝑜,(𝑔 4)𝑜), and the impulse response ℎ𝑜 = ((ℎ1)𝑜, (ℎ2)𝑜, (ℎ3)𝑜,(ℎ4)𝑜).

• 2. Calculate four 1-D DFTs 𝐺𝐹, 𝐺𝐽, 𝐺

𝐾, and 𝐺𝐿 of (𝑔 1)𝑜, (𝑔 2)𝑜, (𝑔 3)𝑜, and (𝑔 4)𝑜, respectively.

• 3. Calculate 1-D DFTs 𝐼𝐹, 𝐼𝐽, 𝐼

𝐾, and 𝐼𝐿 of (ℎ1)𝑜, (ℎ2)𝑜, (ℎ3)𝑜, and (ℎ4)𝑜, respectively.

• 4. Compose the convoluted quaternion signal 𝑕𝑜 from the inverse 𝑂-point DFTs by Eq. 11.

a. For each frequency-point 𝑞 ∈ {0,1,2, . . . , (𝑂 − 1)}, calculate the 4×4 matrix 𝑰 by Eq. 10. b. Calculate the data 𝐻𝐹, 𝐻𝐽, 𝐻𝐾, and 𝐻𝐿 by Eq. 9. c. Calculate four inverse 𝑂-point DFTs of 𝐻𝐹, 𝐻𝐽, 𝐻𝐾, and 𝐻𝐿.

• 5. Compose the quaternion signal 𝑔

𝑜 from the inverse 𝑂-point DFTs by Eq. 15.

a. Calculate the inverse 4×4 matrix 𝑰−1 by Eq. 13. b. Apply the inverse matrix on the vector (𝐻𝐹, 𝐻𝐽, 𝐻𝐾, 𝐻𝐿). c. Calculate four inverse 𝑂-point DFTs of 𝐺𝐹, 𝐺𝐽, 𝐺

𝐾, and 𝐺𝐿.

7/8/2020 16

Figure 2: Four components of the original quaternion signal 𝑔

𝑜.

7/8/2020 17

Figure 3: Four components of the convoluted quaternion signal 𝑕𝑜.

7/8/2020 18

We consider the corresponding binary matrix 𝑩 of signs of elements of 𝑰, 𝑩 = 1 −1 −1 −1 1 1 −1 1 1 1 1 −1 1 −1 1 1 , 𝑳 = −2 2 2 2 . (16) Here, the vector 𝑳 is the sum of elements of the matrix 𝑩 by rows. Therefore, the components of the convolution are plotted after the following amplitude transformation: 𝑕1 𝑜, 𝑕2 𝑜, 𝑕3 𝑜, 𝑕4 𝑜 → − 𝑕1 𝑜/2, 𝑕2 𝑜/2, 𝑕3 𝑜/2, 𝑕4 𝑜/2 The means-square-root error therefore is calculated as 𝑓 𝑔, 𝑕 = 1 𝑂 ෍

𝑜=0 𝑂−1

𝑔

1 𝑜 + 1

2 𝑕1 𝑜

2

+ ෍

𝑙=2 4

𝑔

𝑙 𝑜 − 1

2 𝑕𝑙 𝑜

2

. (17)

7/8/2020 19

Quaternion Convolution and Noise

The frequency characteristics of the input and output signals are considered in forms of quaternion matrices in the (𝐹, 𝐽∗, 𝐾∗, 𝐿∗)-basis, i.e., at each frequency-points 𝑞, we consider that 𝐺 = 𝐹𝐺

𝐹 + 𝐽∗𝐺 𝐽 + 𝐾∗𝐺 𝐾 + 𝐿∗𝐺 𝐿 and 𝐻 = 𝐹𝐻𝐹 + 𝐽∗𝐻𝐽 + 𝐾∗𝐻𝐾 + 𝐿∗𝐻𝐿.

It is important to note and not difficult to verify that equation (9) can be written as

𝐻𝐹 −𝐻𝐽 −𝐻𝐾 −𝐻𝐿 𝐻𝐽 𝐻𝐹 −𝐻𝐿 𝐻𝐾 𝐻𝐾 𝐻𝐿 𝐻𝐹 −𝐻𝐽 𝐻𝐿 −𝐻𝐾 𝐻𝐽 𝐻𝐹 = 𝐼𝐹 −𝐼𝐽 −𝐼𝐾 −𝐼𝐿 𝐼𝐽 𝐼𝐹 −𝐼𝐿 𝐼𝐾 𝐼𝐾 𝐼𝐿 𝐼𝐹 −𝐼𝐽 𝐼𝐿 −𝐼𝐾 𝐼𝐽 𝐼𝐹 𝐺

𝐹

−𝐺

𝐽

−𝐺

𝐾

−𝐺

𝐿

𝐺

𝐽

𝐺

𝐹

−𝐺

𝐿

𝐺

𝐾

𝐺

𝐾

𝐺

𝐿

𝐺

𝐹

−𝐺

𝐽

𝐺

𝐿

−𝐺

𝐾

𝐺

𝐽

𝐺

𝐹

.

Thus, we obtain the equation 𝐻 = 𝑰𝐺, which means that at the frequency-point 𝑞, 𝐻𝑞 = 𝑰𝑞𝐺

𝑞

where all components of this equation are 4×4 matrices. We consider the model with the noise, when 𝑕𝑜 = 𝑔

𝑜 ∗ ℎ𝑜 + 𝑜𝑜 and 𝑜𝑜 is a noise. The

transformation of this equation into the frequency-domain results in the similar equation 𝐻𝑞 = 𝑰𝑞𝐺

𝑞 + 𝑂𝑞,

𝑞 = 0: (𝑂 − 1). (18) Here 𝑂 = 𝐹𝑂𝐹 + 𝐽∗𝑂𝐽 + 𝐾∗𝑂

𝐾 + 𝐿∗𝑂𝐿 and 𝑂 = (𝑂)𝑞.

7/8/2020 20

Consider a quaternion signal 𝑃 = 𝐹𝑃𝐹 + 𝐽∗𝑃𝐽 + 𝐾∗𝑃𝐾 + 𝐿∗𝑃𝐿 which is a filtered signal 𝑃 = 𝑍𝐻, i.e., at the frequency 𝑞, we have 𝑃𝑞 = 𝑍

𝑞𝐻𝑞 = 𝑍 𝑞(𝑰𝑞𝐺 𝑞 + 𝑂𝑞). The filter, or the frequency

characteristics Y is defined as the filter that minimizes the RMS error between 𝑃 and 𝐺. To derive such an optimal filter, we consider the Lagrangian ℒ 𝑍

𝑞 =

𝐺

𝑞 − 𝑍 𝑞(𝑰𝑞𝐺 𝑞 + 𝑂𝑞) 2 = min.

Therefore, the optimal filter is 𝑍 = 𝑰′𝜚𝐺 𝑰 2𝜚𝐺 + 𝜚𝑂 = ഥ 𝑰𝑈𝜚𝐺 𝑰 2𝜚𝐺 + 𝜚𝑂 = ഥ 𝑰𝑈 𝑰 2 + 𝜚𝑂/𝐺 . (21) Here, the noise to signal ratio 𝜚𝑂/𝐺 = 𝜚𝑂/𝜚𝐺 . In the given (𝐹, 𝐽∗, 𝐾∗, 𝐿∗)-basis, the matrix- filter 𝑍 = 𝑍

𝑞 can be written as

𝑍 = 𝜚𝐺 𝑰 2𝜚𝐺 + 𝜚𝑂 ഥ 𝑰𝐹 − 𝐽∗ഥ 𝑰𝐽 − 𝐾∗ഥ 𝑰𝐾 − 𝐿∗ഥ 𝑰𝐿 . (22)

7/8/2020 21

Thus, the inverse problem of restoration has been solved. The result of

• ptimal filtration is 𝑃 = 𝑍𝐻 at the frequency-point 𝑞. The filtered quaternion

signal መ 𝑔 = ( መ 𝑔

𝐹, መ

𝑔

𝐽, መ

𝑔

𝐾, መ

𝑔

𝐿) is calculated as

መ 𝑔

𝐹 = 𝑮−1𝑃𝐹,

መ 𝑔

𝐽 = 𝑮−1𝑃𝐽,

መ 𝑔

𝐾 = 𝑮−1𝑃𝐾,

መ 𝑔

𝐿 = 𝑮−1𝑃𝐿.

The optimality is with respect to the MSR error. In the case when there is no noise, 𝜚𝑂 = 0, we obtain the inverse quaternion filter described above.

7/8/2020 22

Figure 4: The 𝑓-, 𝑗-, 𝑘-, and 𝑙 components of the noisy quaternion signal.

7/8/2020 23

Figure 5: The 𝑓-, 𝑗-, 𝑘-, and 𝑙 components of the filtered quaternion signal 𝑝𝑜.

7/8/2020 24

Figure 6: The 𝑓-, 𝑗-, 𝑘-, and 𝑙 components of the filtered quaternion signal 𝑝𝑜 are plotted together with the corresponding components of the original signal 𝑔

𝑜.

7/8/2020 25

SUMMARY

• The multiple signals were considered in the quaternion space and the cyclic quaternion

convolution was described and calculated in the frequency domain, by using the right-side quaternion discrete Fourier transform (QDFT).

• The traditional approach of multiplication QDFTs does not work in quaternion algebra,

because of non-commutativity of quaternion multiplication.

• The components of the convolution were processed in the frequency domain by special 4×4

matrices that allows us to solve the inverse problem, namely, to restore the degraded quaternion signal in a linear system with the convolution plus a noise.

• The characteristic of the quaternion optimal filter was found.

REFERENCES 1. A.M. Grigoryan, S.S. Agaian, Quaternion and Octonion Color Image Processing with MATLAB, SPIE PRESS, vol. PM279 (2018), (https://doi.org/10.1117/3.2278810). 2. A.M. Grigoryan, E.R. Dougherty, “Bayesian robust optimal linear filters,” Signal Processing, vol. 81, no. 12, pp. 2503–2521 (2001). 3. A.M. Grigoryan, E.R. Dougherty, “Robustness of optimal binary filters,” Journal of Electronic Imaging, vol. 7, no. 1, pp. 117–126 (1998), (https://doi.org/10.1117/1.482633)