Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Some software and hardware implementations of the fast Hartley transform Fu, Yan Kit 1990

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-UBC_1991_A7 F82.pdf [ 5.83MB ]
Metadata
JSON: 831-1.0065455.json
JSON-LD: 831-1.0065455-ld.json
RDF/XML (Pretty): 831-1.0065455-rdf.xml
RDF/JSON: 831-1.0065455-rdf.json
Turtle: 831-1.0065455-turtle.txt
N-Triples: 831-1.0065455-rdf-ntriples.txt
Original Record: 831-1.0065455-source.json
Full Text
831-1.0065455-fulltext.txt
Citation
831-1.0065455.ris

Full Text

SOME SOFTWARE AND HARDWARE IMPLEMENTATIONS OF T H E FAST HARTLEY TRANSFORM By F U , Y A N KIT B.Sc.Eng., University of New Brunswick, 1986 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF M A S T E R OF A P P L I E D S C I E N C E in THE FACULTY OF GRADUATE STUDIES D E P A R T M E N T OF E L E C T R I C A L E N G I N E E R I N G We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y OF BRITISH C O L U M B I A August 1990 © F U Yan Ki t , 1990 In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Electrical Engineering The University of British Columbia 2356 Main Mall Vancouver, Canada Date: Abstract The fast Hartley transform (FHT) is a new tool for converting data between time and frequency domains. In this thesis, some speed-optimized software implementations of the radix-2 and split-radix FHT algorithms are presented initially, and then applied to the problems of convolution, computation of power spectra, image degradation, and image restoration. Subsequent work involved the development of a new bit-reversal algorithm. This algorithm is fast and efficient, and can be used to increase the throughput of the FHT. Finally, several hardware implementations are presented for the discrete Hartley trans-form (DHT) and the FHT with architectures using a single butterfly unit, pipelining and superparallelism. The advantages of each implementation are stressed. The data pro-cessing rates of these hardware implementations are analyzed. i Table of Contents Abstract i List of Tables v List of Figures v i i Acknowledgement v i i i 1 Introduction 1 1.1 Digital signal processing 1 1.2 Overview of the following chapters 2 2 The fast Hartley transform ( F H T ) software implementation 4 2.1 The Hartley transform 4 2.2 The discrete Hartley transform (DHT) 4 2.3 The F H T 5 2.3.1 The radix-2 fast algorithm 5 2.3.2 The split-radix fast algorithm 10 2.4 Comparison with the complex F F T 14 2.5 Comparison with the real F F T 18 2.6 To get the F F T from the F H T 19 2.7 F H T convolution and power spectrum 21 2.8 A n image processing example of the application of the F H T 22 2.9 Summary 26 n 3 The bit-reversal algorithms 27 3.1 Some existing algorithms 27 3.2 A new bit-reversal algorithm 31 3.3 Comparison of the bit-reversal algorithms 34 3.4 Summary 35 4 The F H T hardware implementation 37 4.1 Introduction 37 4.2 Hardware implementation of the D H T 37 4.3 Hardware implementation of the F H T 40 4.3.1 Single butterfly unit F H T architecture 40 4.3.2 Pipeline F H T architecture 45 4.3.3 Superparallel F H T architecture 49 4.4 Data unscrambling in hardware implementation 51 4.5 Comparison of the hardware implementations 52 4.6 Summary 53 5 Conclusion 54 5.1 The work done in this thesis 54 5.2 Future work 56 References 57 Appendices 61 A The radix-2 decimation-in-frequency F H T programs listing 61 A . l Radix-2 F H T implemented with the "three-loop method" 62 A.2 Radix-2 F H T implemented with the "table-lookup method" 64 iii B The split-radix decimation-in-frequency F H T programs listing 66 B . l Split-radix F H T implemented with the "three-loop method" 67 B. 2 Split-radix F H T implemented with the "table-lookup method" 70 C The radix-2 decimation-in-frequency F F T programs listing 73 C. l Radix-2 forward F F T implemented with the "three-loop method" . . . . 74 C. 2 Radix-2 forward F F T implemented with the "table-lookup method" . . . 76 D The split-radix decimation-in-frequency F F T programs listing 78 D. l Split-radix forward F F T implemented with the "three-loop method" . . . 79 D. 2 Split-radix forward F F T implemented with the "table-lookup method" . 82 E Bit-reversal programs listing 85 E. l The implementation of Evans' bit-reversal algorithm 86 E.2 The implementation of Gold and Rader's bit-reversal algorithm 88 E.3 The implementation of Rodriguez's bit-reversal algorithm 89 E.4 The implementation of Buneman's bit-reversal algorithm 90 E.5 The implementation of the new bit-reversal algorithm 91 E. 6 The implementation of the new algorithm using table-lookup method . . 92 F Image processing example using the F H T 93 F. l Using the F H T in the processing of image degradation and restoration . 94 iv L i s t of T a b l e s 2.1 Number of real arithmetic operations for the fast algorithms 15 2.2 Comparison of the fast algorithms in CPU execution time (in seconds) on a Sun 3/60 workstation, where r-2 stands for the radix-2 algorithm and s-r stands for the split-radix algorithm 15 3.3 Bit-reversing of a 16-point array 28 3.4 Exclusive-oring a bit-reversed index in column three and the related operand in column four produces the next bit-reversed index in the data sequence. 32 3.5 CPU execution time (in seconds) of the bit-reversal algorithms running on a Sun 3/60 workstation 34 4.6 Comparison of the three FHT hardware implementations for N — 1024 and memory access time of 60 ns 53 v List of Figures 2.1 Radix-2 decimation-in-frequency and decimation-in-time butterflies. . . . 7 2.2 Flow diagram of radix-2 decimation-in-frequency decomposition from a 16-point FHT into two 8-point FHT blocks 8 2.3 Flow diagram of an 8-point radix-2 decimation-in-frequency FHT 8 2.4 Decomposition pattern of a 64-point radix-2 FHT 9 2.5 Flow diagram of split-radix decimation-in-frequency decomposition from a 16-point FHT into one 8-point FHT block and two 4-point FHT blocks. 12 2.6 Flow diagram of an 8-point split-radix decimation-in-frequency FHT. . . 12 2.7 Decomposition pattern of a 64-point split-radix FHT 13 2.8 FHT requires less CPU time than the complex F F T in the radix-2 soft-ware implementation in both the "three-loop (lp)" and "table-lookup (tb)" methods 16 2.9 FHT requires less CPU time than the complex FFT in the split-radix software implementation in both the "three-loop (lp)" and "table-lookup (tb)" methods 16 2.10 Image degradation and restoration processes 23 2.11 Original image 24 2.12 Degraded image 24 2.13 Restored image 25 3.14 Gold-Rader's bit-reversal algorithm (after Rabiner and Gold) 28 3.15 The table-lookup bit-reversal algorithm (after Buneman) 30 vi 3.16 Improved Gold-Rader's algorithm (after Rodriguez) 30 3.17 The new bit-reversal algorithm 33 3.18 Plotting of Table 3.5 35 4.19 Hardware structure of a 4-point DHT 38 4.20 Architecture of a single butterfly unit FHT processor 41 4.21 Radix-2 FHT butterfly unit hardware structure 41 4.22 Single butterfly unit FHT hardware structure 42 4.23 Data sequence entering the radix-2 FHT butterfly unit 43 4.24 Timing diagram for the single butterfly unit processor 44 4.25 Pipeline architecture of a 16-point FHT 46 4.26 Pipeline hardware structure of a 16-point FHT 46 4.27 Hardware structure of the butterfly unit in the combined last two stages. 47 4.28 Data sequence for each butterfly unit 48 4.29 Timing diagram of the pipeline FHT processor 48 4.30 Superparallel hardware structure of a 16-point FHT 50 4.31 Hardware structure of butterfly unit BUO 50 4.32 Pre-permutation 16-point radix-2 FHT data flow diagram 51 4.33 Data sequence entering the pre-permutation algorithm 52 vii Acknowledgement I would like to express my gratitude to my supervisor Dr. Malcome Wvong for his valuable suggestion and encouragement. I would like to thank Mr. Robert Ross for his technical advice in the use of the departmental computer facilities. Thanks also to Dr. M. Robert Ito for pointing out some interesting properties of the fast Fourier transform, that led to an improvement in this thesis. viii Chapter 1 Introduction 1.1 Digital signal processing Along with the progress of technology, the rapid development of the digital computer has made it a very powerful tool for computation and analysis. Digital signal processing technology takes advantage of the tremendous power of the digital computer to provide much faster and more accurate solutions to scientists and engineers than in the past. Much tedious work in the fields of engineering, biomedicine, physics, astronomy, etc. previously based on the analogue signal processing technology, now can be easily handled with digital signal processing techniques. Digital signal processing is usually carried out for a mixture of data in both the time and frequency domains. The time domain allows processing in real-time systems, while the frequency domain facilitates filtering and spectrum analysis. During regular digital signal processing, data are usually transformed between these two domains frequently. In fact, the transformation is very time consuming, the throughput of most processes depends heavily on the efficiency of the way that the transformation is implemented. The fast Fourier transform (FFT) and the fast Hartley transform (FHT) are solutions to this bottleneck. The FFT is a widely used tool in digital signal processing, and provides a fast nu-merical method to convert data from time domain to frequency domain, and vice versa. The classical Fourier transform is designed for the general case when the data sequence 1 Chapter 1. Introduction 2 is complex, that is, it has real and imaginary parts. For the case when the data sequence is real, its Fourier transform is of complex conjugate nature, this redundancy has been taken advantage of, and several FFT algorithms which make use of this redundancy have been devised [l]-[5], resulting in computer time savings of approximately fifty percent. The FHT is an alternative to the FFT for real-valued sequences, and is well-known for the characteristics that, it has only real arithmetic, and the forms of its forward and inverse transforms are identical. The FHT has a speedup factor of two over the complex FFT [6]-[17], in the case when the complex FFT processes a real sequence as complex by adding an redundant zero-valued imaginary part. The FHT is useful in real data processing in different areas, for example, in one-dimension: spectral analysis, in two-dimension: image processing, and also in multi-dimensional digital signal processsing. 1.2 O v e r v i e w o f t h e f o l l o w i n g c h a p t e r s This thesis presents some software and hardware implementations of the one-dimensional FHT. Chapter 2 gives an introduction to the FHT. The radix-2 and split-radix FHT algo-rithms are implemented in software using several approaches. The FHT software imple-mentaions are evaluated and compared with the complex and real FFT's. Application of the FHT is shown through an image processing example, where an image is first degraded and then restored. Some properties of the FHT are discussed. Chapter 3 discusses the scrambled output problem inherent in the fast radix algo-rithms. A new bit-reversal algorithm is proposed and its efficiency is compared with several other bit-reversal algorithms. Chapter 4 shows some hardware implementations of the DHT and the FHT. Archi-tectures such as single butterfly unit, pipelining and superparallelism are compared. The Chapter 1. Introduction 3 data processing rates are analyzed. Chapter 5 draws conclusions on the implementation results. Chapter 2 The fast Hartley transform (FHT) software implementation 2.1 The Hartley transform In 1942 R.V.L. Hartley introduced a pair of real integral transforms [18] as an alternative to the Fourier transform pair. This pair of formulas have strict reciprocal characteristics for converting signals between time and frequency domains, and are known as the Hartley transform and its inverse. If t is a time variable and / is a frequency variable, the Hartley transform1 [6] for a real function h(t) is /CO h(t) ens (2w ft) dt (2.1) -co and the inverse Hartley transform is /CO H(f)cas(2nft)df (2.2) -co where cas 0 = cos 9 + sin 6. 2.2 The discrete Hartley transform (DHT) In 1983 R.N. Bracewell introduced the discrete Hartley transform (DHT) [19]. The DHT uses discrete time and frequency variables instead of continuous variables as in the Hartley transform, and thus allows numerical analysis on digital computers. 1The original l/y/2n factors of the integrals in Hartley's definition are omitted, when f requency / is used as the variable instead of angular frequency u. 4 Chapter 2. The fast Hartley transform (FHT) software implementation 5 The DHT of a real data sequence h(t) of size iV is defined as N-l H{f) = £ h{t) cas (2TTft/N) (2.3) t=o where / = 0,..., N - 1. The inverse DHT is defined as M O 4 ^ W ) c a s W W ( 2- 4) where * = 0,... ,iV - 1. 2.3 The F H T The DHT requires N2 multiplications and N(N — 1) additions for an iV-element data sequence, then the total arithmetic operations equals 2N2 — N or of order 0(N2). The number of arithmetic operations which grows at a rate of N2 is highly undesirable for the digital computer, and led to the development of the FHT. The FHT is a set of techniques (algorithms) used on digital computers to improve the efficiency of the DHT. In 1984 Bracewell introduced the decimation-in-time radix-2 FHT [7] of order 0(N \og2 JV), and soon after that, many other FHT's of the same order have been pre-sented, for example, Cooley-Tukey radix-2(decimation-in-frequency), radix-4, radix-8, split-radix, prime factor, Winograd, Walsh-Hadamard, etc. [9], [10], [20]-[24]. In fact, most of the fast algorithms applied to the discrete Fourier transform (DFT) basically can also be applied to the DHT because of their similarities. 2.3.1 The radix-2 fast algorithm One of the earliest and most popularly used fast algorithm is the Cooley-Tukey radix-2 algorithm. The radix-2 algorithm has a simple flow structure that can be easily imple-mented in software. Chapter 2. The fast Hartley transform (FHT) software implementation 6 The radix-2 FHT is listed in the following, with the general notations n and k instead of t and / to eliminate the specific reference to time and frequency, since the same FHT algorithm can be applied to both forward and inverse transforms, with only a scaling factor difference between them. The radix-2 decimation-in-time FHT [7] is H(k) = H2n(k) + H2n+1(k) cos (^fc) + H2n+1(N - k) sin (-^fc) (2.5) where k = 0,..., N — 1, and H2n is the DHT of the even-numbered data and H2n+\ is the DHT of the odd-numbered data. The radix-2 decimation-in-frequency FHT [10] is 2 r 'N M /27T, H(2k) = £ [*(") + x ( y + n ) ] c a s { ^ 2 k n ) ( 2 - 6 ) H(2k + 1) = E {[*(») " x ( Y + n)] c o s (^ n ) + Kf " n) " x( iV " n)]sin (%n)}cas (ji2kn) (2-7) where k = 0,..., -y — 1. The radix-2 decimation-in-time and decimation-in-frequency approaches are distin-guished according to the way they process the data sequence, either by doing the mul-tiplication of the trigonometric coefficients firstly or lastly in the butterfly2 as shown in Figure 2.1. However, these approaches are symmetric in flow structures and require the same number of arithmetic operations [25] [26], where the number of nontrivial real ad-ditions, A, and the number of nontrivial real multiplicatons, M, required for an TV-point radix-2 FHT [10] are: A = (Z/2)N\og2N- (3/2)JV + 2 (2.8) M = N\og2N -W + ± (2.9) 2 A butterfly is a basic computation unit in the FHT or the FFT flow structure, such that the implementation of the fast algorithm can be simplified by using the butterfly repetitively. Chapter 2. The fast Hartley transform (FHT) software implementation 7 cos n0— cos n9 — decimation-in-frequency decimation-in-time 9 = 2 J t /N Figure 2.1: Radix-2 decimation-in-frequency and decimation-in-time butterflies. For simplicity, only the decimation-in-frequency approach is discussed in the following. The operation of the radix-2 decimation-in-frequency FHT is first to split the input data into two FHT blocks, where one contains the even-numbered points and another contains the odd-numbered points. Similarly, these two FHT blocks are then split into four FHT blocks, and repeatedly until one-point FHT blocks are obtained. A decimation-in-frequency radix-2 fast Hartley transform has the flow diagram as shown in Figure 2.2 and Figure 2.3. Figure 2.2 shows the first stage of a 16-point FHT decomposed into two 8-point FHT blocks. Figure 2.3 shows the second, third and fourth stages of decomposition in one of the 8-point FHT block. The total number of stages required for the radix-2 FHT equals (log2 N), such that, a 16-point FHT requires 4 stages, a 1024-point FHT requires 10 stages, etc. This decomposition property of the radix-2 FHT is shown clearly in Figure 2.4. Software implementation of the FHT in a general purpose computer has an efficiency closely related to the fast algorithm and implementation technique. The CPU time of a FHT processing is mainly spent on the multiplications and additions of floating-point numbers, and also the computation of the trigonometric functions, the former is chiefly determined by the fast algorithm used, for instant, the split-radix FHT (will be discussed Chapter 2. The fast Hartley transform (FHT) software implementation 8 0 = 2K/ 16 Figure 2.2: Flow diagram of radix-2 decimation-in-frequency decomposition from a 16-point FHT into two 8-point FHT blocks. 6 = 271/8 Figure 2.3: Flow diagram of an 8-point radix-2 decimation-in-frequency FHT. Chapter 2. The fast Hartley transform (FHT) software implementation 9 0 1 2 3 4 5 2 8 2 2 1 6 2 7, 8 2 ? 32 2 7, 8 2 7, 1 6 2 2 8 2 2 6 4 2 ?, 8 2 2 1 6 2 ? 8 2 7 32 4 2 2 8 2 7, 1 6 7, 2 8 2 2 4 2 Figure 2.4: Decomposition pattern of a 64-point radix-2 FHT. in the next subsection) requires less floating-point multiplications and additions than the radix-2 FHT, the latter, is determined by the implementation technique. Some efficient and interesting methods for implementing the radix-2 decimation-in-frequency FHT are presented in the following. The first method (as shown in Appendix A) uses three loops to implement the recur-sive property of the radix-2 algorithm, and is thus termed the "three-loop method". The first loop jumps by stage, and it loops through only (log2 N — 2) stages, because the last two stages are not involved in the multiplication of the trigonometric functions, it is more efficient to treat them separately with special coding. The second loop jumps by data point, and the trigonometric function values are computed in this loop. Since within the same stage, the trigonometric function values required in different blocks are identical, so that, when the third loop jumps by block, the same set of trigonometric function values can be reused in each block, and thus increases Chapter 2. The fast Hartley transform (FHT) software implementation 10 the processing speed. Inside the third loop, a four-point butterfly is used to process four data points at a time. A separate part for sequence with size less than four is added at the beginning of the program to complement the four-point butterfly. The second method (also shown in Appendix A) is an improved version of the first method, and the repetitive computation of the identical trigonometric functions in differ-ent stages is avoided. This is done by first creating a trigonometric function table to store all the precalulated sine and cosine function values, so that later on, when the program requires a trigonometric function value, it will look up the table instead of calculating it repeatedly in different stages. This method is termed the "table-lookup method". One drawback of this method is that a memory space of (N/2) words is required for the trigonometric function table. However, for a process that does fixed length FHT many times, the table needs to be created only once, and the overhead to create the table will become negligible. 2.3.2 The split-radix fast algorithm Among the radix algorithms, the split-radix algorithm is the most efficient one. The split-radix algorithm was first designed for the DFT, which combines the advantages of both the radix-2 and radix-4 fast algorithms, and has its even-numbered part performed in radix-2 and odd-numbered part performed in radix-4. The resulting split-radix algorithm has fewer nontrivial real multiplications and additions than the radix-2, radix-4, and even radix-8 algorithms [27] [4]. The split-radix FHT [9] [20] possesses the same advantages of the split-radix FFT. The even-numbered part of the split-radix FHT algorithm is (2-10) Chapter 2. The fast Hartley transform (FHT) software implementation 11 where k — 0,..., y — 1, and the odd-numbered part is H(Ak + 1) = £ { [ ( * ( » ) - * + » ) ) + ( * ( j - n ) - 1 ( x ~ n ) ) l c o s ( I " ) " [(* ( T + N ) - x (ir+N)) - (* (T " N ) " - N ) ) ] S I N (TF1) 1 cas (2.11) #(4*:+ 3) = S { [ ( s ( n ) ~ x ( Y + N ) ) " (T " n ) " x ( x " n ) ) ] c o s ( l 3 n ) + [ ( * ( f " n ) " x ( x + n ) ) + ( x ( f " n ) ~ x { N ~ n ) ) ] s i n ( 7 7 3 n ) } •cas (2.12) where = 0,..., — 1. The split-radix FHT has a total number of nontrivial real additions, A, and a total number of nontrivial real multiplicatons, M , for an iV-point FHT [9], such that: A = (2/3)N\og2N-(19/9)N+ 3-r(-l)lo*2N -(1/9) (2.13) M = (4/3)A^log2^-(14/9)A^ + 3 + (-l) l o g 2 7 V-(5/9) (2.14) Figures 2.5 and 2.6 show the flow diagrams of the split-radix FHT. The implemen-tation of the split-radix algorithm is more difficult than the radix-2 algorithm. Since the split-radix algorithm is a combination of the radix-2 and the radix-4 algorithms, in a radix-2 stage, a FHT block is decomposed into two half-sized FHT blocks, but in a radix-4 stage, a FHT block is decomposed into four quarter-sized FHT blocks, that means one radix-4 stage is equivalent to two radix-2 stages. This is the reason why the radix-4 is faster than the radix-2 and why it can only accept data sequence with size equal to multiple of four. Chapter 2. The fast Hartley transform (FHT) software implementation 12 e = 2 7 t / 1 6 Figure 2.5: Flow diagram of split-radix decimation-in-frequency decomposition from a 16-point FHT into one 8-point FHT block and two 4-point FHT blocks. Figure 2.6: Flow diagram of an 8-point split-radix decimation-in-frequency FHT. Chapter 2. The fast Hartley transform (FHT) software implementation 13 1 2 3 4 5 ! 8 4 L 2 _ 1 6 2 2 4 2 32 4 L ~ 2 -8 4 l_2_ 2 2 8 4 2 64 2 2 8 4 2 1 6 2 2 4 2 4 L2_ 8 4 La-1 6 2 2 4 2 4 l_2-Figure 2.7: Decomposition pattern of a 64-point split-radix FHT. In Figure 2.7, we can see that the 64-point FHT has its even part decomposed by the radix-2 into a 32-point FHT block and its odd part decomposed by the radix-4 into two 16-point FHT blocks. This combination of different radices makes the proceeding of the even and odd part decomposition unsymmetrical. So that the "three-loop method" applied to the radix-2 requires some modification in order to fit into the split-radix. The first method (as shown in Appendix B) uses the similar approach as in the "three-loop method" in Appendix A. A special part is used for the last three stages which are not involved in the multiplication of the trigonomertic function, in order to optimize the processing speed. The first loop jumps by stage for a total (log2 N — 3) stages. The second loop jumps by point, and the trigonometric function values are calculated in this loop, such that when the third loop jumps by the FHT blocks, each block in the same stage can reuse the trigonometric function values obtained in the second loop and so reduces the C P U Chapter 2. The fast Hartley transform (FHT) software implementation 14 time. The third loop is more complicated than that in the radix-2 implementation, because of the decomposition of the unsymmetrical even-numbered and odd-numbered parts. An elegant method in [28] is utilized here to achieve the task. Inside the third loop, an eight-point butterfly is used to process the data, therefore, a pre-examination part is presented at the beginning of the program to process the FHT with sequence size less than eight. This method avoids the trigonometric functions being calculated repeatedly in each block in the same stage. A further speed-up can be obtained by using a lookup table (as shown in Appendix B). This is similar to that of the radix-2, where a trigonometric function table is created to store all the precalculated sine and cosine values, and thus avoid the repeated calculation of the identical trigonometric functions in different stages. This "table-lookup method" requires a table size of (N/2) words. 2.4 Comparison with the complex F F T For comparison of the FHT with the complex FFT, the radix-2 complex FFT and the split-radix complex FFT [28] [29] based on the three-loop and table-lookup methods are implemented and as shown in Appendix C and Appendix D 3 , respectively. Table 2.1 lists the number of nontrivial real multiplications and additions required for each algorithm, where for the radix-2 complex FFT, A = 3N\og2N -2N + 2 (2.15) M = 2N\og2N -47V + 4 (2.16) and for the split-radix complex FFT [29] [30], A = (8/3)^log 27Y-(16/9)^ + 2 - ( - l ) l o g 2 ; v . ( 2 / 9 ) (2.17) 3Only the forward transform programs are listed. The inverse transform programs can be obtained easily from the forward transform programs by taking complex conjugate on the twiddle factors. A method to use a single forward FFT for both forward and inverse transforms is shown in [31] [32]. Chapter 2. The fast Hartley transform (FHT) software implementation 15 N FHT complex FFT Radix-2 Split-radix Radix-2 Split-radix adds mults adds mults adds mults adds mults 2 2 0 2 0 4 0 4 0 4 8 0 8 0 18 4 16 0 8 26 4 22 2 58 20 52 8 16 74 20 64 12 162 68 144 32 32 194 68 166 42 418 196 372 104 64 482 196 416 124 1026 516 912 288 128 1154 516 998 330 2434 1284 2164 744 256 2690 1284 2336 828 5634 3076 5008 1824 512 6146 3076 5350 1994 12802 7172 11380 4328 1024 13824 7172 12064 4668 28674 16388 25488 10016 2048 30722 16388 26854 10698 63490 36868 56436 22760 4096 67586 36868 59168 24124 139266 81924 123792 50976 Table 2.1: Number of real arithmetic operations for the fast algorithms. Power of two N FHT complex FFT three-loop table-lookup three-loop table-lookup r-2 s-r r-2 s-r r-2 s-r r-2 s-r 8 256 0.10 0.08 0.08 0.07 0.22 0.20 0.17 0.15 9 512 0.22 0.20 0.17 0.15 0.48 0.42 0.37 0.32 10 1024 0.48 0.42 0.37 0.33 1.02 0.88 0.82 0.68 11 2048 1.03 0.88 0.83 0.70 2.20 1.88 1.77 1.48 12 4096 2.17 1.88 1.78 1.53 4.65 3.97 3.82 3.15 13 8192 4.57 3.95 3.87 3.22 9.73 8.30 8.02 6.68 14 16384 9.65 8.30 8.15 6.82 20.40 17.45 17.02 14.18 15 32768 20.25 17.42 17.35 14.48 42.65 36.43 36.00 29.93 16 65536 42.43 36.43 36.90 30.53 89.23 75.97 75.85 63.17 Table 2.2: Comparison of the fast algorithms in CPU execution time (in seconds) on a Sun 3/60 workstation, where r-2 stands for the radix-2 algorithm and s-r stands for the split-radix algorithm. Chapter 2. The fast Hartley transform (FHT) software implementation 16 o.oiH 1 1 1 1 8 10 12 14 16 Power of two Figure 2.8: FHT requires less CPU time than the complex FFT in the radix-2 software implementation in both the "three-loop (lp)" and "table-lookup (tb)" methods. Power of two Figure 2.9: FHT requires less CPU time than the complex FFT in the split-radix software implementation in both the "three-loop (lp)" and "table-lookup (tb)" methods. Chapter 2. The fast Hartley transform (FHT) software implementation 17 M = (4/3)A^log 27V-(32/9)iV-f4-(-l) l o g 2 i V-(4/9) (2.18) The easiest method to apply the complex FFT on a real sequence, is to expand the real sequence to complex by adding redundant zero-valued imaginary part. Based on this method, the execution time4 of the FHT and the complex FFT of a real data sequence is computed and compared in Table 2.2, the graphs in Figures 2.8 and 2.9 based on Table 2.2 further clarify how the execution time increases with N. One important point not shown in Table 2.2 is that in the table-lookup method, the FHT requires a table size of (N/2), where the complex FFT requires a table size of N, double that of the F H T . It is seen that, under the same fast algorithm and implementation method, the FHT is about twice as fast as the complex FFT. It is obvious that, the complex FFT based on the above method wastes one-half of the CPU time, on computing the redundant zero-valued imaginary data, and is then very inefficient. There are two other commonly used methods [31] to deal with real sequence in the complex FFT, and without the use of the redundant zero-valued imaginary part. The first method is to transform an N-point real sequence by means of an N/2-point complex FFT. First of all, an N-point real sequence is divided into two N/2-point sequences, where one contains the even-numbered data, and another contains the odd-numbered data of the original sequence. These two N/2-point sequences are input as the real and imaginary parts to an N/2-point complex FFT. The transform of the N-point sequence can be obtained by combining the real and imaginary parts of the complex F F T output with some algebraic computation. The second method allows the transform of two N-point real sequences simultaneously by using one N-point complex FFT. In this case, the first real sequence is input as the real part and the second sequence is input as the imaginary part of an N-point complex F F T . 4Does not include the input, output and bit-reversing time. Chapter 2. The fast Hartley transform (FHT) software implementation 18 The transforms of these two N-point sequences can be decomposed from the N-point complex FFT output with some additional computation. These two methods have the similar speedup factor as the FHT. However, as indicated in [13] [16], these methods cannot be used to retransform their own output, since they take only real data. Some modified approaches should be used for the inverse transform to convert the complex conjugate symmetric sequence back to a real sequence. So that the FHT is more convenient here for real sequences, since the identical forward transform program can be used for the inverse transform with only a scaling, and does not require extra computation for combining or sorting out the result. 2.5 Comparison with the real FFT Some specially designed FFT's for real sequence [l]-[5], utilize the redundancy and sym-metry of the complex FFT's to reduce the number of multiplications by exactly one-half, and the number of additions by more than one-half of that of the complex FFT's. So that, these real FFT's should run at around the same speed as the FHT under the same fast algorithm implementation. Also, the real FFT requires the same size N storage space as the FHT, and because the real FFT and the FHT are very similar, the other overhead costs should be almost the same. The real FFT transforms a real sequence into a complex conjugate symmetric se-quence, the inverse of the real FFT transforms a complex conjugate symmetric sequence back into a real sequence, such that, two algorithms are required for the application which needs both forward and inverse transforms. The FHT is then more convenient than the real FFT, since only one algorithm is sufficient for both forward and inverse transforms. Chapter 2. The fast Hartley transform (FHT) software implementation 19 2.6 To get the F F T from the F H T The Hartley transform, H(f), can be split into its even part, Heven(f) and odd part, Hodd(f) [6], such that H(f) = Heven(f) + Hodd(f) (2.19) The even part is defined as HeUf) = \[H(f) + H(-f)} (2-20) /CO x(t)cos(2irft)dt (2.21) -oo and the odd part is defined as Hodd{f) = \\HU)-H(-f)] ' (2.22) /CO x(t) sin(2n ft) dt (2.23) -OO The two integrals in Equation (2.21) and (2.23) are known as the Fourier cosine transform and the Fourier sine transform respectively. It is apparent that /oo x(t)[cx>s(2vft) _ j sin(27r/0] dt (2.24) -co /oo x(t) t-^}tdt (2.25) -co is in fact the Fourier transform, F(f). That is F(f) = Heven(f) - JHodd(f) (2.26) In other words, Heven(f) is equal to the real part, Freai(f) of the Fourier transform and Hodd(f) 1 S equal to the sign-reversed imaginary part, FimagU) of the Fourier transform, then the Hartley transform can be written as H{f) = FrealU) ' ^ U ) (2-27) Chapter 2. The fast Hartley transform (FHT) software implementation 20 Equations (2.26) and (2.27) show that from either the Fourier transform or the Hartley transform, one can easily obtain the other transform. Similarly, in the discrete case, for a data sequence of N elements, the FFT of a real function can be obtained from the even and odd part of the FHT. Since H(f) is considered periodic with period N, then Heven(f) = ±[H(f) + H(N - f)] (2.28) HouU) = \[H{f)-H(N-f)} (2.29) and from 2.26 F(f) = Heven(f)-JHodd(f) H(N - f) + H(f) , H(N - /) - H(f) = 2 + J 2 ( 2 - 3 0 ) where / = 0,..., N - 1, and H(N) = H(0). However, Equation (2.30) is only appropriate for real-input data, since neither H(N — f) nor H(f) can deal with complex numbers, so that, a method to obtain the complex FFT by using the FHT is developed in the following. Suppose y(t) is a complex function such that y(t) = a(i) + j b(t), where a(t) and b(t) are real functions, then the FFT of y(t) is Yp(f). By using the linearity of the FFT, the FFT of a(t) + j b(t) is Ap(f) + j BF(f), where AF(f) and BF(f) are complex. If AF(f) and BF(f) are expressed in their real and imaginary parts, then y>(/) = AF(f) + j BF(f) = [AFr(f) -f ; AF,(f)} + j [BFr(f) + j BFi(f)\ = [AFr(f)-BFi(f)}+J[AFi(f) + BFr(f)} (2.31) Chapter 2. The fast Hartley transform (FHT) software implementation 21 since AFr{f) = - f) + AH(f)} (2.32) AFi(f) = \\AH{N - f) - AH(f)\ (2.33) where AH(f) is the FHT of a(t), and BFr(f) = \[B„(N - f) + BH(f)} (2.34) BFi(f) = \[BH(N - f) - BH(f)} (2.35) where BH(f) is the FHT of b(t), Equation (2.31) becomes YF(f) = \[AH(N - f) + A„(f) - BH(N - f) + BHU)\ + j l-[AH(N - f) - AH(f) + BJJ(N -f) + Bfj(f)] (2.36) With Equation (2.36), one can obtain the FFT of complex-input data through the FHT. 2.7 F H T convolution and power spectrum Two of the most common applications of the FHT is for convolution and to obtain power spectrum for filtering. The FHT convolution theorem [7] states that if h(t) is the convolution of two real functions, hi(t) and h2(t), that is h(t) = hi(t) * h2(t) = ^ h^t') h2(t - t') (2.37) t'=o then H(f) = Hx(f)H2e(f) + H,(N - f)H2o(f) (2.38) where H(f), H\(f) and H2(f) are the FHT's of h(t), hi(t) and h2(t), respectively, and H2(f) = H2e(f) + if 2 o(/), the sum of its even and odd parts. Chapter 2. The fast Hartley transform (FHT) software implementation 22 From Equations (2.28) and (2.29), Equation (2.38) becomes H(f) = \[Hx{f)H*U) + Hx{f)H*{N - f) + Ht(N- f)H2(f) - Hi(N - f)H2(N - /)] (2.39) However, when either h\(t) or h2(t) is symmetrical (even function), that is, either /ii(£) = hi(N — t) or h2(t) = h2(N — t), then the convolution theorem in Equation (2.39) can be simplified to a single real multiplication of the two FHT's: H(f) = H1(f)H2(f) (2.40) The power spectrum, Pa(f), of a real function h(t), calculated with the FHT is Ps(f) = \H(f)\2 (2.41) = \{[H(f)Y + [H(N - /)] 2} (2.42) where H(f) is the FHT of h(t). 2.8 A n image processing example of the application of the F H T An application of the FHT is presented in this section through an image processing example. In this example, an image is first degraded, and then restored. The theorems for addition, convolution, power spectrum, etc. on the FHT are applied. Since the image is a two-dimensional signal, in order to use the fast one-dimensional algorithm in Appendix B, the input image is transformed into one dimension by joining each row of the two-dimensional image matrix serially with the previous row5. This is similar to T V image processing. 5This method suffices only in this example. The general methods to obtain the two-dimensional F H T with the one-dimensional FHT algorithm are derived in [13]. A two-dimensional FHT algorithm is developed in [33]. Chapter 2. The fast Hartley transform (FHT) software implementation 23 Degradation Restoration H g g w r n Figure 2.10: Image degradation and restoration processes. The image processing is carried out by passing the image signal through a linear-motion blur operator and with the addition of some white noise, the degraded image signal is then restored by a Wiener filter. The processes are as shown in Figure 2.10. If an image x(t), passes through a linear-motion blur operator m(i) and interfered by a white noise n(t), the degraded image g(t) can be expressed as: g(t) = x(t) * m(t) + n(t) (2.43) where "*" is the convolution operator. A simple method to restore the image is to pass the degraded image through a Wiener filter. The restored image r(t) can be expressed as: r(t) = w(t)*g(t) (2.44) where w(t) is the Wiener filter function in time domain. The Wiener filter function in frequency domain, W(f) is M(N - f) W(f) = (2.45) \M(f)\* + NSR where M(f) is the linear-motion blur operator in frequency domain, and NSR is the noise to signal ratio, which is defined as the ratio of the power spectrum of noise to the power spectrum of signal. Chapter 2. The fast Hartley transform (FHT) soft ware implementation 24 Figure 2.12: Degraded image. Chapter 2. The fast Hartley transform (FHT) software implementation 25 Figure 2.13: Restored image. The degradation starts with a portrait image as shown in Figure 2.11. Each pixel in the image has a gray level between 0 to 255. A pixel in gray level 0 is black, in gray level 255 is white, and in a gray level between 0 and 255 has different gray scale. The image is first input to a linear-motion blur operator, which stretches each pixel in the image to a distance of nine pixels to the right, then the blurred image is further degraded with white noise of random gray level from 0 to 10. The degraded image is shown in Figure 2.12. The degraded image is restored with a Wiener filter, which does the inverse of the blur operator and compensates the noise level. The restored images as shown in Figure 2.13 has a root-mean-square error of 8.4 gray levels with respect to the original image. Chapter 2. The fast Hartley transform (FHT) software implementation 26 2.0 Summary The radix-2 and split-radix FHT's are implemented in the "three-loop method" and "table-lookup method", and compared with the FFT on real sequence. It is found that although the FFT can be as fast as the FHT, the FHT is more convenient in performing the inverse transform on the complex conjugate symmetric data sequence back to a real sequence. The methods to obtain the FFT from the FHT and vice versa are discussed. The applications of the FHT in convolution and computation of power spectra are shown through an image processing example. Chapter 3 The bit-reversal algorithms In the radix FHT algorithms, the data after processing are stored back in the original place where the data were obtained from, in order to economize on data memory usage. However, this "in-place" storing method generates a transform with output data sequence in bit-reversed order. For the applications which require the transform sequence in the correct order, an extra permutation stage should be used to unscramble the data. For the decimation-in-time approach, the permutation stage is placed in front of the transform, and is called the "pre-permutation". For the decimation-in-frequency approach, the permutation stage is placed after the transform, and is called the "post-permuation". However, the permutation stages in both approaches are exactly the same and can be performed with the bit-reversal algorithm. 3.1 Some existing algorithms One of the most widely used bit-reversal algorithm is introduced by Gold and Rader. An implementation of this algorithm [32] is repeated in pseudo-code form in Figure 3.14. This algorithm uses a "bit-reversed counter" to generate bit-reversed index sequence which contains the reversed bits of the index sequence from zero to (TV—2). A comparsion is made between each index and the corresponding bit-reversed index generated, and a swap of the values in these two indices is required whenever the bit-reversed index is greater than the original index in value. The swapping required for a 16-point data sequence is shown in Table 3.3. 27 Chapter 3. The bit-reversal algorithms subroutine reverse(x[ ] , N) reversed.index = 0 N2 = N / 2 for index = 0 to N - 2 i f (index < reversed.index) then swap(x[index], x[reversed.index]) k • N2 while (k <= reversed.index) reversed.index = reversed.index - k k = k / 2 end while reversed.index = reversed.index + k next index end subroutine Figure 3.14: Gold-Rader's bit-reversal algorithm (after Rabiner and Gold). Index in binary reversed swap 0 0000 0000 1 0001 1000 YES 2 0010 0100 YES 3 0011 1100 YES 4 0100 0010 5 0101 1010 YES 6 0110 0110 7 0111 1110 YES 8 1000 0001 9 1001 1001 10 1010 0101 11 1011 1101 YES 12 1100 0011 13 1101 1011 14 1110 0111 15 1111 1111 Table 3.3: Bit-reversing of a 16-point array. Chapter 3. The bit-reversal algorithms 29 The operation of the bit-reversed counter starts from bit-reversed index zero, the counter generates the next bit-reversed index by incrementing one in the most significant bit of the present bit-reversed index, and a rightward carry is carried out if necessary. This procedure is repeated in the counter until the second last bit-reversed index is obtained, where the last bit-reversed index is not required to be generated simply because no swapping is involved1. To increment by one in the most significant bit, means adding a value equivalent to (N/2), where carry is added rightwardly. A rightward carry can be obtained by checking for the first zero bit starting from the most significant bit. If a "one" is found then change it to zero, and check the next bit on the right hand side until a zero is found, then change the zero to one. To change a one to a zero can be done by substracting (N/2), or (N/4), or (N/8), etc. for each one bit that is found. Similarly, a zero can be changed to a one by adding (N/2), or (N/4), or (N/8), etc. Since the last index which has all one-bits need not be checked, there must be at least one zero bit in each index, so that overflow will not occur. Buneman [8] shows an equivalent bit-reversal algorithm but with a lookup table to store the (N/2), (N/4), (N/8), etc. values to reduce the CPU time spent on the repeated integer divisions. Since in some high-level computer languages, the integer division cannot be replaced by simple right-bitwise shifting as in the assembly language, instead, the integer division is performed in time consuming long division. The Buneman's algorithm is shown in pseudo-code form in Figure 3.15. Rodriguez [34] presents an improvement to the Gold and Rader's algorithm by setting an upper bound on the index sequence so as to eliminate the CPU time on checking the indices beyond the upper bound, where the upper bound means the last index of swapping, which has a value equal to (TV — 1 — s/R * N), where R is "one" if the power 1In fact, if swapping is determined under the condition that the bit-reversed index is greater than the original index, then the second last bit-reversed index needs not be generated too, because the second last index must have a value greater than its reversed bits. Chapter 3. The bit-reversal algorithms subroutine reverse(x[ ], power, N) P = i ; for index = 0 to power table[index] = p p = p + p next index reversed.index = 0 for index = 1 to N - 2 p = power do { p = p - 1 reversed.index = reversed.index - table[p] } while (reversed.index >= 0) reversed.index = reversed.index + table[p + 1] i f (index > reversed.index) then swap(x[index], x[reversed.index]) next index end subroutine Figure 3.15: The table-lookup bit-reversal algorithm (after Buneman). subroutine reverse(x[ ], power, N) upper.bound = N - 1 - 2 ** floor{(power + 1) / 2} /* floor{s} = the greatest integer smaller or equal to s */ reversed.index = 0 N2 = N / 2 for index = 1 to upper.bound k = N2 while (k <= reversed.index) reversed.index = reversed.index - k k = k / 2 end while reversed.index = reversed.index + k i f (index < reversed.index) swap(x[index], x[reversed.index]) next index end subroutine Figure 3.16: Improved Gold-Rader's algorithm (after Rodriguez). Chapter 3. The bit-reversal algorithms 31 of two is even or "two" if the power of two is odd. The Rodriguez's improved algorithm is shown in pseudo-code form in Figure 3.16. Evans [35] introduces a more efficient but also more complex bit-reversal algorithm than all the algorithms shown above. This algorithm uses a precalculated seed table to enhance the speed in generating the bit-reversed sequence, and employs some sophisti-cated rules to govern the swapping. An implementation of this algorithm is shown in Appendix E . l . Additional work by Evans in [36] claims that with the addition of an extra rule, this algorithm can run slightly faster in most computers. 3.2 A new bit-reversal algorithm Gold and Rader designed the bit-reversal algorithm based on the bit-reversed counter some two decades ago. The Gold-Rader's bit-reversed counter utilizes integer division on the operations of incrementation and rightward carry, since probably at that time, no suitable high-level language could be used to handle the bitwise operations, and the assembly language programming is tedious. By taking advantage of the bitwise operators of the C language, a new bit-reversed counter is developed. Table 3.4 shows the operation of a bit-reversed counter of N — 16 based on the use of the bitwise "exclusive-or" operator. By inspection, one can see that the reversed bits of the odd-numbered indices can be generated by exclusive-oring the reversed bits of the previous even-numbered indices with the operand "1000", whereas the reversed bits of the even-numbered indices can be obtained by exclusive-oring the reversed-bits of the previous odd-numbered indices with the operand "1100", "1110" and "1111", respectively. The rule for choosing the appropriate "exclusive-or" operand is that, always choose Chapter 3. The bit-reversal algorithms 32 Index in binary reversed xor operand xor result 0 0000 0000 1000 1000 1 0001 1000 1100 0100 2 0010 0100 1000 1100 3 0011 1100 1110 0010 4 0100 0010 1000 1010 5 0101 1010 1100 0110 6 0110 0110 1000 1110 7 0111 1110 1111 0001 8 1000 0001 1000 1100 9 1001 1001 1100 0101 10 1010 0101 1000 1101 11 1011 1101 1110 0011 12 1100 0011 1000 1011 13 1101 1011 1100 0111 14 1110 0111 1000 1111 15 1111 1111 Table 3.4: Exclusive-oring a bit-reversed index in column three and the related operand in column four produces the next bit-reversed index in the data sequence. the operand with width equals log2 TV, and with one or contiguous "1" bits in the most significant bits and some contiguous or no "0" bit in the least significant bits, such that the operand has a value just greater than the value of the bit-reversed index. For instant, in the 16-point sequence example, the operands to be chosen are in the format: "1000", "1100", "1110" and "1111", where the operand to be used for "1001", the reversed bits of index 9, in order to generate the reversed bits of index 10 is "1100". The operand required for "1101", the reversed bits of index 11, in order to gernerate the reversed bits of index 12 is "1110". An alternative rule governing the use of the operands is to use the operand "1000", starting from index (2° — 1) and every (21)-index times, the operand "1100" is used every (22)-index times starting from index (21 — 1), the operand "1110" is used every (23)-index times starting from index (22 — 1), and the operand "1111" is used only once starting Chapter 3. The bit-reversal algorithms 33 subroutine reverse(x[ ], N) lst.operand = N / 2 2nd_operand = (lst.operand » 1) ft (N / 2) reversed.index = N / 2 for index = 1 to N - 2 /* swapping for the odd-numbered index */ if (index < reversed.index) then swap(x[index], x[reversed_index]) /* generate reversed bits for the next even-numbered index */ k = 2nd_operand while (k <= reversed_index) k = (k » 1) $ (N / 2) /* next higher operand */ end while reversed.index = reversed.index " k index = index + 1 /* swapping for the even-numbered index */ i f (index < reversed.index) then swap(x[index], x[reversed.index]) /* generate reversed bits for the next odd-numbered index */ reversed.index = reversed.index " lst.operand next index end subroutine Figure 3.17: The new bit-reversal algorithm. from index (23 — 1). In general, the operand with width equals log2 N and n "1" bits in the most significant bits, where n is from 1 to log2 N, and "0" in the other bits, is used firstly at index position (2 n - 1 — 1) and repeated every (2")-index times. The operands to be used for data sequence in other length can be predicted in a similar way by using either one of the rules above. Since the implementation of the second rule may employ more complicated looping technique and hence more CPU time, the first rule is chosen for implementation, and is presented in pseudo-code form in Figure 3.17. The technique to separately process the even-numbered indices and the odd-numbered indices has a significant speedup effect. Chapter 3. The bit-reversal algorithms 34 N A B C D E lk 0.008 0.012 0.013 0.017 0.010 2k 0.015 0.025 0.027 0.032 0.020 4k 0.032 0.050 0.052 0.063 0.040 8k 0.063 0.102 0.103 0.128 0.082 16k 0.128 0.203 0.210 0.253 0.163 32k 0.258 0.407 0.418 0.507 0.328 A: Evans' algorithm. B: Gold and Rader's algorithm. C: Rodriguez's algorithm. D: Buneman's algorithm. E: The new algorithm. Table 3.5: CPU execution time (in seconds) of the bit-reversal algorithms running on a Sun 3/60 workstation. A C-language implementation based on this pseudo code is shown in Appendix E.5. A second implementation using a lookup table to store the exclusive-or operands is shown in Appendix E.6. A test by using the UNIX-C compiler and a Sun 3/60 workstation shows that the implementation in Appendix E.5 can run faster, however, in other machines or C-compilers, where the indexing is comparatively quick, the implemenation in Appendix E.6 may be more efficient. 3.3 Comparison of the bit-reversal algorithms Table 3.5 and Figure 3.18 show the time required for different sizes of data for the above five order O(N) algorithms. The result is obtained by implementing the algorithms in the C language (as shown in Appendix E.l-5), and running on a Sun 3/60 workstation equipped with a floating point processor. The new algorithm is one of the fastest, and is only second to Evans' bit-reversal algorithm. However, Evans' algorithm requires a complex implementation and the use of a seed table, so that the new algorithm is very Chapter 3. The bit-reversal algorithms 35 Figure 3.18: Plotting of Table 3.5. attractive. The Gold and Rader's algorithm is found to be simple and efficient, it is suitable for the implementation in the high-level languages lacking the bitwise operators. The algorithms by Buneman and Rodriguez do not show the expected effect, it is because the integer division2 in the Sun workstation with the UNIX-C compiler can run as fast as a right-bitwise shifting operation, the indexing and other overheads in the improvement unintentionally slow down the algorithms. 3.4 Summary A new bit-reversal algorithm is presented. This new algorithm is found to be very efficient, since it is fast, has minimum memory demands and can be implemented easily in the C language. Several other bit-reversal algorithms are introduced and compared with the new algorithm. A fast bit-reversal algorithm is desirable and important in the FHT, because it accounts for a few percent of the total run time of a complete radix 2The integer division is replaced by right-bitwise shifting in the programs in Appendix E, just for a consistent comparison of the execution speed between different algorithms Chapter 3. The bit-reversal algorithms 36 transform. The usage of the bit-reversal algorithm in the hardware implementation of the FHT will be discussed in the next chapter. Chapter 4 The FHT hardware implementation 4.1 Introduction The software implementation of the fast algorithms on a general purpose computer or on a general purpose computer equipped with a programmable digital signal processor (DSP) chip is generally speed-limited for data processing from low to medium sampling rate, say, for a maximum of less than a million samples per second1. For some very high speed processing environment, for example, radar signal and real-time digital image processing, custom-designed hardware implementation of the transform is highly preferred. In this chapter, the hardware implementation of the DHT is proposed, and the FHT hardware implementation based on several advanced architectures is developed and evaluated. 4.2 Hardware implementation of the DHT A four-point DHT based on the expression in Equation 2.3 is implemented and is shown in Figure 4.19. This hardware structure can be easily modified for data sequence of other size by appending or reducing devices in a similar layout pattern. The operation of this hardware implementation is synchronized with clock signals, namely CLK-A, CLK-B, CLK-C and CLK-D. The control clock signals of the system are generated with a ring counter, where, each of the four outputs of the ring counter sends a clock signal 1As published in [37], a benchmark of 2.36ms for a 1024-point complex radix-2 FFT was recorded with the TMS320C30 DSP; this is equivalent to a processing rate of 0.434 million samples per second. The split-radix FHT implemented with the "table-lookup method" in Appendix B, requires 39.4ms for a 1024-point real data sequence on a Sun SPARC station 1 workstation. 37 Chapter 4. The FHT hardware implementation 38 INPUT 1' CLK-A CLK-B CLK-C CLK-D ^ x - ^ - ' ^ x 3: t t l ROM 2 r * X CLK-A 3t L I r-i L 3-i7| I X CLK-B CLK-C , I CLK-D ' 4 x V 4 ' ' ' ' 3-L CLK-A ROM 4 OUTPUT L: latch 3-L: tri-state latch Figure 4.19: Hardware structure of a 4-point DHT. sequentially to trigger a group of devices. The devices to be triggered synchronously, including latches and multipliers, but excluding the adders, are grouped together by the dotted-diagonal lines in Figure 4.19. The input line is a multi-bit data bus. Data are entered through the input line to the system serially at each clock signal. Four input latches controlled by clock signals are connected to the input line to select, store and pass the corresponding data points. In each clock signal, only one input latch can pass in one data point. This data point is retained in the input latch until the next time it is triggered and another new data point is sampled. The multipliers sample data from the input latches and trigonometric Chapter 4. The FHT hardware implementation 39 coefficients from the read-only memory (ROM). The adders are driven by the data from the multipliers, and since the two multipliers feeding data to the adder are tiggered in two consecutive clock signals, a latch is required to retain the output of the first multiplier and let the two data points arrive at the adder simultaneously. After three additions in each column, the transform output is passed to the output line in each clock signal. The system operation starts at CLK-A. The input latch at the upper-left corner in Figure 4.19 passes in the first data point and R0M1 is addressed and provides the cas coefficient; these two data are passed to the top multiplier in the first column. At CLK-B, the multiplier is fired and the product is passed to a latch. At CLK-C the latch is selected and the datum inside is passed to an adder; meanwhile, the adder accepts another datum from the second multiplier in the first column. Since the adder is data driven, the adder outputs the sum to the second latch in column one without waiting for another clock signal. The second latch is then selected by CLK-D and the datum stored inside is passed to the next adder. Similarly, the final sum from the first column is passed to the output line through a tri-state output latch which will be triggered by the corresponding clock signal. Just in front of the last column of multipliers, there are some other latches connected to the input latches. These are used to hold the input data, because when one of the multipliers in the last column is triggered, the input latch connected to that multiplier is being triggered simultaneously and the datum in the input latch will be lost, so that an intermediate-stage latch is required to hold the datum. This DHT implementation requires six clock pulses of overhead, and it can convert one data point to the desired domain in one clock cycle. It has a simple structure and control, and is flexible in expansibility. However, for an TV-point transform, TV2 word-Chapter 4. The FHT hardware implementation 40 sized ROM, N(N — 1) adders and N2 multipliers2 are required. For a large N, the device count is very high, and turns out to be a high hardware implementation cost. In addition, each group of devices are triggered only in the corresponding clock signal, which means a low efficiency in the use of the devices. As an alternative, hardware implementation using the fast algorithm is then desired. Several FHT hardware implementation approaches are shown in the following section. 4.3 Hardware implementation of the FHT Unlike the FFT algorithms, most of the FHT algorithms, including the radix algorithms, have comparatively irregular geometry in their data-flow structure. Figures 2.2 and 2.5 are good examples to illustrate this point. The conventional FFT hardware implementation methods [38]-[40] utilize the regular geometry of the FFT; when these methods are applied to the FHT, more complex control and data sequence rearranging are required. This is probably the reason why the FHT hardware implementation is still not so popular today. In this section, the radix-2 algorithm is used as an example of implementation because of its popularity and simplicity. Higher radix algorithm implementation is possible with similar technique, but a more complicated control circuit is required. 4.3.1 Single butterfly unit FHT architecture Figure 4.20 shows a simple single butterfly unit architecture of the radix-2 FHT processor. The architecture shown can be configured with simple hardware. There is only one memory unit, one control unit and one butterfly unit. 2In the actual implementation, the multipliers in the first column and the first row are not required, since a cas value equivalent to one is to be multiplied in these multipliers, data can be passed directly to the latches or adders instead. Chapter 4. The FHT hardware implementation 41 Memory •4 input Unit • output (FHT) Control i i Unit < Butterfly Unit Figure 4.20: Architecture of a single butterfly unit FHT processor. input registers output registers cos(nG) tin(nO) Figure 4.21: Radix-2 FHT butterfly unit hardware structure. The memory unit stores the input data, the intermediate results, the transform output and the trigonometric coefficients. The control unit generates memory addresses, memory clock cycles, butterfly unit clock cycles and other control signals. The butterfly unit performs real additions, substractions and multiplication. Figure 4.21 shows the hardware implementation of the 16-point radix-2 decimation-in-frequency FHT butterfly unit based on the butterfly flow structure in Figure 2.1. The radix-2 FHT butterfly unit can process four real data in each computation cycle, Chapter 4. The FHT hardware implementation 42 memory array 1 memory array 3 input output memory array 2 cosine and sine coefficients Figure 4.22: Single butterfly unit FHT hardware structure. Figure 4.22 shows the hardware inplementation of the single butterfly unit architec-ture. The operation of this implemenation can be explained with the 16-point radix-2 FHT in Figure 2.2 as an example. The 16-point transform has three sets of data points in its first stage, namely, the first set: points 1, 7, 9 and 15, the second set: points 2, 6, 10 and 14, and the third set: points 3, 5, 11 and 13, which are ready to be input to the buffer registers, a, b, c and d of the butterfly unit as in Figure 4.21. The remaining points, 0, 4, 8 and 12 do not require a complete operation of the butterfly unit, that is, no multiplication is required. However, it is not economical in this architecture to provide a separate butterfly unit for the processing of this special case. A solution to this is by entering these points to the butterfly unit also as a set, and inputting a one or cos(0), and a zero or sin(0) from the ROM to the multipliers, the multiplications can virtually be eliminated. Similarly, the data points in the last two stages that do not require complete operation of the butterfly unit, can be processed with the same method to eliminate the multiplication. Figure 4.23 shows the detail of the data point sequence entering the butterfly unit. The corresponding trigonometric coefficients used in the multiplications are listed along with the data sequence, where the data sets without trigonometric coefficients shown will be multiplied by cos(0) and sin(0). Chapter 4. The FHT hardware implementation 43 stage 1 stage 2 stage 3 stage 4 d: a: b c: 0 1 2 3 4 7 6 5 8 9 1011 12 15 14 13 0 1 8 9 2 3 10 11 4 5 12 13 6 7 14 15 0 4 8 12 1 5 9 13 2 6 10 14 3 7 11 15 0 4 8 12 2 6 10 14 1 5 9 13 3 7 11 15 cos 6 sin 6 8 = 27C/N cos 28 cos 36 cos 8 cos 9 sin 26 sin 36 sin 8 sin 6 Figure 4.23: Data sequence entering the radix-2 FHT butterfly unit. The data sequence in each stage can be fetched with the addresses generated by a pre-defined sequence counter. Therefore, for the 16-point FHT example, four counters are required to work sequentially for each transform. The memory unit contains three memory arrays. Memory array 1 and 2 are input buffer memory, while memory array 3 is an output buffer memory. The transform process starts with serial data input to memory array 1. In the example, sixteen data points for the first transform will be buffered. If there is more than one transform, an input switching to memory array 2 will occur, and the next sixteen data points for the second transform will enter to memory array 2. When memory array 2 is busy in receiving the serial input data, the data already stored in memory array 1 are read into the buffer registers of the butterfly unit. The computation output of the butterfly unit is written back to memory array 1 in-place, except that the final transform result from the last stage is written to memory array 3. The output buffer memory is required to provide enough time for a slower data transferring to the next processing stage, and to allow the serial input data of the next transform to enter memory array 1 immediately in order to keep the butterfly unit running uninterruptedly. The butterfly unit requires (^log2N) butterfly cycles for an TV-point transform. A Chapter 4. The FHT hardware implementation 44 Butterfly Unit: | B1 | B2 | B3 | B4 | B5 , B6 | B7 | B8 i • • • • Memory Array 1: i R l IR IW, R , W, R , W ,R , W ,R , W ,R , W ,R , W ,R , . . . . Memory Array 2: | IN | IN | IN i IN | IN i IN | IN | IN | IN . . . .  Memory Array 3: * * * ' Butterfly Unit: " ' ' I B 12 | B13 | B 14 | B 15 | B 16 | B 1 | B2 | B3 | B4 , BS . . . . Memory Array 1: ' ' ' |W|R |W|R | |R i |R | | IN | IN | IN , IN | IN , IN • . . • Memory Array 2: i IN , IN , IN , IN , R , ,R ,W,R , W,R , W,R , W,R • • • •  Memory Array 3: ' U l w l l w l l w l l w l . . . . R: 4 read cycles, W: 4 write cycles, B: butterfly cycle, IN: 1 input cycle. Figure 4.24: Timing diagram for the single butterfly unit processor. . 16-point transform which has four stages then requires sixteen passes in the butterfly unit. Figure 4.24 shows the timing cycles starting from memory array 2 begins to receive the second transform data. One can see that, the use of the overlapped memory and butterfly cycle, and the use of the output memory buffer increase the processing throughput and enhance the butterfly unit to run at full efficiency. The control unit generates read and write memory cycles, butterfly cycles, device select control signals, read memory addresses, etc. Since the read memory addresses are not in sequential order, custom-designed counters are required. Figure 4.24 shows that, when a datum is read into the buffer register a in Figure 4.21 in a read memory cycle, it will wait for the equivalence of twelve memory cycles, and at the thirteenth memory cycle, the output at buffer register a! is written back to the input buffer memory in-place. That is, two sets of data points are read into the butterfly unit before one set of output is written back to the input buffer memory. Therefore, in order Chapter 4. The FHT hardware implementation 45 to avoid the extra hardware for generating the write addresses, eight address registers are used in the control unit to retain the read addresses for twelve memory cycles. The processing rate of this hardware structure depends on the memory access time and the processing time of the butterfly unit. In this design, the butterfly unit processing time matches the time for four read and four write cycles to optimize the efficiency. Suppose, if memory with access time of 60ns is used, then the butterfly unit can be designed to a maximum procesing time of 480ns or vice versa. The procesing time of a transform with length N, can be expressed as 8r(-j) log2 N, where r is the memory access time. Therefore, if memory of 60ns is used, a 16-point transform will take 7.68/JS, or a processing rate of 2.08 million samples per second. If N = 1024, the processing time for one transform will take 1.23ms, or a processing rate of 0.83 million samples per second. The memory sometimes is the most expensive device in a circuit. If a 16-bit number representation is used, a 16-point transform will need only 96 bytes of buffer memory, a 1024-point transform will need only 6 Kbytes of buffer memory, so that, this architecture is very economical. 4.3.2 Pipeline FHT architecture Figure 4.25 shows the FHT pipeline architecture, which is similar to the single butterfly unit architecture, but instead has one butterfly unit for each stage of the FHT. Figure 4.26 shows the pipeline hardware structure of a 16-point FHT with two but-terfly units in the first two stages and one in the last two stages. The input buffer, output buffer and the butterfly units (BU's) in the first two stages are the same as that in the single butterfly unit structure in the last subsection, except that the outputs of the butterfly units are not stored back to the input buffer, but instead are passed to the next butterfly unit through intermediate buffer memory. Since the last two stages of Chapter 4. The FHT hardware implementation 46 input 1 output 1 Memory Unit Figure 4.25: Pipeline architecture of a 16-point F H T . intermediate buffer intermediate buffer BU2 output buffer output Figure 4.26: Pipeline hardware structure of a 16-point F H T . Chapter 4. The FHT hardware implementation 47 Figure 4.27: Hardware structure of the butterfly unit in the combined last two stages. the radix-2 FHT do not require multiplication of trigonometric coefficients, and they can process the same four data points continuously, it is more economical to combine these two stages into one with a different butterfly unit BU2 as shown in Figure 4.27. The operation of the pipeline FHT structure is similar to the single butterfly unit structure. Each butterfly unit accepts four data points sequentially from the input buffer or the intermediate buffer, and writes the result to the next intermediate buffer or the output buffer. The input buffer and the intermediate buffer both have two memory arrays, one array being used to receive data, and the other to pass data to the next butterfly unit. A switching of input and output between these two arrays will occur when a new transform starts. Figure 4.28 shows the data sequence entering each butterfly unit, with the trigono-metric coefficients listed along. The sequence in each stage can be generated by a read address counter. Since read and write are in subsequent steps, a write address counter is avoided by using eight address registers in each stage to retain the read addresses for eight memory cycles (eight read cycles before the first write cycle), and the addresses in these address registers will be used as the write addresses. Figure 4.29 shows the timing diagram of the pipeline structure operation. The data Chapter 4. The FHT hardware implementation 48 B U B U B U 2 0 1 2 3 4 7 6 5 8 9 1011 12 15 14 13 cosO sinO cos 6 sinO cos 28 sin 26 cos 36 sin 36 0 1 8 9 2 3 1011 4 5 1213 6 7 1415 1 1 t t t cos 0 sinO t cos 0 sinO cos8 sin6 cosB sin 6 0 4 8 12 1 5 9 13 2 6 1014 3 7 11 15 6 = 27t/N Figure 4.28: Data sequence for each butterfly unit. input A: R R| R| R INPUT R R R R INPUT R R • • • buffer B: INPUT R | R | R R INPUT R R R R INPUT • • BU: B| B|B B|B B B B B B B B B B B B B intermediate A: |W|W W | W R R R R w w w w R R R R • « • buffer 1 B: W W W W R R R R W W W W • • • BU: B B B B B B B B B B B intermediate A: W W W W R R R R w W • • • buffer 2 B: W W W W R R • • » BU2: 1 B 1 B B B B 1 output |W W W W 1 ' •* buffer R: 4 read cycles, W: 4 write cycles, B: 1 butterfly cycle Figure 4.29: Timing diagram of the pipeline FHT processor. Chapter 4. The FHT hardware implementation 49 points are input to or output from the butterfly unit sequentially as in the single butterfly unit structure, and the processing time of a transform can be expressed as 4 T ( ^ ) . If memory with access time of 60ns is used, a 16-point transform will take 0.96/xs, or a processing rate of 16.7 million samples per second, while a 1024-point tranform will take 61.4/zs, or a processing rate also of 16.7 million samples per second, since the processing rate of the pipeline structure depends only on the memory read and write access time. In this hardware structure, some very high speed buffer registers are required for the butterfly units, to allow a fast transition of data in the registers before receiving the next four data points. For example, if 60ns memory is used for the buffer memory, 20ns or less access time buffer register should be used. Otherwise, a wait cycle should be inserted between each four memory cycles. 4.3.3 Superparallel F H T architecture In order to obtain an extremely high speed of FHT processing, a superparallel FHT processor is required. Figure 4.30 shows such a system structure for a 16-point FHT. It is in fact the expanded structure of the pipeline FHT in Figure 4.26, enhanced with four butterfly units for each stage. A butterfly unit BUO as shown in Figure 4.31 is used to replace butterfly unit BU in the part without the need of multiplication. In this superparallel structure, the intermediate buffer memory can be eliminated, since it is possible to have hard-wired connections between stages. However, the super-parallel structure requires a large number of butterfly units for a large N, for example, a long FHT transform of 1024 points, requires 2304 butterfly units. This results in a high hardware cost, and because the interconnections between stages are irregular, the input control and inter-stage connections become very complex. By using the superparallel structure, a transform in each memory cycle is possible to obtain. That is, for any N, the processing time will be 60ns if buffer memory of 60ns Chapter 4. The FHT hardware implementation 50 b' Figure 4.31: Hardware structure of butterfly unit BUO. Chapter 4. The FHT hardware implementation 51 e = 2 7t/N Figure 4.32: Pre-permutation 16-point radix-2 FHT data flow diagram. access time is used. This means, theoretically, a processing rate of 267 million samples per second can be obtained for TV = 16, or 17 billion samples per second for TV = 1024. However, this depends on the data input rate as well as the butterfly processing rate. 4.4 Data unscrambling in hardware implementation The transform results obtained by using the methods shown above are scrambled in bit-reversed order; reordering is then required. Since TV is fixed in the hardware structures, the bit-reversing can be done by copying data in the output memory buffer in a bit-reversed order to the output port or to the next processing stage. The bit-reversed copying sequence can be generated by a bit-reversed counter. Another method is to employ pre-permutation algorithm. Figure 4.32 shows the pre-permutation 16-point radix-2 FHT algorithm. If data are entered into the butterfly unit in the reversed-order as shown in Figure 4.33, using the decimation-in-time butterfly as in Chapter 4. The FHT hardware implementation 52 stage 1 0 2 1 3 4 6 5 7 8 10 9 11 12 1413 15 stage 2 0 2 1 3 8 10 9 11 4 6 5 7 1214 13 15 stage 3 8 1 9 12 5 13 10 3 11 14 7 15 1.1 cos 8 cosG sin 6 sin 9 stage 4 0 8 4 12 2 14 6 10 1 9 5 13 3 15 7 11 ttt cos 6 sinO 0 = 2 7 C / N cos 29 sin 29 cos 39 sin 39 Figure 4.33: Data sequence entering the pre-permutation algorithm. Figure 2.1 in the butterfly unit, and the hardware stages in the pipeline and superparallel structures are going in the reversed direction, the transform result will be in the correct order. 4.5 Comparison of the hardware implementations The DHT hardware implementation has a simple control logic and does not require the bit-reversal processing, because the data flow is very regular. Even though it involved more arithmetic than the FHT, it may be optimum for a relatively small N. The DHT hardware implementation can transform a point in each clock signal, therefore, it can run as fast as the FHT pipeline implementation. A comparison of the FHT architectures in the processing rates, the number of butterfly units and the amounts of memory space required is shown in Table 4.6. Chapter 4. The FHT hardware implementation 53 Architecture Number of butterfly units Speed (M samples/sec.) Memory (16-bit word) Single butterfly unit 1 0.83 R A M : 6 Kbytes ROM: 1 Kbytes Pipelining 9 16.7 R A M : 38 Kbytes ROM: 2 Kbytes Superparallelism 2304 17,067 R A M : 4 Kbytes ROM: 7.2 Kbytes Table 4.6: Comparison of the three FHT hardware implementations for TV = 1024 and memory access time of 60 ns. 4.6 Summary A hardware implementation of the DHT is presented in this chapter. It has a simple con-trol logic and flexible expansibility. However, because the number of arithmetic devices required is proportional to TV2, it is only suitable for small TV. Solution to this is sought by the FHT implementation. Several implementation methods of the FHT are discussed, in which, the most economical method, namely, the single butterfly unit architecture, requires the fewest devices but the method is slow, whereas, the fastest method, namely, the superparallel architecture, requires the most devices. A compromise is the pipeline architecture. Chapter 5 Conclusion 5.1 The work done in this thesis Some properties of the fast Hartley transform (FHT) and its application in convolution and power spectrum were reviewed. The method to obtain a complex FFT through the FHT was derived. In the software implementation, some speed-optimized FHT programs have been de-veloped. The radix-2 and split-radix FHT programs were optimized to reduce the number of trigonometric coefficient calculations, the technique of computing the last few stages separately further improved the program execution speed. Moreover, a table-lookup ver-sion of both the radix-2 and split-radix FHT was implemented. The table-lookup version was proved to have even higher efficiency. These programs were compared with some complex FFT programs, which were implemented in an identical way. Tests of these programs indicated that the FHT is about twice as fast as the complex FFT assigned with redundant zero-valued imaginary part. However, when some special methods are used for the complex FFT to handle the real sequence, or when the real FFT algorithm is used, the FHT does not have the speedup effect. However, the FHT is more convenient in doing the inverse transform, because only one algorithm is required for both forward and inverse transforms. On the other hand, the real sequence on the complex F F T handled with the special methods and on the real FFT, both require a separate algorithm to convert the complex conjugate symmetric transform sequence back to a real 54 Chapter 5. Conclusion 55 sequence. The application of the FHT was shown with an image processing example, where the FHT was used in the convolution and computation of power spectrum. The example used a linear-motion blur operator and white noise to degrade a protrait image, and used a Wiener filter to restore the image. The expected result was obtained as if it was done with the FFT. A new bit-reversal algorithm was proposed based on the C-language bitwise opera-tors. This new algorithm and several other bit-reversal algorithms were implemented and evalulated. The new algorithm is fast and simple in implementation. In the hardware implementation, several new hardware structures were constructed based on some advanced architectures. The hardware implementation is new to the DHT and the FHT, since, from the best of the author's knowledge, no other hardware implementation of the DHT or the FHT has previously been published. The DHT implementation requires a number of arithmetic devices proportional to N 2 , but it has a simple control circuit, and a processing speed of 16.7 million samples per second for a 60ns clock cycle, so that, the DHT hardware implementation is attractive for small TV. The single butterfly unit implementation is the simplest, requires the least hardware, and has a fairly high speed, namely, a data processing speed of 0.83 million samples per second for TV = 1024 and a 60ns buffer memory access time. The pipelining implementation is the most efficient one for either large or small TV, which has a very high processing speed, namely, 16.7 million samples per second for TV = 1024 and a 60ns buffer memory access time, and only requires (log2 TV — 1) butterfly units. The superparallelism implementation is the fastest, theoretically can perform a trans-form in one memory cycle. It requires a hardware count of ^f(log2 TV — 1) butterfly units, Chapter 5. Conclusion 56 which turns out to be a large number if N is large, but it is acceptable when N is small or when a really high processing speed is necessary. 5.2 Future work A regular geometry FHT algorithm is sought to ease the software and hardware imple-mentations. To implement the FHT in software with a regular geometry algorithm, will reduce the coding effort and allow a simplier and more efficient program. To implement the FHT in hardware with a regular geometry, will greatly simplify the control circuit, and reduce the effort in designing the address generation circuit. The pipeline architec-ture can then be improved to accept two or even four input data points in each memory cycle, and hence increase the processing speed to two or four times. The superparal-lel architecture will also benefit since regular geometry allows constant geometry and repeatable interconnection design between stages. In the hardware implementation, as N increases, the hardware size will increase cor-respondingly, such that, for a large N, the chip count will be very high. It is more economical to have a custom-designed VLSI chip in case of mass production. In fact, a VLSI implementation project is currently being scheduled at UBC. In the near future, the FHT processor will very likely become a great competitor of the FFT processors. Two dimensional and multi-dimensional FHT's have not been implemented or eval-uated in this thesis; further development in this area will surely benefit two-dimensional and multi-dimensional digital signal processing. References [1] G.D. Bergland, "A fast Fourier transform algorithm for real-valued series," Comm. ACM, vol. 11, no. 10, pp. 703-710, October 1968. [2] J.-B. Martens, "Discrete Fourier transform algorithms for real valued sequences," IEEE Tran. ASSP, vol. 32, no. 2, pp. 390-396, April 1984. [3] R. Kumaresan and P.K. Gupta, "A prime factor FFT algorithm with real valued arithmetic," Proc. IEEE, vol. 73, no. 7, pp. 1241-1243, July 1985. [4] P. Duhamel, "Implementation of 'split-radix' FFT algorithms for complex, real, and real-symmetric data," IEEE Trans. ASSP, vol. 34, no. 2, pp. 285-295, April 1986. [5] H.V. Sorensen, S.L. Jones, M.T. Heideman, and S. Burrus, "Real-valued fast Fourier transform algorithms," IEEE Trans. ASSP, vol. 35, no. 6, pp. 849-863, June 1987. [6] R.N. Bracewell, The Hartley Transform. NY: Oxford Universtiy Press, March 1986. [7] R.N. Bracewell, "The fast Hartley transform," Proc. IEEE, vol. 72, no. 8, pp. 1010-1018, August 1984. [8] 0. Buneman, "Conversion of FFT's to fast Hartley transform," SIAM J. Sci. Stat. Comput., vol. 7, no. 2, pp. 624-638, April 1986. [9] H.V. Sorensen, D.L. Jones, C.S. Burrus and M.T. Heideman, "On computing the discrete Hartley transform," IEEE Trans. ASSP, vol. 33, no. 4, pp. 1231-1238, October 1985. 57 References 58 [10] H.-J. Meckleburg and P.K. Gupta, "Fast Hartley transform algorithm," Electronics Letters, vol. 21, pp. 341-343, April 1985. [11] J. Prado, "Comments on 'The fast Hartley transform'," Proc. IEEE, vol. 73, no. 12, pp. 1862-1863, December 1985. [12] G.E.J. Bold, "A comparison of the time involved in computing fast Hartley and fast Fourier transforms," Proc. IEEE, vol. 73, no. 12, pp. 1863-1864, December 1985. [13] R.N. Bracewell, 0. Buneman, H. Hao, and J. Villasenor, "Fast two-dimensional Hartley transform," Proc. IEEE, vol. 74, no. 9, pp. 1282-1283, September 1986. [14] H. Hao and R.N. Bracewell, "A three-dimensional DFT algorithm using the fast Hartley transform," Proc. IEEE, vol. 75, no. 2, pp. 264-266, February 1987. [15] M.A. O'Neill, "Faster than fast Fourier," BYTE, pp. 293-300, April 1988. [16] R.N. Bracewell and J. Villasenor, "Hartley vs. Fourier," BYTE, vol. 13, no. 12, pp. 24-30, November 1988. [17] T. Le-Ngoc and M.T. Vo, "Implementation and performance of the fast Hartley transform," IEEE MICRO, pp. 20-27, October 1989. [18] R.V.L. Hartley, "A more symmetrical Fourier analysis applied to transmission prob-lems," Proc. IRE, vol. 30, pp. 144-150, March 1942. [19] R.N. Bracewell, "The discrete Hartley transform," J. Opt. Soc. Am., vol. 73, no. 12, pp. 1832-1835, December 1983. [20] S.-C. Pei and J.-L. Wu, "Split-radix fast Hartley transform," Electronics Letters, vol. 22, no. 1, pp. 26-27, January 2, 1986. References 59 [21] R.N. Bracewell, "Alternative to split-radix Hartley transform," Electronics Letters, vol. 23, no. 21, pp. 1148-1149, October 1987. [22] H.S. Hou, "The fast Hartley transform algorithm," IEEE Trans. Computers, vol. 36, no. 2, pp. 147-156, February 1987. [23] C.-Y. Hsu and J.-L. Wu, "Fast computation of discrete Hartley transform via Walsh-Hadamard transform," Electronic Letters, vol. 23, no. 9, pp. 466-468, Aprl 1987. [24] C.-Y. Hsu and J.-L. Wu, "The Walsh-Hadamard/discrete Hartley transform," Int. J. Electronics, vol. 62, no. 5, pp. 746-755, 1987. [25] C P . Kwong and K.P. Shiu, "Structured fast Hartley transform algorithms," IEEE Trans. ASSP, vol. 34, no. 4, pp. 1000-1002, August 1986. [26] H. Hao, "On fast Hartley transform algorithms," Proc. IEEE, vol. 75, no. 7, pp. 961-962, July 1987. [27] P. Duhamel and H. Hollmann, "Split radix FFT algorithm," Electronics Letters, vol. 20, no. 1, pp. 14-16, January 5, 1984. [28] C.S. Burrus and T.W. Parks. DFT/FFT and Convolution Algorithms. NY: Wiley, 1985. [29] H.V. Sorensen, M.T. Heideman and C.S. Burrus, "On computing the split-radix FFT," IEEE Trans. ASSP, vol. 34, no. 1, pp. 152-156, February 1986. [30] P. Duhamel and H. Hollmann, "Existence of a 2" FFT algorithm with a number of multiplications lower than 2 n + 1," Electronics Letteres, vol. 20, no. 17, pp. 690-692, August 1984. References 60 [31] E.O. Brigham, The Fast Fourier Transform. Englewood Cliff, NJ: Prentice-Hall, 1974. [32] L.R. Rabiner and B. Gold, Theory and Application of Digital Signal Processing. Englewood Cliffs, NJ: Prentice-Hall, 1975. [33] R. Kumaresan and P.K. Gupta, "Vector-radix algorithm for a 2-D discrete Hartley transform," Proc. IEEE, vol. 74, no. 5, pp. 755-757, May 1986. [34] J.J. Rodriguez, "An improved FFT digit-reversal algorithm," IEEE Trans. ASSP, vol. 37, no. 8, pp. 1298-1300, August 1989. [35] D.M.W. Evans, "An improved digit-reversal permutation algorithm for the fast Fourier and Hartley transforms," IEEE Trans. ASSP, vol. 35, no. 8, pp. 1120-1125, August 1987. [36] D.M.W. Evans, "A second improved digit-reversal permutation algorithm for fast transforms," IEEE Trans. ASSP, vol. 37, no. 8, pp. 1288-1291, August 1989. [37] Texas Instruments, "New 'C30 FFT benchmark," Details on Signal Processing, Issue 20, pp. 5, fall 1989. [38] Z.M. Ali, "A high-speed FFT processor," IEEE Trans. Communications, vol. 26, no. 5, pp. 690-696, May 1978. [39] E.E. Swartzlander, Jr. and K.W. Young, "A radix-4 delay communtator for fast Fourier transform processor," IEEE J. Solid-State Circuits, vol. 19, no. 5, pp. 702-709, October 1984. [40] K. Yamashita, et al, "A wafer-scale 170000-gate FFT processor with built-in test circuits," IEEE J. Solid-State Circuits, vol. 23, no. 2, pp. 336-342, April 1988. Appendix A The radix -2 decimation-in-frequency F H T programs listing 61 Appendix A. The radix-2 decimation-in-frequency FHT programs listing A . l R a d i x - 2 F H T i m p l e m e n t e d w i t h t h e " t h r e e - l o o p m e t h o d " #include <math.h> #define BUTTERFLY(a, b) {double temp = a; a = temp + b; b = temp - b;> /****************************************************** * Radix-2 decimation-in-frequency FHT * * (three-loop method) * * filename: r2.1.c * * * * Written by Y.K. Fu, Dec. 1989 * * * * Input: * * x -- data arrray with s i z e equals power of two. * * power -- power of two. * * siz e -- number of input points. * * Output: * * scrambled discrete Hartley transform stored i n x. * ******************************************************/ void FHT(double x [ ] , unsigned power, unsigned size) i unsigned n, block, stage, t l , t2, t3, t4, h a l f , q r t r , step; double Red, c, s, A, x1, x2, x3, x4; s t a t i c double TuoPi=6.28318530717958647692; /* check data s i z e < 4 */ i f (power < 2) { i f (power == 1) { BUTTERFLY(x[0], x[1]); return; > else return; > /* f i r s t (power - 2) stages */ half = s i z e ; for (stage = 2; stage < power; stage++) { /* jump by stage step = half; half » = 1; qrtr = half » 1; Rad = TwoPi / step; for (n = 1; n < q r t r ; n++) < /* jump by point A = Rad * n; c = cos(A); s = sin(A); for (block = n; block < s i z e ; block += step) { /* jump by block i f (n == 1) < t l = block - 1; t2 = t1 + q r t r ; t3 = t2 + q r t r ; t4 = t3 + q r t r ; BUTTERFLY(x[t1], x [ t 3 ] ) ; BUTTERFLY(x[t2], x [ t 4 ] ) ; > t3 = block + half; t2 = t3 - n - n; t4 = t2 • half; x1 = x[block] - x [ t 3 ] ; x2 = x[t2] - x[ t 4 ) ; x [block] += x [ t 3 ] ; x[t2) += x [ t 4 ] ; x[t3] = x1 * c + x2 * s; x[t4] = x1 * s - x2 * c; > Appendix A. The radix-2 decimation-in-frequency FHT programs listing 63 / * Last two stages * / for (block a 0; block < size; block += 4) { t1 = block • 1; t2 = block + 2; t3 = block + 3; x1 = x [block] • xCt2]; x2 = x[t1] + xlt3]; x3 = x[block] - x[t2]; x4 = x[t1] - x[t3); x[block] = x1 + x2; x[t1] = x1 - x2; x[t2] = x3 • x4; x[t3] = x3 - x4; > > Appendix A. The radix-2 decimation-in-frequency FHT programs listing A.2 R a d i x-2 F H T i m p l e m e n t e d w i t h t h e " t a b l e - l o o k u p m e t h o d #include <math.h> #include <stdio.h> #define BUTTERFLY(a, b) {double temp = a; a = temp + b; b = temp - b;> /****************************************************** * Radix-2 decimation-in-frequency FHT * * (table-lookup method) * * filename: r2.2.c * • * * Written by Y.K. Fu, Dec. 1989 * * * * Input: * * x data arrray with s i z e equals power of two. * * power -- power of two. * * size -- number of input points. * * Output: * * scrambled discrete Hartley transform stored i n x. * ******************************************************j void FHT(double x [ ] , unsigned power, unsigned size) < unsigned n, block, stage, m, t , t 1 , t2, t3, t4, half, q r t r , step; double Rad, *c, *s, A, x l , x2, x3, x4; s t a t i c double TwoPi=6.28318530717958647692; /* check data s i z e < 4 */ i f (power < 2) { i f (power == 1) { BUTTERFLY(x[0], x[1]); return; > else return; /* Create lookup tables */ i f ( s i z e > 4) { Rad = TwoPi / s i z e ; q r t r = si z e » 2; n = sizeof(double); c = (double * ) c a l l o c ( q r t r , n); s = (double * ) c a l l o c ( q r t r , n); i f (s == NULL) { putsC'** run out of memory i n creating lookup tables * * " ) ; e x i t ( 1 ) ; > for (m = 1; m < q r t r ; m++) { A = Rad * m; c[m] = cos(A); s[m] = sin(A); > /* f i r s t (power - 2) stages */ half = s i z e ; for (stage = 2; stage < power; stage++) < /* jump by stage step = half; half » = 1; qrtr « half » 1; t = stage - 2; Appendix A. The radix-2 decimation-in-frequency FHT programs listing 65 for (n = 1; n < q r t r ; n++) { /* jump by point */ m = n « t ; for (block = n; block < s i z e ; block += step) { /* jump by block */ i f (n == 1) < t1 = block - 1; t2 = t1 + q r t r ; t3 = t2 + q r t r ; t4 = t3 + q r t r ; BUTTERFLY(x[t1], x [ t 3 ] ) ; BUTTERFLY(x[t2], x [ t 4 ] ) ; > t3 = block + h a l f ; t2 = t3 - n - n; t4 = t2 + half; x1 = x[block] - x [ t 3 ] ; x2 = x[t2] - x t t 4 ] ; x [block] += x [ t 3 ] ; x[t2] +="x[t4];" x[t3] = x1 * c[m] + x2 * s[m]; x[t4] = x1 * s[m] - x2 * c[m]; > > > /* Last two stages */ for (block = 0; block < s i z e ; block •= 4) { t1 = block • 1; t2 = block + 2; t3 = block + 3; x1 = x [block] + x [ t 2 ] ; x2 = x t t l ] + x [ t 3 ] ; x3 = x[block] - x [ t 2 ] ; x4 = x[t1] - x [ t 3 ] ; x[block] = x1 + x2; x[t1] = x l - x2; x[t23 = x3 + x4; x[t3] = x3 - x4; > > A p p e n d i x B T h e split-radix decimation-in-frequency F H T programs listing 66 Appendix B. The split-radix decimation-in-frequency FHT programs listing B.l Split-radix FHT implemented with the "three-loop method" #include <math.h> #define BUTTERFLYCa, b) {double temp = a; a = temp + b; b = temp - b;> ^****************************************************** * S p l i t - r a d i x decimation-in-frequency FHT * * (three-loop method) * * filename: s r . l . c * * * * Written by Y.K. Fu, Dec. 1989 * * . * * The looping technique i s adopted from the FORTRAN * * subroutine FFT written by C.S. Burrus, D^ ec. 1984 i n * * "DFT/FFT and Convolution Algorithms", NY: Wiley, * * pp.141-142, 1985 * * * * Input: * * x -- data arrray with s i z e equals power of two. * * power -- power of two. * * siz e -- number of input points. * * Output: * * scrambled discrete Hartley transform stored i n x. * ******************************************************^ void FHT(double x [ ] , unsigned power, unsigned size) { unsigned q r t r , octa, stage, n, block, blocksize, s t a r t , step, stepsize, t11, t12, t13, t14, t21, t22, t23, re-double Rad, c, s, c3, s3. A, A3, xO, x l , x2, x3, x4; s t a t i c double TwoPi=6.28318530717958647692, Root2=1.41421356237309514547; /* Check data s i z e < 8 V i f (power < 3) { i f (power == 2) € x1 = x[0J + x[2]; x2 = x t l ) • x[3]; x3 = x[0] - x[2]; x4 = x t l l - x[3]; x[0] = x1 + x2; xt1] = x l - x2; x[2] = x3 + x4; x[3] = x3 - x4; return; > i f (power == 1) { BuTTERFLY(x[0], x t l ] ) ; return; > else return; /* f i r s t (power - 3) stages */ blocksize = siz e « 1; for (stage = 3; stage < power; stage**) { /* jump by stage */ stepsize = blocksize; blocksize » = 1; qrtr - blocksize » 2; octa = q r t r » 1; Rad = TwoPi / blocksize; Appendix B. The split-radix decimation-in-frequency FHT programs listing 68 for (n = 1; n < octa; n++) ( /* jump by point */ A = Rad * n; A3 = A * 3; c = cos(A); s = sin(A); c3 = cos(A3); s3 = sin(A3); s t a r t = n; step = stepsize; do { /* jump by block */ for (block = s t a r t ; block < s i z e ; block •= step) { i f (n == 1) i t l 1 = block - 1; t12 » t11 + q r t r ; t13 = t12 + q r t r ; t14 = t13 + q r t r ; x3 = x[t11] - x[t13); x4 = x[t!2] - x[t14]; x[t11] += x[t13]; x[t12] += x[t14]; x[t13] ~ x3 + x4; xtt14] = x3 - x4; t21 = t11 + octa; t22 = t21 + q r t r ; t23 = t22 + q r t r ; t24 = t23 + q r t r ; x1 = x[t21J; x2 = x[t22]; x[t21] = x1 • x[t23]; x[t22] = x2 • x[t241; x[t23) = (x1 - xtt23J) * Root2; xlt24] = (x2 - x[t24)) * Root2; t12 = block + q r t r ; t13 = t12 + q r t r ; t14 = t13 + q r t r ; t21 = t12 - n - n; t22 - t21 + q r t r ; t23 = t22 • q r t r ; t24 = t23 + q r t r ; x2 = x[block] - x[t13]; xO = x[t21] - x[t23]; x1 = x2 + xO; x2 -- xO; xO = x[t12] - x[t14); x4 = x[t22] - x[t24]; x3 = x4 - xO; x4 += xO; x [block] += xtt13J; x[t21] += x[t23]; xtt12] •= x[t14]; x[t22] += x[t24]; x[t13] = x1 * c + x3 * s; x[t23] = x1 * s - x3 * c; x[t14] = x2 * c3 + x4 * s3; x[t24] = x2 * s3 - x4 * c3; > s t a r t = (step « 1) - blocksize + n; step « = 2; > while (s t a r t < s i z e ) ; } > Appendix B. The split-radix decimation-in-frequency FHT programs listing /* Last 3 stages */ /* 8-point FHT blocks */ start = 0; step = 16; do { for (block = s t a r t ; block < s i z e ; block += step) { t12 = block + 2; t13 = t12 + 2; t14 = t13 + 2; x3 = x[block] - x[t13]; x4 = x[t12] - x[t14]; x [block] = x [block] + x[t13]; x[t12] = x[t12] • x[t14]; x[t13] = x3 + x4; x[t14] = x3 - x4; t21 = block + 1; t22 = t21 + 2; t23 = t22 + 2; t24 = t23 + 2; x1 = x[t21]; x2 = x[t22]; x[t21] = x1 + x[t23]; x[t22] = x2 + x[t24]; x[t23] = (x1 - x[t23]) * Root2; x[t24] = (x2 - x[t24)) * Root2; BUTTERFLY(x[t13], x[t23]); BUTTERFLY(x[t14], x[t24]); > s t a r t = (step « 1) - 8; step « = 2; > while (s t a r t < s i z e ) ; /* 4-point FHT blocks */ start = 0; step = 8; do < for (block = s t a r t ; block < s i z e ; block += step) t t12 = block + 1; t13 = t12 + 1; t14 = t13 + 1; x1 = x [block] + x[t13); x2 = x[t12] + x[t14]; x3 - x[block] - x[t13]; x4 = x[t12] - x[t14]; x[block] = x l + x2; x[t12] = x1 - x2; x[t13] = x3 + x4; x[t14] = x3 - x4; > start = (step « 1) - 4; step « = 2; > while (s t a r t < s i z e ) ; Appendix B. The split-radix decimation-in-frequency FHT programs listing B.2 Split-radix F HT implemented with the "table-lookup method' #include <math.h> #include <stdio.h> #define BUTTERFLY(a, b) {double temp = a; a = temp + b; b = temp - b;> / *************************»**************************** * S p l i t - r a d i x decimation-in-frequency FHT * * (table-lookup method) * * filename: sr.2.c * * * * Written by Y.K. Fu, Dec. 1989 * * • * The looping technique i s adopted from the FORTRAN * * subroutine FFT written by C.S. Burrus, Dec. 1984 i n * * "DFT/FFT and Convolution Algorithms", NY: Wiley * * pp.141-142, 1985 * * * * Input: * * x -- data arrray with s i z e equals power of two. * * power -- power of two. * * size -- number of input points. * * Output: * * scrambled discrete Hartley transform stored i n x. * ******************************************************/ void FHT(double x [ ] , unsigned power, unsigned size) { unsigned q r t r , octa, stage, n, m, t , block, blocksize, s t a r t , step, stepsize, t i l , t12, t13, t14, t21, t22, t23, t24; double Rad, *c, *s, *c3, *s3. A, A3, xO, x1, x2, x3, x4; s t a t i c double TwoPi=6.28318530717958647692, Root2=1.41421356237309514547; /* Check data s i z e < 8 */ i f (power < 3) { i f (power == 2) C x1 = x[0] + x[2]; x2 = x[1] + x[3]; x3 = x[0] - x[2); x4 = x t D - x[3]; xtO] a x1 + x2; x[1] = x1 - x2; x[2] = x3 + x4; x[3] = x3 - x4; return; > i f (power == 1) ( BUTTERFLY(x[0], x[1]>; return; > else return; /* Create lookup tables */ i f ( size > 8) { Rad = TwoPi / s i z e ; octa = s i z e » 3; n = sizeof(double); c a (double * ) c a l l o c ( o c t a , n); s a (double * ) c a l l o c ( o c t a , n); c3 a (double * ) c a l l o c ( o c t a , n); s3 = (double * ) c a l l o c ( o c t a , n); i f (s3 == NULL) < putsC'** run out of memory i n creating lookup tables * * " ) ; e x i t d ) ; > Appendix B. The split-radix decimation-in-frequency FHT programs listing for (m = 1; m < octa; m++) { A = Rad * m; A3 = A * 3; cCm] = cos(A); s[m] = sin(A); c3[m] = cos(A3); s3[m] = sin(A3); > > /* f i r s t (power - 3) stages */ blocksize = s i z e « 1; for (stage = 3; stage < power; stage++) { /* jump by stage */ stepsize = blocksize; blocksize » = 1; qr t r = blocksize » 2; octa = q r t r » 1; t = stage - 3; for (n = 1; n < octa; n++) { /* jump by point */ sta r t = n; step = stepsize; m = n « t ; do { /* jump by block */ for (block * stBrt; block < s i z e ; block += step) { i f (n==1) I t11 s block - 1; t12 » t11 + q r t r ; t13 = t12 + q r t r ; t H = t13 + q r t r ; x3 = x[t11] - x[t13]; x4 = x[t12] - x [ t H ] ; x[t11) •= x[t13); xtt12] += x t t H ] ; xtt13] = x3 + x4; x [ t H ] = x3 - x4; t21 = t11 + octa; t22 e t21 + q r t r ; t23 = t22 + q r t r ; t24 s t23 + q r t r ; x1 = x[t21]; x2 = x[t223; xtt21] = x l + x[t23]; x(t22] = x2 + xtt241; x[t23] = (x1 - xtt23]) * Root2; x[t24] = (x2 - x[t24]) * Root2; > t12 = block + q r t r ; t13 = t12 + q r t r ; t H = t13 • q r t r ; t21 = t12 - n - n; t22 = t21 + q r t r ; t23 = t22 + q r t r ; t24 = t23 + q r t r ; x2 = x[block] - x[t13); xO = x[t21] - x[t23); x l = x2 + xO; x2 -= xO; xO = x[t12] - x [ t H ) ; x4 = xtt22] - x[t24]; x3 = x4 - xO; x4 •= xO; x [block] += x[t13); x[t21] •= x[t23]; x[t12] += x [ t H ] ; x[t22] += x[t24]; x[t13] = x1 * c[m] + x3 * s[m]; x[t23] = x1 * s[m] - x3 * c[m]; Appendix B. The split-radix decimation-in-frequency FHT programs listing x[t14] = x2 * c3[m] + x4 * s3[m]; x[t24] = x2 * s3tm] - x4 * c3[m]; > start - (step « 1) - btocksize + n; step « = 2; > while (s t a r t < s i z e ) ; > > /* Last 3 stages */ /* 8-point FHT blocks */ start = 0; step = 16; do { for (block « s t a r t ; block < s i z e ; block += step) { t12 = block + 2; t13 = t12 + 2; t14 = t13 + 2; x3 = x[block] - x[t13]; x4 = x[t12] - x[t14]; x [block] = x [block] + x[t13]; x[t12] = x[t12] + x[t14]; x[t13] = x3 + x4; x[t14] = x3 - x4; t21 = block • 1; t22 = t21 + 2; t23 = t22 + 2; t24 = t23 + 2; x1 = x[t21]; x2 = x[t22J; x[t21] = x1 + x[t23]; x[t22] = x2 + x[t24]; x[t23] = ( x l - x[t23]) * Root2; x[t24] = (x2 - x[t24]) * Root2; BUTTERFLY(x[t13], x[t23]); BUTTERFLY(xtt14], x[t24]); > start = (step « 1) - 8; step « = 2; > while (s t a r t < s i z e ) ; /* 4-point FHT blocks */ start = 0; step = 8; do < for (block = s t a r t ; block < s i z e ; block += step) { t12 = block + 1; t13 = t12 + 1; t14 = t13 + 1; x1 = x[block] • x[t13]; x2 = x[t12] + x[t14]; x3 = x [block] - x[t13]; x4 = x[t12] - x[t14]; x[block] = x l + x2; x[t12] = x1 - x2; x[t13] = x3 + x4; x[t14] = x3 - x4; > start = (step « 1) - 4; step « = 2; > while (s t a r t < s i z e ) ; Appendix C The radix-2 decimation-in-frequency FFT programs listing 73 Appendix C. The radix-2 decimation-in-frequency FFT programs Usting C l Radix-2 forward FFT implemented with the "three-loop method #include <math.h> #define BUTTERFLY(a, b) {double temp = a; a = temp + b; b = temp • b;> /****************************************************** * Radix-2 decimation-in-frequency FFT * * (three-loop method) * * filename: r 2 f f t . c * * . ..... * * Written by Y.K. Fu, Dec. 1989 * * Input: * * x, real part and y, imaginary part of a data * * sequence with s i z e equals power of two. * * power -- power of two. * * size -- number of input points. * * Output: * * x, real part and y, imaginary part of the * * scrambled discrete Hartley transform. * ************************************************* void FFT(double x [ ] , double y [ ] , unsigned power, unsigned size) { unsigned stage, block, n, hal f , step, t1, t2; double A, Rad, c, s, x_tmp, y_tmp; s t a t i c double TwoPi = 6.28318530717958647692; half = s i z e ; for (stage = 1; stage < power; stage++) { /* jump by stage */ step = ha l f ; half » = 1; Rad = TwoPi / step; for (n = 1; n < half; n++) { /* jump by point */ A = Rad * n; c = cos(A); s = sin(A); for (block = n; block < s i z e ; block += step) { /* jump by block */ i f (n == 1) { t1 = block - 1; t2 = t1 + half; BUTTERFLY(x[t1], x [ t 2 ] ) ; BUTTERFLY(y[t1], y [ t 2 ] ) ; > t2 = block + half; /* (x1 + j y l ) - (x2 + jy2) = (x1 - x2) + j(y1 - y2) */ x_tmp = x[block] - x [ t 2 ] ; y~tmp = y[block] - y [ t 2 ] ; /* (x1 + jy1) + (x2 + jy2) = (x1 + x2) • j(y1 + y2) */ x[block) += x [ t 2 ] ; y [block] += y [ t 2 ] ; /* (x_tmp + jy_tmp) * (cosA - jsinA) */ /* = (x_tmp * cosA • y_tmp * sinA) + j(y_tmp * cosA - x_tmp * sinA) */ x[t2] = x_tmp * c + y_tmp * s; y[t2] = y~tmp * c - x tmp * s; > Appendix C. The radix-2 decimation-in-frequency FFT programs listing I* last stage V for (block = 0; block < s i z e ; block += 2) t t1 = block + 1; BUTTERFLY(x[block], x[t1]); BUTTERFLY(y[block], y[t1]); Appendix C. The radix-2 decimation-in-frequency FFT programs listing 76 C .2 R a d i x - 2 f o r w a r d F F T i m p l e m e n t e d w i t h the " t a b l e - l o o k u p m e t h o d " #include <stdio.h> #include <math.h> #define BUTTERFLY(a, b> {double temp - a; a = temp + b; b = temp - b;> /****************************************************** * Radix-2 decimation-in-frequency FFT * * (table-lookup method) * * filename: r2 f f t . 2 . c * * * * Written by Y.K. Fu, Dec. 1989 * * - * * Input: * * x, real part and y, imaginary part of a data * * sequence with s i z e equals power of two. * * power -- power of two. * * size -- number of input points. * * Output: * * x, real part and y, imaginary part of the * * scrambled discrete Hartley transform. * ******************************************************/ void FFT(double x [ ] , double y [ ] , unsigned power, unsigned size) { unsigned stage, block, n, m, ha l f , step, t , t l , t2; double A, Rad, *c, *s, x_tmp, y tmp; s t a t i c double TwoPi = 6.28318530717958647692; /* Create lookup tables */ Rad = TwoPi / s i z e ; half = si z e » 1; n = sizeof(double); c = (double * ) c a l l o c ( h a l f , n); s = (double * ) c a l l o c ( h a l f , n); i f (s == NULL) { puts("** run out of memory i n creating lookup tables * * " ) ; e x i t ( 1 ) ; > for (m = 1; m < half; m*+) { A = Rad * m; c[m] = cos(A); slm] = sin(A); > half = s i z e ; for (stage =1; stage < power; stage++) { /* jump by stage */ step = half; half » = 1; t = stage - 1; for (n = 1; n < half; n++) { /* jump by point */ m = n « t ; for (block = n; block < s i z e ; block *- step) { /* jump by block */ i f (n == 1) { t1 = block - 1; t2 = t1 • half; BUTTERFLY(x[t1], x [ t 2 ] ) ; BUTTERFLY(y[t1], y [ t 2 ] ) ; > t2 = block + h a l f ; Appendix C. The radix-2 decimation-in-frequency FFT programs listing 77 /* <x1 + jy1) - (x2 + jy2) = <x1 - x2) + j(y1 - y2) */ x_tmp = x[block] - x [ t 2 ] ; y~tmp = y[block] - y [ t 2 ] ; /* (x1 + jy1) + <x2 + jy2) = (x1 + x2) + j(y1 + y2) */ x [block] •= x [ t 2 ] ; y [block] += y [ t 2 ] ; /* (x_tmp + jy_tmp) * (cosA - jsinA) */ /* = (x_tmp * cosA + y_tmp * sinA) + j(y_tmp * cosA - x_tmp * sinA) */ x[t2] = x_tmp * c[m] + y_tmp * s[m]; y[t2] = y_tmp * c[m] - x tmp * s[m]; > > /* last stage */ for (block »167^b^^Fizey"&l^k-+=_2'>"l: t1 = block + 1; BUTTERFLY(X[block] , X[t1]); BUTTERFLY(y[block), y [ t 1 ] ) ; > > Appendix D The split-radix decimation-in-frequency F F T programs listing 78 Appendix D. The split-radix decimation-in-frequency FFT programs listing D.l Split-radix forward FFT implemented with the "three-loop method #include <math.h> #define BUTTERFLY(a, b) {double temp = a; a = temp + b; b = temp - b;> /****************************************************** * S p l i t - r a d i x decimation-in-frequency FFT * * (three-loop method) * * filename: s r f f t . c * * . . ._* * This program i s modified from the FORTRAN * * subroutine FFT written by C.S. Burrus, Dec. 1984 i n * * "DFT/FFT and Convolution Algorithms", NY: Wiley, * * pp.141-142, 1985 * * * * Rewritten i n C with extensive modification * * by Y.K. Fu, Dec. 1989 * * * * Input: * * x, real part and y, imaginary part of a data * * sequence with s i z e equals power of two. * * power -- power of two. * * siz e -- number of input points. * * Output: * * x, real part and y, imaginary part of the * * scrambled discrete Hartley transform. * *************************************************** / void FFT(double x [ ] , double y [ ] , unsigned power, unsigned size) { unsigned q r t r , stage, n, s t a r t , block, blocksize, step, stepsize, tO, t1, t2, t3; double c, c3, s, s3, A, A3, x1, x2, y1, y2, Re1, Re2, Im1, Im2, Rad; s t a t i c double TwoPi=6.28318530717958647692; /* check data s i z e < 4 */ i f (power < 2) { i f (power == 1) { BUTTERFLY(xtO], x[1]); BUTTERFlY(y[0], y[1]); return; > else return; > /* f i r s t (power - 2) stages */ blocksize - siz e « 1; for (stage = 2; stage < power; stage++) { /*jump by stage */ stepsize = blocksize; blocksize » = 1; qrtr = blocksize » 2; Rad = TwoPi / blocksize; for (n = 1; n < q r t r ; n++) { /* jump by point */ A = n * Rad; A3 = A * 3; c = cos(A); s = sin(A); c3 = cos(A3); s3 = sin(A3); start = n; step = stepsize; Appendix D. The split-radix decimation-in-frequency FFT programs listing 80 do { /* jump by block */ for (block = s t a r t ; block < s i z e ; block += step) { i f (n == 1) < tO « block - 1; t l = tO + q r t r ; t2 = t l + q r t r ; t3 = t2 • q r t r ; x1 n X [ t 0 ] - x [ t 2 ] ; x2 « x t t l ) - x [ t 3 ] ; x[t0] += x [ t 2 ] ; x t t l ] •= x [ t 3 ] ; y l = yttO] - yl t 2 ] ; y2 = y[t1] - y [ t 3 ] ; y[t0] += y [ t 2 ] ; y[t1] += y [ t 3 ] ; x[t2] = x1 + y2; y[t2] = y1 - x2; x[t3] = x1 - y2; y[t3] = y1 + x2; > t1 = block • q r t r ; t2 = t l + q r t r ; t3 = t2 • q r t r ; x1 = x [block] - x [ t 2 ] ; x [block] += x [ t 2 ] ; x2 = x[t1] - x [ t 3 ] ; x[t1] •= x [ t 3 ] ; y1 = y[block] - y [ t 2 ] ; y [block] += y t t 2 ) ; y2 = y[t1) - y [ t 3 ] ; y[t1] •= y [ t 3 ) ; /* ( x l + j y l ) - j(x2 + jy2) = (x1 + y2) + j(y1 - x2) V Re1 = x1 + y2; Im1 = y1 - x2; /* (x1 + j y D + j(x2 + jy2) = (x1 - y2) + j(y1 + x2) */ Re2 = x l - y2; Im2 = y1 • x2; /* (Re1 + jlmD * (cosA - jsinA) */ /* • (Re1 * cosA + Im1 * sinA) + j(Im1 * cosA - Re1 * sinA) */ x[t2] = Rel * c + Im1 * s; y[t2] - Im1 * c - Re1 * s; /* (Re2 • jlm2) * (cos3A - jsin3A) */ /* = (Re2 * cos3A + Im2 * sin3A) + j(Im2 * cos3A - Re2 * sin3A) */ x[t3] = Re2 * c3 + Im2 * s3; y[t3] = Im2 * c3 - Re2 * s3; > /* new block s t a r t address */ s t a r t = (step « 1) - blocksize + n; step « = 2; > while (s t a r t < s i z e ) ; > > Appendix D. The split-radix decimation-in-frequency FFT programs listing 81 /* last 2 stages */ /* 4-point FFT blocks */ start = 0; step = 8; do i for (block = s t a r t ; block < s i z e ; block += step) t t1 = block + 1; t2 = t l + 1; t3 = t2 + 1; x1 = xCblock] - x[t21; x2 = x[t1] - x [ t 3 ] ; x[block] += xEtZJ; x[t1] += x [ t 3 ] ; y1 = y[block) - y [ t 2 ] ; y2 = y[t1] - y [ t 3 ] ; y [block] += y [ t 2 ] ; y[t1] += y [ t 3 ] ; x[t2] = x1 + y2; y[t2] = y1 - x2; x[t3] = x1 - y2; y[t3] = y i • x2; > s t a r t = (step « 1) - 4; step « = 2; > while (s t a r t < s i z e ) ; /* 2-point FFT blocks */ start = 0; step = 4; do { ' ' ~ for (block = s t a r t ; block < s i z e ; block •= step) t t1 = block + 1; BUTTERFLY(x[block], x [ t 1 ] ) ; BUTTERFLY(y[block], y [ t 1 ] ) ; > start = (step « 1) - 2; step « = 2; > while (s t a r t < s i z e ) ; Appendix D. The split-radix decimation-in-frequency FFT programs listing 82 D.2 Split-radix forward F F T implemented with the "table-lookup method" #include <stdio.h> #include <math.h> #define BUTTERFLY(a, b) {double temp = a; a = temp + b; b = temp • b;> ^****************************************************** * S p l i t - r a d i x decimation-in-frequency FFT * * (table-lookup method) * * filename: s r f f t . 2 . c * * - * * This program i s modified from the FORTRAN * * subroutine FFT written by C.S. Burrus, Dec. 1984 i n * * "DFT/FFT and Convolution Algorithms", NY: Wiley, * * pp.141-142, 1985 * • . • * Rewritten i n C with extensive modification * * by Y.K. Fu, Dec. 1989 * * . * * Input: * * x, real part and y, imaginary part of a data * * sequence with si z e equals power of two. * * power -- power of two. * * siz e -- number of input points. * * Output: * * x, real part and y, imaginary part of the * * scrambled discrete Hartley transform. * ******************************************************/ void FFT(double x [ J , double y [ ] , unsigned power, unsigned size) { unsigned q r t r , stage, n, m, s t a r t , block, blocksize, step, stepsize, t , tO, t 1 , t2, t3; double *c, *c3, *s, *s3. A, A3, x1, x2, y1, y2, Re1, Re2, Im1, Im2, Rad; s t a t i c double TwoPi=6.28318530717958647692; /* check data s i z e < 4 */ i f (power < 2) { i f (power == 1) { BUTTERFLY(xCO), x[1]); BUTTERFLY(y[0j, y[1]); return; > else return; /* Create lookup tables V Rad = TwoPi / s i z e ; q r t r = siz e » 2; n = sizeof(double); c = (double * ) c a l l o c ( q r t r , n); s = (double * ) c a l l o c ( q r t r , n); c3 = (double * ) c a l l o c ( q r t r , n); s3 = (double * ) c a l l o c ( q r t r , n); i f (s3 == NULL) { putsC'** run out of memory i n creating lookup tables * * " ) ; e x i t ( 1 ) ; > for (m = 1; m < q r t r ; m++) { A = Rad * m; A3 = A * 3; c[m] = cos(A); s[m] = sin(A); c3[m] = cos(A3); s3[m] = sin(A3); > Appendix D. The split-radix decimation-in-frequency FFT programs listing 83 /* F i r s t (power - 2) stages */ blocksize = s i z e « 1; for (stage = 2; stage < power; stage++) ( /•jump by stage */ stepsize = blocksize; blocksize » = 1; qr t r = blocksize » 2; t = stage - 2; for (n = 1; n < q r t r ; n++) < st a r t = n; step = stepsize; m = n « t ; /* jump by point */ do { /* jump by block */ for (block = s t a r t ; block < s i z e ; block += step) < i f (n == 1) { to = block - 1; t1 = tO + q r t r ; t2 = t1 + q r t r ; t3 = t2 + q r t r ; x1 = x[tO) - x( t 2 ) ; x2 = x t t l ] - x [ t 3 ] ; xttO] += x [ t 2 ] ; x t t l ] += x [ t 3 ] ; y1 = y[tO] - y [ t 2 ] ; y2 = y(t1] - y£t3]; yltO] += y [ t 2 ] ; y[t1] += y [ t 3 ] ; x[t2) = x l + y2; y[t2] = y1 - x2; x[t3] = x1 - y2; y[t3] = y1 + x2; > t1 = block + q r t r ; t2 = t1 + q r t r ; t3 = t2 + q r t r ; x1 = x[block] - x [ t 2 ] ; x [block] += x [ t 2 ] ; x2 • x[t1] - x [ t 3 ] ; x[t1] += x [ t 3 ] ; y1 = y[block] - y [ t 2 ] ; ylblock] += y[ t 2 ] ; y2 = y[t1] - y t t 3 ) ; y[t1] += y [ t 3 ] ; /* (x1 + j y D - j(x2 + jy2) = (x1 + y2) + j(y1 - x2) */ Rel = x1 + y2; Im1 = y1 - x2; /* (x1 + jyD + j(x2 * jy2) = (x1 - y2) + j(y1 + x2) */ Re2 = x1 - y2; Im2 = y1 + x2; /* (Re1 + jlmD * (cosA - jsinA) V /* = (Re1 * cosA + Im1 * sinA) + j(Im1 * cosA - Re1 * sinA) */ x[t2] = Re1 * c[m] • Im1 * s[m]; y[t2] = Im1 * c[m] - Re1 * s[m]; Appendix D. The split-radix decimation-in-frequency FFT programs listing I* (Re2 + jlm2) * (cos3A - jsin3A) */ /* = (Re2 * cos3A «• Im2 * sin3A) + j(Im2 * cos3A - Re2 * sin3A) */ x[t33 = Re2 * c3tm] + Im2 * s3[m]; y[t3] = Im2 * c3[m] - Re2 * s3[m]; > /* new block s t a r t address */ s t a r t = (step « 1) - blocksize + n; step « = 2; > while (s t a r t < s i z e ) ; > > /* last 2 stages */ /* 4-point FFT blocks */ start = 0; step = 8; do { for (block = s t a r t ; block < s i z e ; block += step) C t1 = block + 1; t2 = t1 + 1; t3 = t2 + 1; x1 = x[block] - x [ t 2 ] ; x2 = x[t1] - x [ t 3 ] ; x[block] •= x [ t 2 ] ; x[t1] += x [ t 3 ] ; y1 = y[block] - y [ t 2 ] ; y2 = y[t1] - y [ t 3 ) ; y[block] •= y [ t 2 ] ; ytt1] •= y [ t 3 ) ; x[t2] = x1 + y2; y[t2] = y i - x2; x[t3] = x1 - y2; y[t3] = y1 + x2; > s t a r t = (step « 1) - 4; step « = 2; > while (s t a r t < s i z e ) ; /* 2-point FFT blocks */ start = 0; step = 4; do < for (block = s t a r t ; block < s i z e ; block += step) { t1 = block + 1; BUTTERFLY(x[block], x [ t 1 ] ) ; BUTTERFLY(y[block], y [ t 1 ] ) ; > start = (step « 1) - 2; step « = 2; > while (s t a r t < s i z e ) ; A p p e n d i x E Bit-reversal programs listing 85 endix E. Bit-reversal programs listing The implementation of Evans' bit-reversal algorithm #include <stdio.h> #define BASE 2 #define SUAP(a, b) (double temp = a; a = b; b = temp;} ^**************************************************************** * filename: f _ b i t r . c * * Evans1 bit-reversal permutation (for FHT) * * * * Written by Y.K. Fu, Dec. 1989 * * * * Ref: D.M.W. Evans, "An improved d i g i t - r e v e r s a l permutation * * algorithm for the fast Fourier and Hartley transforms," IEEE * * Trans. ASSP, vol.ASSP-35, no.8, pp.1120-1125, August 1987. * * * * Input: * * x -- bit-reversed data sequence. * * power power of two, where 2"power=size. * * size -- length of data swquence. * * Output: * * x -- reordered data sequence. * A * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * i void reverse(double x [ ] , unsigned power, unsigned size) <. unsigned SeedWidth, OffsetWidth, o f f s e t , GpNum, GpSize, OStemp, RevOffset, RevOffsetO, width, *seed, TableSize; /************************* * define seed table s i z e * *************************/ /* get offset width and seed table width */ SeedWidth = OffsetWidth = power » 1; /* i f power i s odd, seed i s one-bit wider */ if (power & 1) { size » = 1; SeedWidth += 1; > GpSize = size » OffsetWidth; /******************************************** * create seed table according to seed width * /* allocate table space */ seed = (unsigned *)calloc(1 « SeedWidth, sizeof(unsigned)); if (seed == NULL) < puts("** run out of memory i n creating seed table * * " ) ; e x i t d ) ; > /* i n i t i a l i z e the two basic elements */ TableSize = BASE; seedlO] = 0; seed[1] = 1; /* expand seed table */ for (width = 1; width < SeedWidth; width++) { /* append zeros to the exi s t i n g seeds (become f i r s t half then) */ /* and create the second-half seeds */ for (offset = 0; offset < TableSize; offset++) { /* append a zero to the end of the seed */ seed [offset] « = 1; /* the second-half table equals the f i r s t - h a l f table plus one * seed[offset + TableSize] = seedtoffset] + 1; > /* for each extra b i t i n width, the table s i z e w i l l be double */ TableSize « = 1; Appendix E. Bit-reversal programs listing 87 /******************************* * switch bit-reversed elements * *******************************/ for (offset = 1; offset < GpSize; offset++) < /* group zero */ RevOffsetO = seed[offset] « OffsetWidth; /* Theorem 1 */ SWAP(x[offset], x[RevOffsetO]>; /* other groups */ OStemp = o f f s e t ; for (GpNum = 1; GpNum < seed[offset]; GpNum++) { /* Theorem 3 */ RevOffset = RevOffsetO + seed[GpNum]; /* Theorem 2 */ OStemp += GpSize; SWAP(x[OStemp), x[RevOffset]); > > > Appendix E. Bit-reversal programs listing E.2 The implementation of Gold and Rader's bit-reversal algorithm 88 tfdefine SWAPCa, b) {double temp = a; a = b; b = temp;} * GoId-Rader's bi t - r e v e r s a l algorithm * * filename: f_r2.c * * * * Written by: Y.K. Fu, March 1990 * * * * Ref: L.R. Rabiner and B. Gold, Theory and Application of * * D i g i t a l Signal Processing, NJ: Prentice-Hall, pp.367, 1975. * * * * Input: * * x -- bit-reversed data sequence. * * power -- power of two, not used here, just to match the * * argument format. * * N -- length of data sequence. * * Output: * * x -- reordered data sequence. * ***************************************************** void reverseCdouble x [ ] , unsigned power, unsigned N) { int i , j , k, M, N2; j =^0; M = N - 1; N2 = N » 1; for <i = 0; i < M; i++) { i f ( i < j ) SWAP(xti], x t j l ) ; /* generate bit-reversed index j */ k = N2; while (k <= j ) i j -= k; k » = 1; > j |= k; > > Appendix E. Bit-reversal programs listing E.3 The implementation of Rodriguez's bit-reversal algorithm 89 #define SUAP(a, b) {double temp = a; a = b; b = temp;> ^*************************************************************^ * Rodriguez's improved (Gold and Rader's) bit-reversal algorithm * * filename: f r3.c * • T — * * Written by: Y.K. Fu, March 1990 * • * * Ref: J.J. Rodriguez, "An improved FFT di g i t - r e v e r s a l algorithm," * * IEEE Trans. ASSP, vol.ASSP-37, no.8, pp1298-1300, August 1989. * * --* * Input: * * x -- bit-reversed data sequence. * * power -- power of two, where 2"power=N. * * N -- length of data sequence. * * Output: * * x -- reordered data sequence. * *************************************************************** void reverse(double x [ ] , unsigned power, unsigned N) { unsigned i , j , k, N2, l a s t ; last = N - (1 « (power + 1 » 1)); j = 0; N2 = N » 1; for ( i = 1; i < l a s t ; i++) { /* generate bit-reversed index j */ k = N2; while (k <= j ) { j "= k; k » = 1; > j |= k; i f ( i < j ) SWAP(xti], x t j l ) ; > Appendix E. Bit-reversal programs listing E.4 The implementation of Buneman's bit-reversal algorithm 90 #define SUAP(a, b) {double temp = a; a = b; b = temp;} /********************************************************** * Buneman's bit- r e v e r s a l algorithm * * filename: f_r4.c * * .- * * Written by Y.K. Fu, March 1990 * * * * Ref: 0. Buneman, "Conversion of FFT's to fast Hartley * * transforms," SIAM J. S c i . Stat. Comput., vol.7, no.2, * * pp.624-638, A p r i l 1986. * * * * Input: * * x -- bit-reversed data sequence. * * power -- power of two, not used here, just * * to match the argument format. * * N -- length of data sequence. * * Output: * * x -- reordered data sequence. * **********************************************************/ void reverse(double x [ ] , unsigned power, unsigned N) I int i , j , k, M, R[16]; j « 1; for (k = 0; k <= power; k++) { R[k] = j ; j += j ; > M = N - 1; j = 0; for ( i = 1; i < M; i++) t /* generate bit-reversed index j V k = power; do { k -= 1; j -= RCk]; } while <j >= 0); j += R[k+1]; /* i f check for ( i < j ) , i i s only from 1 to N - 3. */ i f ( i > j ) SWAP(x[i], x [ j ] ) ; } > Appendix E. Bit-reversal programs listing E.5 The implementation of the new bit-reversal algorithm 91 #define SUAP(a, b) {double temp = a; a = b; b = temp;} /**************************************************************** * A new bit-reverseI algorithm based on the use of the C * * language bitwise operators. * * filename: f_r5.c * * • * Written by: Y.K. Fu, March 1990 * * * * This program generates the bit-reversed indices with a b i t - * * reversal counter, and swap the elements x[index] and * * x[bit-reversed index] only when (index < bit-reversed index). * * • * Input: * * x -- bit-reversed data sequence. * * power -- power of two, not used here, just to match the * * argument format. * * N -- length of data sequence. * * Output: * * x -- reordered data sequence. * *************************************************************** void reverse(double x [ ] , unsigned power, unsigned N) { int i , j , k, M, N2, S; M = N - 2; j = N2 = N » 1; S = (N2 » 1) " N2; /* unscramble array x */ for ( i = 1; i < M; i++) { i f ( i < j ) SWAP(x[i], x [ j ] ) ; k = S' while'(k <= j ) k = (k » 1) * N2; j •= k; i += 1; i f ( i < j ) SWAP(x[i], x [ j j ) ; j "= N2; > > Appendix E. Bit-reversal programs listing 92 E.6 The implementation of the new algorithm using table-lookup method #define SWAP(a, b) {double temp = a; a = b; b = temp;} /**************************************************************** * A new bit-reversal algorithm based on the use of the C * * language bitwise operators. * * filename: f r5.2.c * * 7 * * Written by: Y.K. Fu, March 1990 * * * * This program generates the bit-reversed indices with a b i t - * * reversal counter, and swap the elements x[index] and * * x[bit-reversed index] only when (index < bit-reversed index). * • * * Input: * * x -- bit-reversed data sequence. * * power -- power of two. * * N -- length of data sequence. * * Output: * * x -- reordered data sequence. * ************************************************************ void reverse(double x [ ] , unsigned power, unsigned N) { int h, i , j , M, N2, S[16]; M = N - 2; SCO] = j = N2 = N » 1; /* generate table for "xor" operation */ for ( i = 1; i < power; i++) S[i] = (S[i-1) » 1) " N2; /* unscramble array x */ for ( i = 1; i < M; i++) { i f ( i < j ) SWAP(xCi), x [ j ] ) ; h = 1; while (S[h] <= j ) ++h; j "= S[h]; i •= 1; i f ( i < j ) SWAP(x[i], x [ j ] ) ; j •= N2; > > Appendix F Image processing example using the F H T 93 Appendix F. Image processing example using the FHT 94 F . l Using the F H T in the processing of image degradation and restoration /****************************************************************** * • * Image processing example using FHT * * * * Filename: image.fht.c * * Written by: Y.K. Fu, Feb. 1990 * * * ******************************************************************* * * * b * * f --->| h |--->(+)---> g g --->| W |---> r * * n * * * * Degradataion Restoration * ******************************************************************* * * * Notations * * * * F = FHT(f) f: o r i g i n a l image • * H = FHT(h) h: blur operator • * B = FHT(b) b: blurred image * * N = FHT(n) n: white noise * * G = F * H + N * * g = IFHT(G) g: blurred and noisy image • • W = Wiener f i l t e r function • • R = G * W * * * r = IFHT(R) r: restored image * • ******************************************************************* * INPUT: binary image f i l e with number of points equals power of * • two. • • * * OUTPUT: degraded binary image f i l e and restored binary image * * f i l e . • * * ******************************************************************/ ^include <math.h> /* Define s i z e of data and work arrays */ #define DataSize 16384 /* number of points */ tfdefine Worksize (DataSize « 1) /* Program constants */ #define H OPERATOR 0.1111111111 #define ZERO 0.0 /* White noise generator constants */ tfdefine LONG 2147483648.0 /* the largest random number */ #define LEVEL 10.0 /* maximum noise gray level */ #define seed 654321 /* random number seed */ /* Macro d e f i n i t i o n s */ #define S0R(y> ((y)*(y)) /* square */ iKdefine ABS(y) (y < ZERO) ? -(y) : (y) /* absolute value */ Ux F. Image processing example using the FHT I* Function declarations */ void FHT(double x [ ] , unsigned power, unsigned s i z e ) ; void input (double x [ ] , double f [ ] , unsigned s i z e ) ; void output(double x [ ] , unsigned s i z e ) ; void white(double n [ ] , double MaxLvl, unsigned s i z e ) ; void convolution(double a t ] , double b [ ] , unsigned s i z e ) ; void rmse(double x t l , double f [ ] , unsigned s i z e ) ; void f i l t e r ( d o u b l e x [ ] , double X2[], double h [ ] , double n[], unsigned s i z e ) ; void tables(unsigned s i z e ) ; mainO t unsigned power, i ; double fCDataSize), /* a copy of the o r i g i n a l image */ x[WorkSize], /* image(work) array */ X2[Worksize), /* a copy of the FHT of f */ h[WorkSize], /* blur operator array */ n[WorkSize]; /* white noise array */ / A * * * * * * * * * * * * * * * * * * * * * * * * * i n i t i a l i z a t i o n *******************************/ /* Compute power of two */ i = Worksize; for (power=0; i>1; power++) i » = 1; /* I n i t i a l i z e arrays */ for (i=DataSize; i<WorkSize; i++) x [ i ] = h[i] = n[i] = ZERO; for (i=0; i O a t a S i z e ; i++) h[i] = ZERO; /* Read i n image and store i t i n a linear array, x */ /* image gray l e v e l : 0 (black) - 255 (white) */ input(x, f, OataSize); /* Create uniform linear motion blur operator, h */ /********************************************* 1 I — > | | | | | | | | | V 9 * * I * *********************************************/ h[0] = h[1] = h[2] = h[3] = h[4] = h[5] = h[6] = h[7] = h[8] = H_OPERATOR; /* Create lookup tables for the FHT */ tables(WorkSize); /*************************** Qegradation *********************************/ /* Get F, FHT of image array f and H, FHT of the blur operator h */ FHT(x, power, WorkSize); FHT(h, power, WorkSize); /* Keep a copy of the transformed image for the NSR calculation */ for (i=0; i<WorkSize; i++) X2[i] = x [ i ] ; /* Get blurred image, B = F * H */ convolution^, h, WorkSize); /* Generate white noise, n */ white(n, LEVEL, DataSize); /* Get N, FHT of n */ FHT(n, power, WorkSize); /* Add noise to the blurred image, G = B + N */ for (i=0; i<WorkSize; i++) x [ i ) •= n [ i ] ; Appendix F. Image processing example using the FHT I* Scale result "/ for (i=0; i<UorkSize; i++) x [ i ] /= WorkSize; /* Show blurred and noisy image g, IFHT of G */ FHT(x, power, WorkSize); for (i=DataSize; i<WorkSize; i++) x [ i ] = ZERO; output(x, DataSize); /***************************** Restoration ******************************/ /* Get G, FHT of g V FHT(x, power, WorkSize); /* Restore image by f i l t e r i n g , R = G * W */ f i l t e r ( x , X2, h, n, WorkSize); /* Get r, I FHT of R */ FHT(x, power, WorkSize); /* Show restored image */ output(x, DataSize); /* Calculate root-mean-square error i n restoration */ rmse(x, f, DataSize); > /******••••**********•*•••**•*• Subrout i nes ********************************/ /******************************************************** * Wiener f i l t e r : W = reverse(H) / ( sqr(|H|) + NSR ) * * Where, * * reverse(H(i)) = H(Worksize - i ) * * sqr(|H(i)|) = Power(H(i)) * * = (sqr(H(i>) + sqr(H(WorkSize - i ) ) / 2 * * NSR = PowerNoise / PowersignaI * * = Power(n) / Power(X2) * ********************************************************/ void f i l t e r ( d o u b l e x [ ] , double X2[], double h[], double n[], unsigned size) < unsigned i , j , hsize; double temp, tp; /* Compute W */ hsize = s i z e » 1; h[0] /= (SQR(ht03) + SQR(n(0]> / SQR(X2[0])); for (i=1; i<hsize; i++) { j = si z e - i ; /* temp = sqr(|H(i)|) + NSR */ temp = 0.5 * (SQR(h[i]) + SQR(h[j])) • (SQR(n[i]) • SOR(n[j])) / (SQR(X2ti]) + SQR(X2[j))); tp = h [ i ) ; h t i ) = h[j] / temp; htj] = tp / temp; > /* F * W */ convolution^, h, s i z e ) ; > /* scale the result */ for (i=0; i<size; i++) x t i ) /= s i z e ; Appendix F. Image processing example using the /*********************************************** * Root-mean-square error of the restored image * *******************************•***************/ void rmse(double x [ ] , double f t ] , unsigned size) { unsigned i ; double sum; sum - ZERO; for (i=0; i<size; i++) sum += SQR(x[i] - f [ i ] ) ; p r i n t f ("root mean square error = % f \ n " , sqrt(sum/size)); /************************ * White noise generator * ************************ j void white(double n [ ] , double MaxLvl, unsigned size) unsigned i ; srand(seed); for (i=0; i<size; i++) n t i ) = rand() / LONG * MaxLvl; > /*************************************************** * FHT convolution * * If either a or b i s symmetrical, the convolution * * of a and b can be reduced to a * b. * • • * INPUT: arrays a and b i n frequency domain, * * siz e of arrays. * * OUTPUT: convolution stored i n array a. * ***************************************************/ void convolution(double a [ ] , double b [ ] , unsigned size) i unsigned i , j , hsize; double temp, t 1 , t2; hsize = s i z e » 1; a CO] *= b[0J; for (i=1; i<hsize; i++) < j = s i z e - i ; t1 = b[i ] + bCj]; t2 = bCi] - b [ j ) ; temp = 0.5 * ( a t i l * t1 + a[j] * t 2 ) ; atj] = 0.5 * (a[j] * t1 + aCi] * - t 2 ) ; at i ] = temp; 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0065455/manifest

Comment

Related Items