Sparse fast Clifford Fourier transform

(1)

Sparse fast Clifford Fourier transform

^*

Rui WANG¹, Yi-xuan ZHOU¹, Yan-liang JIN¹, Wen-ming CAO^†‡2,3

(¹School of Communication and Information Engineering, Shanghai University, Shanghai 200444, China) (²College of Information Engineering, Shenzhen University, Shenzhen 518060, China) (³Department of Electrical and Computer Engineering, University of Missouri, Columbia 65211, USA)

†E-mail: [email protected]

Received Dec. 9, 2015; Revision accepted Aug. 23, 2016; Crosschecked Aug. 24, 2017

Abstract: The Clifford Fourier transform (CFT) can be applied to both vector and scalar fields. However, due to problems with big data, CFT is not efficient, because the algorithm is calculated in each semaphore. The sparse fast Fourier transform (sFFT) theory deals with the big data problem by using input data selectively. This has inspired us to create a new algorithm called sparse fast CFT (SFCFT), which can greatly improve the computing performance in scalar and vector fields. The experiments are implemented using the scalar field and grayscale and color images, and the results are compared with those using FFT, CFT, and sFFT. The results demonstrate that SFCFT can effectively improve the performance of multivector signal processing.

Key words: Sparse fast Fourier transform (sFFT); Clifford Fourier transform (CFT); Sparse fast Clifford Fourier transform (SFCFT); Clifford algebra

http://dx.doi.org/10.1631/FITEE.1500452 CLC number: TP391

1 Introduction

Methods have been used to analyze and visualize unstructured vector field data. There are basically two approaches: extracting features to determine the dataset and visualizing an entire dataset. As dataset sizes increase, feature extraction becomes increasingly important. A method called ‘Clifford Fourier transform (CFT)’ was proposed to apply the Fourier transform to vector fields (Ebling and Scheuermann, 2005; Schlemmer et al., 2005; Hitzer and Mawardi, 2008). Hitzer and Sangwine (2013) provided an overview of the modern development of CFT and wavelet transformations. In recent years, a detailed algebraic characterization of the continuous mani- folds of square roots of í1 has been established in all real Clifford algebras (Hitzer et al., 2013). This is

why CFT uses multivector square roots of –1 instead of the complex imaginary unit. Based on this, a general CFT was introduced (Hitzer, 2013). CFT has already become one of the most significant steps in image processing and analysis (Cao and Feng, 2010;

Xu et al., 2011; Cao et al., 2013; Wang et al., 2013), network analysis (Li, 2005; Xie et al., 2008; Wang et al., 2015), color edge detection (Evans et al., 2000;

Mishra et al., 2015), and image registration (Batard et al., 2010). CFT allows for frequency analysis of vector fields and the behavior of vector-valued filters.

In the frequency space, CFT can transform vectors into general Clifford algebra multivectors. Many basic vector-valued patterns (such as source, sink, saddle points, and potential vortices) can be described by a certain number of multivectors in the frequency space (Ebling and Scheuermann, 2005).

On the other hand, explosive growth in data makes big data problems a cause for concern. Some solutions have already been proposed, such as using fast approximate correlation (Mueen et al., 2010), sparse matrices (Gilbert and Indyk, 2010), and the sparse recovery system (Porat and Strauss, 2012).

Frontiers of Information Technology & Electronic Engineering www.jzus.zju.edu.cn; engineering.cae.cn; www.springerlink.com ISSN 2095-9184 (print); ISSN 2095-9230 (online)

E-mail: [email protected]

‡ Corresponding author

* Project supported by the National Natural Science Foundation of China (Nos. 61301027, 61375015, and 11274226)

ORCID: Yi-xuan ZHOU, http://orcid.org/0000-0002-5008-5802

© Zhejiang University and Springer-Verlag GmbH Germany 2017

(2)

Iwen (2010) presented a method to estimate the Fou- rier representation for sparse signals. However, Has- sanieh et al. (2012b) developed a new approach called

‘sparse fast Fourier transform (sFFT)’ to process big data, especially sparse data. sFFT deals with big data problems by using only a small subset of the input data to compute a compressed Fourier transform. As a result, sFFT is much faster than the fast Fourier transform (FFT). sFFT has already been applied to areas such as GHz-wide sensing and decoding (Has- sanieh et al., 2014), light field reconstruction (Shi et al., 2014), and the Global Positioning System (GPS) (Hassanieh et al., 2012a). Our work is to combine sFFT with CFT to operate on multivector signals.

sFFT can choose ‘large’ coefficients for calculation.

Combining it with CFT reduces useless or remote data, and thus the Fourier transform of multivector signals can be computed more efficiently (Schu- macher and Puschel, 2014).

2 Related work 2.1 Clifford algebra

Clifford algebra provides a powerful computing framework without using coordinate information.

Hestenes and Sobczyk (1984) used it as an efficient and versatile computational tool. It integrates vector algebra, matrix algebra, and complex numbers algebra into a coherent mathematical language which not only reduces the computation complexity, but also improves the computation efficiency (Hestenes, 1999).

Let Ü^d be a d-dimensional Euclidean vector space and G_d a 2^d-dimensional real Clifford algebra over Ü. One obtains G3 by using the rules of 3D Clifford algebra with the given basis {1, e₁, e₂, e₃, e₁e₂, e₁e₃, e₂e₃, e₁e₂e₃}:

, 1,2,3,

j j j

1e e (1) , 1,2,3,

j j 1 j

e e (2)

, 1,2,3, .

j k k jj k j kz

e e e e (3)

Note that the multiplication of Clifford algebra is not commutative.

We define an arbitrary multivector P as

i (3 ),

P Į a b E (4) with Į, ȕÜ, a, bÜ³, i₃=e₁e₂e₃, and (i₃)²=í1. We call such components, Į, a, i₃b, and i₃ȕ, as blades of Clifford algebra that are generalizations of the con- cept of scalars and vectors.

We define the grade projector ¢²_j: G₃ĺG₃ as

0 , 1 ,

¢ ² P Į ¢ ² P a (5)

2 i ,3 3 i .3

¢ ² P b ¢ ² P E (6) The Clifford product of two vectors a, bÜ³ results in

, ,

¢ ²

ab a b a b (7) where ,¢ ² is the inner product and is the outer product.

As shown in Eq. (4), a multivector is a linear combination of components in different grades. Thus, by using linearity, the Clifford product of multivectors involves Clifford products of component blades.

The appendix in Ebling and Scheuermann (2003) contains an example of multiplying multivectors in 2D, and Table 2 on page 344 in Hitzer (2012) shows the multiplication table for G₃.

Further, we have

0 , cos ,Z

¢ ² ¢ab a b² a b (8)

2 sin ,Z

¢ ² ab a b a b (9)

where the usual norm for vectors, || ab ₂||

2

( 2) ,

 ab and Ȧ is the angle between a and b.

2.2 Clifford convolution

The convolution theorem states that a convolution in the time domain is equal to a multiplication in the frequency space. Therefore, convolution becomes an important tool in signal processing. Bujack et al.

(2015) presented a convolution theorem for general CFTs with separable mappings being orthogonal, separable, and linear. It can be interpreted as a mutual commutation or anticommutation among the functions.

(3)

Let U, V: Ü^dĺG_d (d=2, 3) be two multivector fields. Then the Clifford convolution can be defined as

(

(V U )( )x

³

\_dV [) (U x[) | d |,[ (10) where xV and ȗU.

It can be seen that the Clifford convolution is a conventional convolution when both fields are scalar.

Further, let U be a scalar field and V a vector field.

The Clifford convolution then allows a scalar multiplication and vector smoothing model.

If both V and U are vector fields, we obtain the following relationship:

( ) ( ) ( ),( ) ( ) ( ).

V [ U x [ V [ x [ V [ U x [

(11) It means that the convolution in the multivector field contains additional information, as the scalar part of the convolution result is a conventional convolution.

Thus, Clifford convolution can be used to better analyze the vector field data.

2.3 Clifford Fourier transform

CFT allows a transformation from the position space to the frequency space, allowing signals to be analyzed in the frequency space. This makes it feasi- ble to use such features as phases and frequencies of the signals. Further, the convolution theorem allows better filter response when a signal is analyzed in the frequency space (Ebling and Scheuermann, 2005).

Similarly, Clifford convolution can be transformed into the frequency space. Thus, CFT can be developed as an extension of the usual Fourier transform for vector fields (Schlemmer et al., 2005; Reich and Scheuermann, 2010).

2.3.1 Clifford Fourier transform in the 2D space Ebling and Scheuermann (2005) proposed the CFT on both scalar and multi-dimensional fields, initiating a novel approach to analyzing vector fields.

Subsequently, the generalization of the Fourier transform in Clifford geometric algebra has been presented (Hitzer and Mawardi, 2008; Hitzer and Sangwine, 2013).

CFT for multivector-valued functions F: Ü²oG₂ can be defined as

2

{ }( )F

³

F( )exp( 2 i S ¢2 , )d ,²

\

u x x u x

F ⁽¹²⁾

where vectors x, uÜ².

Based on Clifford algebra, a 4D multivector field can be written as

0 1 1 2 2 12 12

0 1 1 2 1 2 12 2

0 12 2 1 1 2 2

( ) ( ) ( ) ( ) ( )

= ( ) ( ) ( ) i ( )i ( ( ) ( )i ) ( ( ) ( )i ).

F F F F F

F F F F

1

x x x e x e x e

x x e x e x

x x e x x

(13) Then, considering the linearity of CFT, one obtains

0 12 2

1 1 2 2

{ }( ) ( { ( ) ( )i }( )) ( { ( ) ( )i }( )).

F F F

F F

1

u x x u

e x x u

F F

F (14)

It can be seen that 2D CFT is the linear combi- nation of two traditional complex Fourier transforms.

2.3.2 Clifford Fourier transform in the 3D space CFT for multivector-valued functions F: Ü³oG₃ can be defined as

3

{ }( )F

³

F( )exp( 2 i S ¢3 , )d ,²

\

u x x u x

F ⁽¹⁵⁾

where vectors x, uÜ³.

Based on Clifford algebra, an 8D multivector field can be written as

0 1 1 2 2 3 3 23 23

31 31 12 12 123 123

0 1 1 2 2 3 3 23 3 1

31 3 2 12 3 3 123 3

i

i i i ,

F F F F F F

F F F

F F F F F

F F F

e e e e

e e e

e e e e

e e

(16)

which can be seen as the decomposition of an 8D multivector field into four complex signals. Each complex signal is transformed separately by the standard Fourier transform. Then, considering the linearity of CFT, one obtains

(4)

0 123 3

1 23 3 1

2 31 3 2

3 12 3 3

{ }( ) [ { ( ) ( )i }( )]

[ { ( ) ( )i }( )]

[ { ( ) ( )i }( )] .

F F F

F F

1

u x x u

x x u e

F F

F F F

(17)

This is useful as one can divide multivector algebra into four complex components. Also, the usual Fourier transform can be calculated separately for each direction {1, e1, e2, e3}.

3 Sparse fast Clifford Fourier transform 3.1 Sparse fast Fourier transform

Algorithms for computing the Fourier transform are inefficient because they take time proportional to the output size, and most of the Fourier signal coefficients are small or equal to zero in many applications. For instance, a normal 16×16 block in a video frame is sparse, since it has 228 negligible coefficients on average. It means that nearly 90% of the data are useless. Audio and image data are sparse as well. Compression schemes, such as JPEG and MPEG, have been provided to use this sparsity. The sFFT algorithm uses the sparsity of the signal spectrum and computes an approximated or compressed version of the Fourier transform. The algorithm uses only a small subset of the input data and runs in time proportional to the sparsity or desirable compression. It can be faster than that in time proportional to the signal length (Hassanieh et al., 2012c).

There are three stages of sFFT: identification of frequencies with large Fourier coefficients, accurate estimation of the Fourier coefficients of the frequencies identified in the first stage, and subtraction of the contribution of the partial Fourier representation computed in the first two stages (Gilbert et al., 2014).

Generally, repetitions of the three stages are guaranteed to concentrate a substantial portion of the energy with a high probability. In the first stage, frequencies are sampled randomly to obtain permuted Fourier coefficients. Then filters are used to divide the permuted coefficients into various frequency bands. In the final stage, the energy is estimated in each frequency band.

For the complex affine N-space ^^N, the Fourier component fˆ^^N can be defined as

1

0

ˆ 1 ^N _jexp( 2ʌi / ),

j

f f j N

Z N

¦

Z (18) where 0Ȧ and j<N.

Equivalently, one can simply obtain

1 0

ˆ exp ʌi(2 / ).

N

fj f_Z j N

Z

¦

(19)

Two basic Fourier properties are used for randomly permuting f ^^N and then fˆ. One is the scaling property, stating that for Įj=fcj, we have

ˆ 1

ˆ_j _{c j}

a f (where c^í1 is the inverse of (c mod N)); the other is the modulation property, stating that for aj=exp(2ʌibj/N)āfj, we have aˆ_j fˆ_{j b} . Then we ran- domly select two integers b, c[0, N] and define

^N

a as

2

exp( ʌi / ) .

j cj

a bj N f (20)

Note that for j[0, Ní1], ˆa is a permutation of fˆ , as entry ˆf_Z appears in ˆa as (Ȧc+b) mod N.

In the simplest case, the Fourier coefficients ˆf_Z can be computed for each Ȧ identified in the first stage by using L independent and uniformly distrib- uted random samples fl (L should be far smaller than N) and the estimator:

1

1 exp(

ˆ 2ʌi / ),

L l l

f f l N

Zc L

¦

Z (21) It can be seen that ˆf_Zc is an unbiased estimator for fˆ ._Z Therefore, high precision is guaranteed, as the estimator can approximate ˆf_Z with a high probability.

Assume that

^

^f^ˆ^Z^c^m ^m ^{1,2, ,}^" ^N

`

^{^}^, (22) which is the approximated sparse Fourier transform computed in the previous stage of the current repeti- tion. Here, Ȧ1, Ȧ2, …, ȦN are the frequencies that

(5)

have been identified in the first stage, while

1, 2,..

ˆ ˆ ., ˆ _N

f_Zc f_Zc f_Zc are the estimators whose Fourier coefficients are discovered in the second stage. In future iterations of all the stages, each sampled entry of f and fj can be replaced with

1

ˆ_mexp(2ʌi / ).

N

j m

m

f

¦

f_Z^c Z j N (23) In most cases, the entries of f used in each itera- tion can be predetermined, and then be used to update themselves all. These ‘updated samples’ are used in subsequent repetitions of the three stages (Gilbert et al., 2014).

3.2 Sparse fast Clifford Fourier transform Assume that a signal f has k nonzero coefficients.

We can regard the signal as k-sparse. More often, we calculate a signal’s k largest coefficients; thus, we call it ‘approximately k-sparse’. If its CFT is fˆ , then an approximation ˆf to fˆ ,c of an algorithm is required to be output to meet the following A2/A2 guarantee:

-sparse

2 2

ˆ ˆ min ˆ ,

Ck

c d

f f y f y (24) where y is a k-sparse standard reference. Note that the minimization is over k-sparse signals, and C is an approximation factor. Further, if k is small enough, the output of the sparse fast Clifford Fourier transform (SFCFT) can be expressed succinctly. The key to SFCFT is to filter useless and negligible coefficients, and then implement CFT.

Fig. 1 shows an example of normalizing a k-sparse vector field. Assume that a vector field (Fig. 1b) is sampled over a given cube (Fig. 1a). For each face of the cube, there is a vector ai (i=0, 1, …) on each sample point. Place the tail of each sample vector at the origin. Then normalize these vectors to form trivectors a_i/ 6a_i₁a_i₂. Sum the trivectors over all six faces of the cube (Fig. 1c). Note that such a trivector is not a triple of vectors, but a volume in space. The magnitude of the resulting trivector will be the volume of the sphere approximately. Fig. 1d shows the estimated k-sparse transformed vector field with an equal result. The normalization is carried out using the following steps:

1. Window function

Definition 1 (Window function) For an n- dimensional signal, a constant ׫>0, and a parameter į>0, we define the window function w(׫, į, Ȧ) to be a symmetric vector FÜⁿ with Ȧ[íʌ, ʌ]. The set of nonzero coordinates of vector F is nzc(F)[íȦ/2, Ȧ/2], and Fˆ₀ 1, and 0,Fˆ_i ! for i[í׫n, ׫n].

For each ׫ and į, a standard window function (׫, į, O(log(1/į)/׫)) exists (Smith, 2011), since one can obtain a standard window function through taking a Gaussian with standard deviation ⁴

^log(^ͳ^/^G⁾

and truncating it at Ȧ=O(log(1/į)/׫).

Definition 1 means that the window function w works like a filter, allowing us to concentrate on a subset of the Clifford Fourier coefficients.

2. Spectral permutation: This step allows us to permute the Clifford Fourier spectrum by permuting the time domain as follows:

Lemma 1 Given a constant vector aÜⁿ, a random integer ı, and an integer Ĳ[n], we define a transform Tı,Ĳ: (Tı,Ĳa)i=aıi+Ĳ, which refers to the permutation of the components of a. Then, n

,

, ,

(T_{V W}a)_V_i aˆ_{i j}Z^Wⁱ(T_{V W}a). Proof m, we have

(a) (b)

(c) (d)

Fig. 1 Example of normalizing a k-sparse vector field:

(a) the given cube; (b) normalized vector field; (c) trivectors over all six faces of the cube; (d) estimated k-sparse transformed vector field with an equal result

(6)

n ¹

1 1

( ) ,

1 1

( )

ˆ .

n n

mj m j

j j

i i j j

m m

T a a

a

V W V V V W W V

W V V

Z Z

Z

¦ ¦

a

Varying the permutation can vary the set of coefficients binned to a bucket. Note that for simplicity, our algorithm is analyzed on the condition that n is a power of two.

3. Subsampled CFT: suppose a vector a ^ⁿ and a parameter B dividing n, one obtains

( / )

ˆ_i ˆ_{i n B}.

y a (25)

Lemma 2 ˆy is the B-dimensional CFT of y_i

/ 1

0 .

n B i ai

¦

Proof _{( / )} ¹ ^{( / )} ^{1 /} ¹ ⁽ ^{) /}

0 0 0

ˆ_{i n B} ⁿ _j ^{ij n B} ^B ^{n B} _{Bj m} ^{i Bj m n B}

j t j

a

¦

aZ

¦ ¦

a Z

1 / 1 1

/ /

0 0 0

B n B B ˆ .

imn B imn B

Bj m i

m j m

a Z yZ y

¦ ¦

¦

^a

4. Location loops

Definition 2 (Hash function) Given a parameter B that divides n, we define a hash function hı: [n]ĺ[B]

by hı(i)=round(ıiB/n) and an ‘offset’ oı: [n]ĺ [ín/(2B), n/(2B)] by oı(i)=ıiíhı(i)(n/B), where [n]=

{0, 1, …, ní1}.

The location loops begin with a given parameter d and then output a parameter J that controls to find the coordinate of each large coefficient, and a set I[n]

of coordinates that contains each large coefficient.

5. Estimation loops: given a set I[n], we define G as a window function: w(1/B, į, Ȧ) and y=Gā(Tı,Ĳa);

thus, for iI, estimate ˆa as

ˆ

ˆ ˆ ,

i i

i

i h o

a z G

V V

ZW

c

where zˆ_i yˆ_{i n B}_{( / )}. This estimates each coordinate with a good probability.

Fig. 2 shows an example of loops of the algorithm when the input is sparse. Note that for simplicity, the process is analyzed only in the scalar field. The algorithm outputs Tı,Ĳa and obtains the spectrum permutation by permuting the Clifford algebra samples in the time domain (Fig. 2a). The algorithm computes signal f=GāTı,Ĳa in the time domain

(Fig. 2b). It can be seen that the spectrum of f is large around the large coordinates. The algorithm computes the subsampling of ˆf (Fig. 2c). Finally, the algorithm outputs the coordinates of each large coefficient (Fig. 2d). Note that there is a ‘hash collision’ when two coordinates are too close in the permutation. It results in missing the second coordinate from the left one.

Fig. 2 Results of estimating large regions in the scalar field: (a) permuted signal; (b) convolved signal; (c) actu- ally computed samples; (d) estimated large regions

(7)

3.3 The proposed algorithm

We use Clifford algebra to divide a multivector into a Clifford basis {1, e1, e2, e3, ...} and transform a Clifford basis into complex signals, as in Eqs. (13) and (16). Then we calculate the Fourier transform of each complex signal separately. The Fourier transform is usually dominated by a small number of

‘peaks’, meaning that it is sparse. The primary algorithm works by first ‘locating’ a set of loops that contain most of the peaks and then ‘estimating’ ˆf to obtain ˆ .z Thus, we can calculate the k-sparse Fourier transform as mentioned before. With the property of linearity, we add up each transform to obtain the final transform. The proposed method is provided in Algorithm 1.

4 Experiment analysis

In this section, we compare the SFCFT algorithm to several existing FFT implementations and demonstrate how SFCFT provides a better computational performance when processing multivector signals. We evaluate SFCFT performance in comparison to the fastest Fourier transform in the West (FFTW) (Frigo and Johnson, 2005), sFFT1.0, sFFT2.0 (Hassanieh et al., 2012c), and FCFT (Schlemmer et al., 2005) in the scalar field. FFTW is the fastest public implementation for computing FFT.

sFFT is the implementation for computing sFFT.

Compared with sFFT1.0, sFFT2.0 adds a heuristic to improve the runtime. FCFT is also a fast algorithm for computing CFT. We fix the signal size and signal sparsity, and present the runtime and robustness of the compared algorithms in Figs. 3–5. In higher- dimensional fields, we implement the algorithms on processing multispectral images and present the runtime and spectrum map in Figs. 6 and 7. For ac- curacy, each experiment is run more than 5000 times.

Each experiment is implemented on a 2.40 GHz Core™ i5 processor running Ubuntu 15.01 with a 4 GB RAM.

4.1 Performance in the scalar field 4.1.1 Runtime vs. signal size

Here, the signal sparsity parameter is fixed to k=50, and we compute the runtime of each algorithm

for five different signal sizes: 2¹⁶, 2¹⁷, 2¹⁸, 2¹⁹, and 2²⁰. Table 1 summarizes the runtimes for processing various signal sizes using SFCFT, sFFT1.0, sFFT2.0, FCFT, and FFTW. Fig. 3 shows comparisons among these algorithms.

As expected, the runtime of each algorithm is approximately linear (Fig. 3). However, Fig. 3 demonstrates that for signal size n>2¹⁶, SFCFT is faster than FFTW and FCFT when they have to

Algorithm 1 Sparse fast Clifford Fourier transform Input: a multivector signal

Output: result of SFCFT

1. Extract the Clifford basis components of the multivector 2. Reconstruct the components to obtain vectors

0 123 3

1 23 3

2 31 3

3 12 3

[ ( ) ( ) ] [ ( ) ( ) ] [ ( ) ( ) ] [ ( ) ( ) ]

F x F x i

F x F x i

§ ·

¨ ¸

© ¹

1 2 3

1 e e e

// For each vector, it is appropriate to use sFFT, // as it is sparse

3. Run a number of location loops, returning L sets of coor- dinates I1, I2, …, IL

4. Randomly choose a ı invertible mod n and Ĳ[n], and permute input vector f with permutation Tı,Ĳ: (Tı,Ĳf)j=fıj+Ĳ

5. Use a flat window function G, and compute the filtered and permuted vector y=Gā(Tı,Ĳf)

6. With B dividing n, compute zˆ yˆ_{j n B}_{( / )} for j[B]

7. Keep only the dāk coordinates of the maximum magnitude in ˆz

8. Reverse steps 5–7 as a hash function h_ı: [n]ĺ[B]

9. Count the number of occurrences of each coordinate j found, i.e., sj=|{r|jIr}|

10. Keep only the coordinates that occur in at least half of the location loops:

Iƍ={jI1I2…IL|sj>L/2}

11. Run a number of estimation loops on Iƍ, and return sets of frequency coefficients ˆf_j^r

// These loops work similarly to location loops in // steps 4–6 and the next step is: given a set of coordinates // I, estimate ˆ

fj as fˆ_jc zˆ_h_V_{( )}_jZ^W^j/Gˆ_V_{( )}_j

12. Take median real and imaginary components separately and estimate each frequency coefficient ˆ

fj as median{

ˆ_j f_j^r|r {1,2, ,L}}

f c "

13. Based on the linearity of CFT, add up ˆ

fj to obtain the output

(8)

recover exactly 50 nonzero coefficients. Thus, we can see that the slopes of the lines for FFTW and FCFT are larger than those for sFFT1.0, sFFT2.0, and SFCFT. It means that FFTW and FCFT do not efficiently process signals with large sizes. SFCFT is slower than sFFT1.0 and sFFT2.0, because SFCFT spends a certain amount of time for Clifford algebra operations.

4.1.2 Runtime vs. sparsity (k)

In this experiment, the signal size is fixed to n=2²⁰, and the runtime of each algorithm is calculated.

We choose sparsity parameter k from 50 to 500. The experiment is repeated 50 times for each value of k.

Table 2 summarizes the runtimes for processing various signal sparsities using SFCFT, sFFT1.0, sFFT2.0, FCFT, and FFTW. Fig. 4 illustrates the comparison among these algorithms.

Fig. 4 shows that the runtime of each algorithm is approximately linear. Note that FFTW and FCFT are independent of sparsity k. It can be seen that when a sparse approximation of the Fourier transform is available, SFCFT, sFFT1.0, and sFFT2.0 greatly improve the runtime performances and extend the range of applications. Further, SFCFT, sFFT1.0, and sFFT2.0 are faster than FFTW and FCFT when spar- sity parameter k is small enough (Fig. 4).

4.1.3 Robustness to noise

Here, we choose signal size n=2²⁰ and sparsity parameter k=50. For robustness, we add Gaussian noise to signals and run each experiment with different signal-to-noise ratios (SNRs). SNR is changed by changing the Gaussian noise. For each variance of the Gaussian noise, we regenerate new cases of noises and signals to run each experiment 5000 times. For each run, we compute the error metric as the average error between the best k-sparse approximation of ˆx and the output approximation ˆxc which is restricted to its k largest entries. This allows us to check whether our SFCFT can improve its performance in terms of runtime without reducing its robustness to noise.

The average errors per entry for SFCFT, sFFT1.0, sFFT 2.0, FCFT, and FFTW are shown in Fig. 5. It can be seen that the performances of all these algorithms are stable under noise, and SFCFT is more robust to noise than FFTW and FCFT.

4.2 Performance using 2D grayscale images In this experiment, we apply our algorithm in 2D grayscale image Fourier transform and compare it with FFT and CFT. We choose a ‘Lena’ grayscale image (512×512) to evaluate the performance. Re- sults show that our SFCFT is the fastest, since SFCFT computes data directly in the spatial domain first, and then computes the effective coefficients.

SFCFT sFFT2.0

2¹⁶ 2¹⁷ 2¹⁸ 2¹⁹ 2²⁰

Signal size

Runtime (ms)

sFFT1.0 FFTW FCFT 6

4 2 0

Fig. 3 Runtimes of the compared algorithms for different signal sizes

Table 1 Runtimes of each algorithm at different signal sizes

Signal size

Runtime (ms)

SFCFT sFFT2.0 sFFT1.0 FFTW FCFT 2¹⁶ 0.4031 0.3471 0.3615 0.2028 0.8804 2¹⁷ 0.6234 0.5218 0.5804 0.4405 1.0615 2¹⁸ 0.8098 0.7048 0.7595 1.1626 1.6095 2¹⁹ 1.5343 0.9578 1.2825 2.1543 2.8825 2²⁰ 3.3743 1.8775 2.5088 4.5130 5.5088

Table 2 Runtimes of each algorithm with different sparsity parameters

k Runtime (ms)

SFCFT sFFT2.0 sFFT1.0 FFTW FCFT 50 3.3743 1.8875 2.1088 4.5130 5.7505 100 3.7175 2.9301 3.3965 4.5130 5.7505 200 4.6654 3.8355 4.2505 4.5130 5.7505 500 6.5449 5.3826 5.8371 4.5130 5.7505

SFCFT sFFT2.0 sFFT1.0 FFTW FCFT

50 100 200 500

Sparsity k

Runtime (ms)

8 6 4 2 0

Fig. 4 Runtimes of the compared algorithms for different sparsity parameters

(9)

Further, Figs. 6b–6d show that SFCFT outputs a clearer spectrum map than FFT and CFT. It means that SFCFT presents more image data information. As mentioned before, compared with the traditional Fourier transform algorithm, our SFCFT algorithm has lower error rates. SFCFT transforms the original field into a field where nonzero Fourier coefficients exist. This improves SFCFT’s robustness to noise.

4.3 Performance using color images

In this experiment, we use SFCFT in processing higher-dimensional signals. We apply our algorithm in color-image processing and compare the results with those using FFT and CFT. We choose a ‘Lena’

RGB image (512×512) to evaluate the performance of the algorithm. An RGB image has three color-

channels (i.e., red, green, and blue channels); thus, FFT is calculated for each color-channel separately. It means that three grayscale Fourier transforms are calculated for each channel.

However, SFCFT can extract the multivectors from the color-image data and transform them into complex signals by implementing the Clifford algebra and then calculating their Fourier transformations in a numerical way. A 3D vector field is transformed into a multivector field with only bivector and vector parts that are unequal to zero, as bivector and vector form three complex signals. Thus, we can see from Fig. 7 that the SFCFT algorithm is much faster than FFT and CFT, and its frequency spectrum map is clearer. It means that SFCFT presents more image data information. SFCFT can discover both geometric and spectral information of the multispectral image. The wider the bands of the multispectral image, the more outstanding the SFCFT.

5 Conclusions

In this paper, we have proposed a novel algorithm called the ‘sparse fast Clifford Fourier transform’, which combines the sparse Fourier transform with CFT to process multivector signals. We have Fig. 7 The original RGB image (a) and the frequency spectrum maps with FFT (b), CFT (c), and SFCFT (d) The runtimes of FFT, CFT, and SFCFT are 0.195382, 0.774559, and 0.016852 s, respectively. References to color refer to the online version of this figure

Fig. 6 The original image (a) and the frequency spectrum maps with FFT (b), CFT (c), and SFCFT (d)

The runtimes of FFT, CFT, and SFCFT are 0.029599, 0.054489, and 0.006536 s, respectively

SFCFT sFFT2.0 sFFT1.0 FFTW FCFT

Average error per entry 10^í5 10^í4 10^í3 10^í2 10^í1 10⁰ 10¹

Signal-to-noise ratio (dB)

Fig. 5 Robustness to noise of the compared algorithms

(10)

discussed the application of SFCFT to the scalar field and grayscale and color image processing. The results demonstrate that SFCFT can effectively improve the performance of multivector signal processing. Our method is available in most vector fields; however, it may not be sparse for a large number of vectors concentrated in specific subareas. As next steps, we will measure the uncertainties to indicate the im- portance of matches, and offer some improvements to the preliminary segmentation and transformation computed on each segment.

References

Batard, T., Berthier, M., Saint-Jean, C., 2010. Clifford–Fourier transform for color image processing. In: Bayro- Corrochano, E., Scheuermann, G. (Eds.), Geometric Al- gebra Computing. Springer London, London, UK, p.135-162. https://doi.org/10.1007/978-1-84996-108-0_8 Bujack, R., Scheuermann, G., Hitzer, E., 2015. Demystifica- tion of the geometric Fourier transforms and resulting convolution theorems. Math. Meth. Appl. Sci., 39(7):

1877-1890. https://doi.org/10.1002/mma.3607

Cao, W.M., Feng, H., 2010. Geometric Algebra in Biomimetic Pattern Recognition and Signal Processing. Science Press, Beijing, China (in Chinese).

Cao, W.M., Liu, H., Xu, C., et al., 2013. 3D medical image registration based on conformal geometric algebra. Sci.

China Inform. Sci., 43(2):254-274.

https://doi.org/10.1360/112012-592

Ebling, J., Scheuermann, G., 2003. Clifford convolution and pattern matching on vector fields. IEEE Visualization, p.193-200.

https://doi.org/10.1109/VISUAL.2003.1250372

Ebling, J., Scheuermann, G., 2005. Clifford Fourier transform on vector fields. IEEE Trans. Visual. Comput. Graph., 11(4):469-479. https://doi.org/10.1109/TVCG.2005.54 Evans, C.J., Sangwine, S.J., Ell, T.A., 2000. Colour-sensitive

edge detection using hypercomplex filters. 10th European Signal Processing Conf., p.1-4.

Frigo, M., Johnson, S.G., 2005. The design and implementa- tion of FFTW3. Proc. IEEE, 93(2):216-231.

https://doi.org/10.1109/JPROC.2004.840301

Gilbert, A., Indyk, P., 2010. Sparse recovery using sparse matrices. Proc. IEEE, 98(6):937-947.

https://doi.org/10.1109/JPROC.2010.2045092

Gilbert, A.C., Indyk, P., Iwen, M., et al., 2014. Recent devel- opments in the sparse Fourier transform: a compressed Fourier transform for big data. IEEE Signal Process.

Mag., 31(5):91-100.

https://doi.org/10.1109/MSP.2014.2329131

Hassanieh, H., Adib, F., Katabi, D., et al., 2012a. Faster GPS via the sparse Fourier transform. Proc. 18th Annual Int.

Conf. on Mobile Computing and Networking, p.353-364.

https://doi.org/10.1145/2348543.2348587

Hassanieh, H., Indyk, P., Katabi, D., et al., 2012b. Nearly optimal sparse Fourier transform. Proc. 44th Annual ACM Symp. on Theory of Computing, p.563-578.

https://doi.org/10.1145/2213977.2214029

Hassanieh, H., Indyk, P., Katabi, D., et al., 2012c. Simple and practical algorithm for sparse Fourier transform. Proc.

23rd Annual ACM-SIAM Symp. on Discrete Algorithms, p.1183-1194.

https://doi.org/10.1137/1.9781611973099.93

Hassanieh, H., Shi, L., Abari, O., et al., 2014. GHz-wide sensing and decoding using the sparse Fourier transform.

Proc. IEEE INFOCOM, p.2256-2264.

https://doi.org/10.1109/INFOCOM.2014.6848169 Hestenes, D., 1999. New Foundations for Classical Mechanics.

Springer, New York, USA.

Hestenes, D., Sobczyk, G., 1984. Clifford Algebra to Geomet- ric Calculus. Springer, New York, USA.

Hitzer, E., 2012. Introduction to Clifford’s geometric algebra.

SICE J. Contr. Meas. Syst. Integr., 4(1):1-11.

Hitzer, E., 2013. The Clifford Fourier transform in real Clifford algebras. Int. Conf. on the Applications of Computer Science and Mathematics in Architecture and Civil En- gineering, p.227-240.

Hitzer, E., Mawardi, B., 2008. Clifford Fourier transform on multivector fields and uncertainty principles for dimen- sions n=2 (mod 4) and n=3 (mod 4). Adv. Appl. Clifford Alg., 18(3-4):715-736.

https://doi.org/10.1007/s00006-008-0098-3

Hitzer, E., Sangwine, S.J., 2013. Quaternion and Clifford Fourier transforms and wavelets. In: Hitzer, E., Stephen, J., Sangwine, S.J. (Eds.), Trends in Mathematic. Springer Basel, Basel, Switzerland.

https://doi.org/10.1007/978-3-0348-0603-9

Hitzer, E., Helmstetter, J., Abáamowicz, R., 2013. Square roots of –1 in real Clifford algebras. In: Hitzer, E., Stephen, J., Sangwine, S.J. (Eds.), Trends in Mathematic. Springer Basel, Basel, Switzerland.

https://doi.org/10.1007/978-3-0348-0603-9_7

Iwen, M.A., 2010. Combinatorial sublinear-time Fourier al- gorithms. Found. Comput. Math., 10(3):303-338.

https://doi.org/10.1007/s10208-009-9057-1 Li, H.B., 2005. Conformal geometric algebra—a new frame-

work for computational geometry. J. Comput. Aid. Des.

Comput. Graph., 17(11):2383-2393.

Mishra, B., Wilson, P., Wilcock, R., 2015. A geometric algebra co-processor for color edge detection. Electronics, 4(1):

94-117. https://doi.org/10.3390/electronics4010094 Mueen, A., Nath, S., Liu, J., 2010. Fast approximate correla-

tion for massive time-series data. Proc. ACM SIGMOD Int. Conf. on Management of Data, p.171-182.

https://doi.org/10.1145/1807167.1807188

Porat, E., Strauss, M.J., 2012. Sublinear time, measurement- optimal, sparse recovery for all. Proc. 23rd Annual ACM-SIAM Symp. on Discrete Algorithms, p.1215- 1227. https://doi.org/10.1137/1.9781611973099.96

(11)

Reich, W., Scheuermann, G., 2010. Analyzing real vector fields with Clifford convolution and Clifford–Fourier transform. In: Bayro-Corrochano, E., Scheuermann, G.

(Eds.), Geometric Algebra Computing. Springer, London, UK, p.121-133.

https://doi.org/10.1007/978-1-84996-108-0_7

Schlemmer, M., Hotz, I., Natarajan, V., et al., 2005. Fast Clifford Fourier transformation for unstructured vector field data. Proc. Int. Conf. on Numerical Grid Generation in Computational Field Simulations, p.101-110.

Schumacher, J., Puschel, M., 2014. High-performance sparse fast Fourier transforms. IEEE Workshop on Signal Pro- cessing Systems, p.1-6.

https://doi.org/10.1109/SiPS.2014.6986055

Shi, L., Hassanieh, H., Davis, A., et al., 2014. Light field reconstruction using sparsity in the continuous Fourier domain. ACM Trans. Graph., 34(1), Article 12.

https://doi.org/10.1145/2682631

Smith, J.O., 2011. Spectral Audio Signal Processing. W3K Publishing, London, UK.

Wang, R., Jing, L.B., Tao, L., et al., 2013. Digital watermark- ing algorithm for 3D point cloud model based on Clifford algebra. J. Shanghai Jiao Tong Univ., 47(12):1863-1869.

Wang, R., Zhang, X., Cao, W.M., 2015. Clifford fuzzy support vector machines for classification. Adv. Appl. Clifford Alg., 26(2):1-22.

https://doi.org/10.1007/s00006-015-0616-z

Xie, W., Cao, W.M., Meng, S., 2008. Coverage analysis for sensor networks based on Clifford algebra. Sci. China Inform. Sci., 51(5):460-475.

https://doi.org/10.1007/s11432-008-0048-7 Xu, C., Liu, H., Ouyang, C.J., et al., 2011. Theory and appli-

cation of Clifford pseudo-differential operator on multi- spectral image. Sci. Sin. Inform., 41(12):1423-1435.

https://doi.org/10.1360/zf2011-41-12-1423