- Email: [email protected]
Fast 3-D decimation-in-frequency algorithm for 3-D Hartley transform
Signal Processing 82 (2002) 121–126
www.elsevier.com/locate/sigpro
Fast communication
Fast 3-D decimation-in-frequency algorithm for 3-D Hartley transform O. Alshibami, S. Boussakta ∗ Institute of Integrated Information Systems, School of Electronic and Electrical Engineering, The University of Leeds, Leeds, LS2 9JT, UK Received 1 March 2001
Abstract The three-dimensional discrete Hartley transform (3-D DHT) has been applied in a wide range of 3-D applications such as 3-D power spectrum analysis, 3-D 0ltering, and medical applications, etc. In this paper, a three-dimensional algorithm for fast computation of the three-dimensional discrete Hartley transform is developed. The mathematical concept and derivation is presented and the arithmetic complexity is analysed and compared to the familiar row– column approach. It is found that this algorithm o3ers substantial savings in both the number of multiplications and additions. ? 2002 Elsevier Science B.V. All rights reserved.
1. Introduction The three-dimensional discrete Hartley transform (3-D DHT) [1,6] has been applied in many 3-D applications such as 3-D spectrum analysis, signal and image processing, motion analysis, multidimensional 0ltering, etc. [2,4,5]. The 3-D DHT is usually calculated using the row–column approach [3]. However, proper 3-D algorithms are faster and require to be developed. In this letter, the 3-D vector radix decimation-in-frequency algorithm (3-D V-DIF) for fast calculation of the 3-D DHT is developed. The new algorithm is implemented and its arithmetic complexity is analysed and compared to the familiar row–column approach. It is found that, the 3-D V-DIF o3ers a saving of more than 40% in multiplications and about (20 –25)% in additions.
∗
Corresponding author. Tel.: +44-113–233-2008; fax: +44-113-323-2032. E-mail address: [email protected] (S. Boussakta).
0165-1684/02/$ - see front matter ? 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 1 6 8 4 ( 0 1 ) 0 0 1 0 2 - 5
122
O. Alshibami, S. Boussakta / Signal Processing 82 (2002) 121–126
2. Transform denition The forward 3-D DHT of a real 3-D input data, x(n1 ; n2 ; n3 ), is de0ned by X (k1 ; k2 ; k3 ) =
N 3 −1 1 −1 N 2 −1 N n1 =0 n2 =0
2 2 2 x(n1 ; n2 ; n3 ) cas n1 k 1 + n2 k2 + n3 k3 N N N3 1 2 n =0
3
k1; 2; 3 = 0; 1; 2; : : : ; N1; 2; 3
(1)
where cas( + + ) = cos( + + ) + sin( + + ). The inverse 3-D DHT only di3ers from the forward 3-DHT by a factor 1=N1 N2 N3 . 3. Fast algorithm for the 3-D DHT Owing to the diBculty of developing fast multidimensional algorithms in three and more dimensions and the fact that the multidimensional Hartley transform is not separable, the 3-D DHT is usually computed, through an intermediate transform, using algorithms developed for the 1-D case applied over each dimension in row–column approach [2–5]. In this section, the 3-D vector-radix decimation in frequency algorithm is introduced for fast calculation of the 3-D DHT. The transform size should be equal to powers of 2 × 2 × 2. Unlike the row–column approach, the 3-D vector-radix algorithm involves elements from throughout the entire three-dimensional array rather than operating on individual rows and columns. 3.1. Mathematical development In the 3-D vector-radix decimation-in-frequency algorithm, the 3-D input, x(n1 ; n2 ; n3 ), of size N ×N ×N is partitioned into eight summations each of size N=2 × N=2 × N=2 as follows: X (k1 ; k2 ; k3 ) = X000 (k1 ; k2 ; k3 ) + X001 (k1 ; k2 ; k3 ) + X010 (k1 ; k2 ; k3 ) +X011 (k1 ; k2 ; k3 ) + X100 (k1 ; k2 ; k3 ) + X101 (k1 ; k2 ; k3 ) +X110 (k1 ; k2 ; k3 ) + X111 (k1 ; k2 ; k3 );
(2)
where
N=2−1 N=2−1 N=2−1
Xabc (k1 ; k2 ; k3 ) =
n1 =0 n2 =0
2 ×cas N
N N N x n1 + a ; n2 + b ; n3 + c 2 2 2 n =0 3
N n1 + a 2
N k 1 + n2 + b 2
N k 2 + n3 + c 2
k3
(3)
and abc = {000; 001; 010; 011; 100; 101; 110; 111}. Using the identity cas(m + n) = cas(m) cos(n) + cas(−m) sin(n)
(4)
and considering the even and odd parts of k1 , k2 , and k3 , Eq. (2) can be decomposed into eight equations as:
O. Alshibami, S. Boussakta / Signal Processing 82 (2002) 121–126
123
For k1 -even, k2 -even, and k3 -even, Eq. (2) can be written as N=2−1 N=2−1 N=2−1 N X (2k1 ; 2k2 ; 2k3 ) = x (n1 ; n2 ; n3 ) + x n1 ; n2 ; n3 + 2 n =0 n =0 n =0 1
2
3
N N N ; n3 + x n1 ; n2 + ; n3 + 2 2 2 N N N +x n1 + ; n2 ; n3 + x n1 + ; n2 ; n3 + 2 2 2 N N N N N +x n1 + ; n2 + ; n3 + x n1 + ; n2 + ; n3 + 2 2 2 2 2 2 × cas (2n1 k1 + 2n2 k2 + 2n3 k3 ) N
+x n1 ; n2 +
N − 1; i = 1; 2; 3: 2 Eq. (5) can be identi0ed as an N=2 × N=2 × N=2-point 3-D DHT. For k1 -even, k2 -even, and k3 -odd, Eq. (2) can be written as N=2−1 N=2−1 N=2−1 N X (2k1 ; 2k2 ; 2k3 + 1) = x (n1 ; n2 ; n3 ) − x n1 ; n2 ; n3 2 n =0 n =0 n =0 ki = 0; 1; : : : ;
1
2
3
N N N +x n1 ; n2 + ; n3 − x n1 ; n2 + ; n3 + 2 2 2 N N N +x n1 + ; n2 ; n3 − x n1 + ; n2 ; n3 + 2 2 2 N N N N N +x n1 + ; n2 + ; n3 − x n1 + ; n2 + ; n3 + 2 2 2 2 2 2 ×cas (2n1 k1 + 2n2 k2 + n3 (2k3 + 1)) : N
Using the identity in Eq. (4), and the fact that N=2−1 N=2−1 N 2 2 2nk x(n) cas − 2nk = x − n cas N 2 N n=0 n=0
2
3
N N ; n3 − x n1 ; n2 + ; n3 + 2 2 N N +x n1 + ; n2 ; n3 − x n1 + ; n2 ; n3 + 2 2
+x n1 ; n2 +
(6)
(7)
Eq. (6) can be decomposed as N=2−1 N=2−1 N=2−1 N x(n1 ; n2 ; n3 ) − x n1 ; n2 ; n3 + X (2k1 ; 2k2 ; 2k3 + 1) = 2 n =0 n =0 n =0 1
(5)
N 2 N 2
124
O. Alshibami, S. Boussakta / Signal Processing 82 (2002) 121–126
N N N N N +x n1 + ; n2 + ; n3 − x n1 + ; n2 + ; n3 + 2 2 2 2 2 N N 2 N × cos n3 + x − n1 ; − n2 ; − n3 N 2 2 2 N N N N −x − n1 ; − n2 ; N − n3 + x − n1 ; N − n2 ; − n3 2 2 2 2 N N N −x − n1 ; N − n2 ; N − n3 + x N − n1 ; − n2 ; − n3 2 2 2 N N −x N − n1 ; − n2 ; N − n3 + x N − n1 ; N − n2 ; − n3 2 2 2 − x (N − n1 ; N − n2 ; N − n3 )) sin n3 N 2 ×cas (2n1 k1 + 2n2 k2 + 2n3 k3 ) N
N − 1; i = 1; 2; 3: (8) 2 Again Eq. (8) can be identi0ed as an N=2 × N=2 × N=2-point 3-D DHT. Following the same procedure X (2k1 ; 2k2 +1; 2k3 ); X (2k1 ; 2k2 +1; 2k3 +1); X (2k1 +1; 2k2 ; 2k3 ); X (2k1 + 1; 2k2 ; 2k3 + 1); X (2k1 + 1; 2k2 + 1; 2k3 ); and X (2k1 + 1; 2k2 + 1; 2k3 + 1) can be calculated. Hence in this algorithm, an N × N × N 3-DHT is divided into eight 3-D DHTs of size N=2 × N=2 × N=2 each, plus some additions and multiplications by the twiddle factors. In the next stage, each 3-D DHT of size N=2 × N=2 × N=2 is further divided into eight 3-D transforms of size N=4 × N=4 × N=4 and the process continues until we get transforms of size 2 × 2 × 2. ki = 0; 1; : : : ;
3.2. Arithmetic complexity For in place calculation, two butterFies are combined together as shown in Fig. 1. This calculates 16 points and involves 28 real multiplications and 62 real additions. The whole 3-D DHT calculation needs log2 N stages. Therefore, the calculation of the whole 3-D DHT (including trivial multiplications and additions) requires (7=4)N 3 log2 N real multiplications and (31=8)N 3 log2 N real additions. On the other hand, if the 3-D DHT is calculated using the row–column approach, the total number of real multiplications, without counting the halvings (multiplication by 12 ), will be 3N 3 log2 N and the total number of real additions will be (9=2)N 3 log2 N + 3N 3 . Fig. 2 shows the number of multiplications plus additions per point, involved in the calculation of the 3-D DHT, using the new algorithm and the row–column approach. 4. Conclusion In this letter, the 3-D vector-radix decimation-in-frequency algorithm for fast calculation of 3-D discrete Hartley transform has been developed and implemented. Its arithmetic complexity has been analysed and
O. Alshibami, S. Boussakta / Signal Processing 82 (2002) 121–126
125
Fig. 1. 3-D vector-radix DIF butterFy for 3-D DHT.
Fig. 2. Arithmetic complexity comparison.
compared to similar algorithms. Compared to the row–column approach based on a single butterFy implementation, the new algorithm reduces the number of multiplications by more than 40% and reduces the number of additions by about (20 –25)%. The advantages/disadvantages of both algorithms, in terms of indexing and regularity of memory access, etc., have not been discussed in this paper. References [1] R.N. Bracewell, The Discrete Hartley Transform, Oxford University Press, Oxford, 1986. [2] A. Erdi, E. Yorke, M. Loew, Y. Erdi, M. Sarfaraz, B. Wessels, Use of the fast Hartley transform for three-dimensional dose calculation in radionuclide therapy, Med. Phys. 25 (1998) 2226–2233. [3] H. Hao, R.N. Bracewell, A three-dimensional DFT algorithm using the fast Hartley transform, Proc. IEEE 75 (1987) 264–266.
126
O. Alshibami, S. Boussakta / Signal Processing 82 (2002) 121–126
[4] S.A. Mahmoud, Motion analysis of multiple moving objects using Hartley transform, IEEE Trans. Syst. Man Cybernet. 21 (1991) 280–287. [5] C.H. Paik, M.D. Fox, Fast Hartley transforms for image processing, IEEE Trans. Med. Imaging 7 (1988) 149–153. [6] H.V. Sorensen, D.L. Jones, C.S. Burrus, M.T. Heidman, On computing the discrete Hartley transform, IEEE Trans. ASSP-33 (1985) 1231–1238.
www.elsevier.com/locate/sigpro
Fast communication
Fast 3-D decimation-in-frequency algorithm for 3-D Hartley transform O. Alshibami, S. Boussakta ∗ Institute of Integrated Information Systems, School of Electronic and Electrical Engineering, The University of Leeds, Leeds, LS2 9JT, UK Received 1 March 2001
Abstract The three-dimensional discrete Hartley transform (3-D DHT) has been applied in a wide range of 3-D applications such as 3-D power spectrum analysis, 3-D 0ltering, and medical applications, etc. In this paper, a three-dimensional algorithm for fast computation of the three-dimensional discrete Hartley transform is developed. The mathematical concept and derivation is presented and the arithmetic complexity is analysed and compared to the familiar row– column approach. It is found that this algorithm o3ers substantial savings in both the number of multiplications and additions. ? 2002 Elsevier Science B.V. All rights reserved.
1. Introduction The three-dimensional discrete Hartley transform (3-D DHT) [1,6] has been applied in many 3-D applications such as 3-D spectrum analysis, signal and image processing, motion analysis, multidimensional 0ltering, etc. [2,4,5]. The 3-D DHT is usually calculated using the row–column approach [3]. However, proper 3-D algorithms are faster and require to be developed. In this letter, the 3-D vector radix decimation-in-frequency algorithm (3-D V-DIF) for fast calculation of the 3-D DHT is developed. The new algorithm is implemented and its arithmetic complexity is analysed and compared to the familiar row–column approach. It is found that, the 3-D V-DIF o3ers a saving of more than 40% in multiplications and about (20 –25)% in additions.
∗
Corresponding author. Tel.: +44-113–233-2008; fax: +44-113-323-2032. E-mail address: [email protected] (S. Boussakta).
0165-1684/02/$ - see front matter ? 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 1 6 8 4 ( 0 1 ) 0 0 1 0 2 - 5
122
O. Alshibami, S. Boussakta / Signal Processing 82 (2002) 121–126
2. Transform denition The forward 3-D DHT of a real 3-D input data, x(n1 ; n2 ; n3 ), is de0ned by X (k1 ; k2 ; k3 ) =
N 3 −1 1 −1 N 2 −1 N n1 =0 n2 =0
2 2 2 x(n1 ; n2 ; n3 ) cas n1 k 1 + n2 k2 + n3 k3 N N N3 1 2 n =0
3
k1; 2; 3 = 0; 1; 2; : : : ; N1; 2; 3
(1)
where cas( + + ) = cos( + + ) + sin( + + ). The inverse 3-D DHT only di3ers from the forward 3-DHT by a factor 1=N1 N2 N3 . 3. Fast algorithm for the 3-D DHT Owing to the diBculty of developing fast multidimensional algorithms in three and more dimensions and the fact that the multidimensional Hartley transform is not separable, the 3-D DHT is usually computed, through an intermediate transform, using algorithms developed for the 1-D case applied over each dimension in row–column approach [2–5]. In this section, the 3-D vector-radix decimation in frequency algorithm is introduced for fast calculation of the 3-D DHT. The transform size should be equal to powers of 2 × 2 × 2. Unlike the row–column approach, the 3-D vector-radix algorithm involves elements from throughout the entire three-dimensional array rather than operating on individual rows and columns. 3.1. Mathematical development In the 3-D vector-radix decimation-in-frequency algorithm, the 3-D input, x(n1 ; n2 ; n3 ), of size N ×N ×N is partitioned into eight summations each of size N=2 × N=2 × N=2 as follows: X (k1 ; k2 ; k3 ) = X000 (k1 ; k2 ; k3 ) + X001 (k1 ; k2 ; k3 ) + X010 (k1 ; k2 ; k3 ) +X011 (k1 ; k2 ; k3 ) + X100 (k1 ; k2 ; k3 ) + X101 (k1 ; k2 ; k3 ) +X110 (k1 ; k2 ; k3 ) + X111 (k1 ; k2 ; k3 );
(2)
where
N=2−1 N=2−1 N=2−1
Xabc (k1 ; k2 ; k3 ) =
n1 =0 n2 =0
2 ×cas N
N N N x n1 + a ; n2 + b ; n3 + c 2 2 2 n =0 3
N n1 + a 2
N k 1 + n2 + b 2
N k 2 + n3 + c 2
k3
(3)
and abc = {000; 001; 010; 011; 100; 101; 110; 111}. Using the identity cas(m + n) = cas(m) cos(n) + cas(−m) sin(n)
(4)
and considering the even and odd parts of k1 , k2 , and k3 , Eq. (2) can be decomposed into eight equations as:
O. Alshibami, S. Boussakta / Signal Processing 82 (2002) 121–126
123
For k1 -even, k2 -even, and k3 -even, Eq. (2) can be written as N=2−1 N=2−1 N=2−1 N X (2k1 ; 2k2 ; 2k3 ) = x (n1 ; n2 ; n3 ) + x n1 ; n2 ; n3 + 2 n =0 n =0 n =0 1
2
3
N N N ; n3 + x n1 ; n2 + ; n3 + 2 2 2 N N N +x n1 + ; n2 ; n3 + x n1 + ; n2 ; n3 + 2 2 2 N N N N N +x n1 + ; n2 + ; n3 + x n1 + ; n2 + ; n3 + 2 2 2 2 2 2 × cas (2n1 k1 + 2n2 k2 + 2n3 k3 ) N
+x n1 ; n2 +
N − 1; i = 1; 2; 3: 2 Eq. (5) can be identi0ed as an N=2 × N=2 × N=2-point 3-D DHT. For k1 -even, k2 -even, and k3 -odd, Eq. (2) can be written as N=2−1 N=2−1 N=2−1 N X (2k1 ; 2k2 ; 2k3 + 1) = x (n1 ; n2 ; n3 ) − x n1 ; n2 ; n3 2 n =0 n =0 n =0 ki = 0; 1; : : : ;
1
2
3
N N N +x n1 ; n2 + ; n3 − x n1 ; n2 + ; n3 + 2 2 2 N N N +x n1 + ; n2 ; n3 − x n1 + ; n2 ; n3 + 2 2 2 N N N N N +x n1 + ; n2 + ; n3 − x n1 + ; n2 + ; n3 + 2 2 2 2 2 2 ×cas (2n1 k1 + 2n2 k2 + n3 (2k3 + 1)) : N
Using the identity in Eq. (4), and the fact that N=2−1 N=2−1 N 2 2 2nk x(n) cas − 2nk = x − n cas N 2 N n=0 n=0
2
3
N N ; n3 − x n1 ; n2 + ; n3 + 2 2 N N +x n1 + ; n2 ; n3 − x n1 + ; n2 ; n3 + 2 2
+x n1 ; n2 +
(6)
(7)
Eq. (6) can be decomposed as N=2−1 N=2−1 N=2−1 N x(n1 ; n2 ; n3 ) − x n1 ; n2 ; n3 + X (2k1 ; 2k2 ; 2k3 + 1) = 2 n =0 n =0 n =0 1
(5)
N 2 N 2
124
O. Alshibami, S. Boussakta / Signal Processing 82 (2002) 121–126
N N N N N +x n1 + ; n2 + ; n3 − x n1 + ; n2 + ; n3 + 2 2 2 2 2 N N 2 N × cos n3 + x − n1 ; − n2 ; − n3 N 2 2 2 N N N N −x − n1 ; − n2 ; N − n3 + x − n1 ; N − n2 ; − n3 2 2 2 2 N N N −x − n1 ; N − n2 ; N − n3 + x N − n1 ; − n2 ; − n3 2 2 2 N N −x N − n1 ; − n2 ; N − n3 + x N − n1 ; N − n2 ; − n3 2 2 2 − x (N − n1 ; N − n2 ; N − n3 )) sin n3 N 2 ×cas (2n1 k1 + 2n2 k2 + 2n3 k3 ) N
N − 1; i = 1; 2; 3: (8) 2 Again Eq. (8) can be identi0ed as an N=2 × N=2 × N=2-point 3-D DHT. Following the same procedure X (2k1 ; 2k2 +1; 2k3 ); X (2k1 ; 2k2 +1; 2k3 +1); X (2k1 +1; 2k2 ; 2k3 ); X (2k1 + 1; 2k2 ; 2k3 + 1); X (2k1 + 1; 2k2 + 1; 2k3 ); and X (2k1 + 1; 2k2 + 1; 2k3 + 1) can be calculated. Hence in this algorithm, an N × N × N 3-DHT is divided into eight 3-D DHTs of size N=2 × N=2 × N=2 each, plus some additions and multiplications by the twiddle factors. In the next stage, each 3-D DHT of size N=2 × N=2 × N=2 is further divided into eight 3-D transforms of size N=4 × N=4 × N=4 and the process continues until we get transforms of size 2 × 2 × 2. ki = 0; 1; : : : ;
3.2. Arithmetic complexity For in place calculation, two butterFies are combined together as shown in Fig. 1. This calculates 16 points and involves 28 real multiplications and 62 real additions. The whole 3-D DHT calculation needs log2 N stages. Therefore, the calculation of the whole 3-D DHT (including trivial multiplications and additions) requires (7=4)N 3 log2 N real multiplications and (31=8)N 3 log2 N real additions. On the other hand, if the 3-D DHT is calculated using the row–column approach, the total number of real multiplications, without counting the halvings (multiplication by 12 ), will be 3N 3 log2 N and the total number of real additions will be (9=2)N 3 log2 N + 3N 3 . Fig. 2 shows the number of multiplications plus additions per point, involved in the calculation of the 3-D DHT, using the new algorithm and the row–column approach. 4. Conclusion In this letter, the 3-D vector-radix decimation-in-frequency algorithm for fast calculation of 3-D discrete Hartley transform has been developed and implemented. Its arithmetic complexity has been analysed and
O. Alshibami, S. Boussakta / Signal Processing 82 (2002) 121–126
125
Fig. 1. 3-D vector-radix DIF butterFy for 3-D DHT.
Fig. 2. Arithmetic complexity comparison.
compared to similar algorithms. Compared to the row–column approach based on a single butterFy implementation, the new algorithm reduces the number of multiplications by more than 40% and reduces the number of additions by about (20 –25)%. The advantages/disadvantages of both algorithms, in terms of indexing and regularity of memory access, etc., have not been discussed in this paper. References [1] R.N. Bracewell, The Discrete Hartley Transform, Oxford University Press, Oxford, 1986. [2] A. Erdi, E. Yorke, M. Loew, Y. Erdi, M. Sarfaraz, B. Wessels, Use of the fast Hartley transform for three-dimensional dose calculation in radionuclide therapy, Med. Phys. 25 (1998) 2226–2233. [3] H. Hao, R.N. Bracewell, A three-dimensional DFT algorithm using the fast Hartley transform, Proc. IEEE 75 (1987) 264–266.
126
O. Alshibami, S. Boussakta / Signal Processing 82 (2002) 121–126
[4] S.A. Mahmoud, Motion analysis of multiple moving objects using Hartley transform, IEEE Trans. Syst. Man Cybernet. 21 (1991) 280–287. [5] C.H. Paik, M.D. Fox, Fast Hartley transforms for image processing, IEEE Trans. Med. Imaging 7 (1988) 149–153. [6] H.V. Sorensen, D.L. Jones, C.S. Burrus, M.T. Heidman, On computing the discrete Hartley transform, IEEE Trans. ASSP-33 (1985) 1231–1238.