Arndt - Algorithms for Programmers (523138), страница 9
Текст из файла (страница 9)
In the normalization loop you precompute1.0/n and multiply as divisions are much slower than multiplications. [FXT: fht fft convolution andsplit radix fft convolution in fft/fftcnvl.cc]Auto (or self) convolution is defined ash=hτ:=a~aX(2.6)ax ayx+y≡τ (n)The corresponding procedure should be obvious. [FXT: fht convolution and fht convolution0 infht/fhtcnvl.cc]CHAPTER 2. CONVOLUTIONS43In the definition of the cyclic convolution (2.1) one can distinguish between those summands where thex + y ‘wrapped around’ (i.e.
x + y = n + τ ) and those where simply x + y = τ holds. These are (followingthe notation in [36]) denoted by h(1) and h(0) respectively. Thenh= h(0) + h(1)(2.7)whereh(0)=Xax bτ −xx≤τh(1)=Xax bn+τ −xx>τThere is a simple way to separate h(0) and h(1) as the left and right half of a length-2 n sequence. Thisis just what the acyclic (or linear) convolution does: Acyclic convolution of two (length-n) sequences aand b can be defined as that length-2 n sequence h which is the cyclic convolution of the zero paddedsequences A and B:A:={a0 , a1 , a2 , .
. . , an−1 , 0, 0, . . . , 0}(2.8)Same for B. Thenhτ2Xn−1:=Ax Bτ −xτ = 0, 1, 2, . . . , 2 n − 1(2.9)x=0Xax byX=ax by +0≤x<nx+y≡τ (2n)x,y<2nXax bywhere the right sum is zero because ax = 0 for n ≤ x < 2n. NowXXXax by =ax bτ −x +ax b2n+τ −x =: Rτ + Sτ0≤x<n(2.10)n≤x<2n(2.11)x>τx≤τwhere the rhs. sums are silently understood as restricted to 0 ≤ x < n.For 0 ≤ τ < n the sum Sτ is always zero because b2n+τ −x is zero (n ≤ 2n + τ − x < 2n for 0 ≤ τ − x < n);(0)the sum Rτ is already equal to hτ . For n ≤ τ < 2n the sum Sτ is again zero, this time because it(1)extends over nothing (simultaneous conditions x < n and x > τ ≥ n); Rτ can be identified with hτ 0(0 ≤ τ 0 < n) by setting τ = n + τ 0 .As an illustration consider the convolution of the sequence {1, 1, 1, 1} with itself: its linear self convolutionis {1, 2, 3, 4, 3, 2, 1, 0}, its cyclic self convolution is {4, 4, 4, 4}, i.e.
the right half of the linear convolutionelement-wise added to the left half.By the way, relation 2.3 is also true for the more general z-transform, but there is no (simple) backtransform, so we cannot turnZ −1 [Z [a] Z [b]]a~b =(2.12)(the equivalent of 2.5) into a practical algorithm.A convenient way to illustrate the cyclic convolution of to sequences is the following semi-symbolicaltable:+-|0:0123456789 10 1112 13 14 150123456789 10 1112 13 14 15CHAPTER 2. CONVOLUTIONS1:2:3:1232343454:5:6:7:456756784565677 88 99 109 10 11 1210 11 12 1311 12 13 146 77 88 99 108 9 10 119 10 11 1210 11 12 1311 12 13 1412 13 14 1513 14 15 014 15 0 115 0 1 28:9:10:11:8 9 10 119 10 11 1210 11 12 1311 12 13 1412 13 14 1513 14 15 014 15 0 115 0 1 212:13:14:15:12 13 14 1513 14 15 014 15 0 115 0 1 20123678441234234534560123123423453456456756786 77 88 99 1013 14 1514 15 015 0 10120123123423453456456756786 77 88 99 108 9 10 119 10 11 1210 11 12 1311 12 13 14The entries denote where in the convolution the products of the input elements can be found:+-|0:1:2:3:012011348...3...24...5 <--= h[5] contains a[2]*b[1]9 <--= h[9] contains a[2]*b[b]Acyclic convolution (where there are 32 buckets 0..31) looks like:+-|0:1:2:3:012345670123123423453456456756786 77 88 99 104:5:6:7:456756786 77 88 99 1089 10 1112 13 14 158 9 10 119 10 11 1210 11 12 1311 12 13 14121314151314151614151617151617188 9 10 119 10 11 1210 11 12 1311 12 13 1412131415131415161415161715161718161718191718192018192021192021228:9:10:11:8 9 10 119 10 11 1210 11 12 1311 12 13 1412131415131415161415161715161718161718191718192018192021192021222021222321222324222324252324252612:13:14:15:12131415161718191718192018192021192021222021222321222324222324252324252624252627252627282627282927282930131415161415161715161718.
. . the elements in the lower right triangle do not ‘wrap around’ anymore, they go to extra buckets.Note that bucket 31 does not appear, it is always zero.The equivalent table for a (cyclic) correlation is+-|0:1:2:3:04:5:6:7:45678:9:10:11:123450 15 14 131 0 15 142 1 0 153 2 1 012131415111213146710 911 1012 1113 123456234512340 15 14 131 0 15 142 1 0 153 2 1 08 79 810 911 1067895678456734562345123489 10 118 79 810 911 1012131415111213146789567810 911 1012 1113 120 15 14 131 0 15 142 1 0 153 2 1 012 13 14 1545673456234512348 79 810 911 1067895678121314151112131410 911 1012 1113 12CHAPTER 2. CONVOLUTIONS12:13:14:15:121314151112131410 911 1012 1113 128 79 810 911 10678956784545673456234512340 15 14 131 0 15 142 1 0 153 2 1 09 10 1112 13 14 15while the acyclic counterpart is:+-|0:1:2:3:04:5:6:7:4567123456780 31 30 291 0 31 302 1 0 313 2 1 02829303127282930262728292526272824252627232425262223242521222324202122231920212218192021171819202829303127282930262728292526272824252627232425262223242521222324282930312728293026272829252627283456234512340 31 30 291 0 31 302 1 0 313 2 1 08:9:10:11:8 79 810 911 1067895678456712:13:14:15:121314151112131410 911 1012 1113 123456234512340 31 30 291 0 31 302 1 0 313 2 1 08 79 810 911 106789567845673456234512340 31 30 291 0 31 302 1 0 313 2 1 0Note that bucket 16 does not appear, it is always zero.Two-dimensional convolution (here 4x4) looks like+-|0:1:2:3:01234567012312302301301245675674674574564:5:6:7:45675674674574568:9:10:11:8 9 10 119 10 11 810 11 8 911 8 9 1012:13:14:15:121314151314151214151213151213148 9 10 119 10 11 810 11 8 911 8 9 1089 10 118 9 10 119 10 11 810 11 8 911 8 9 1012 13 14 15121314151314151214151213151213141213141513141512141512131512131401231230230130124567567467457456121314151314151214151213151213140123123023013012012312302301301245675674674574568 9 10 119 10 11 810 11 8 911 8 9 1089 10 1112 13 14 15while 4x4 correlation would be+-|0:1:2:3:01234567032110322103321047655476654776544:5:6:7:12151413131215141413121515141312032110322103321047655476654776548 9 10 1111 8 9 1010 11 8 99 10 11 812151413131215141413121515141312032110322103321047655476654776548 9 10 1111 8 9 1010 11 8 99 10 11 81215141313121514141312151514131203211032210332108:9:10:11:12:13:14:15:47655476654776548 9 10 1111 8 9 1010 11 8 99 10 11 8121514131312151414131215151413128 9 10 1111 8 9 1010 11 8 99 10 11 8CHAPTER 2.
CONVOLUTIONS2.246Mass storage convolution using the MFAThe matrix Fourier algorithm is also an ideal candidate for mass storage FFTs, i.e. FFTs for data setsthat do not fit into physical RAM1 .In convolution computations it is straight forward to save the transpositions by using the MFA followedby the TMFA. (The data is assumed to be in memory as row0 , row1 , . . . , rowR−1 , i.e. the way array datais stored in memory in the C language, as opposed to the Fortran language.) For the sake of simplicityauto convolution is considered here:Idea 2.1 (matrixfft convolution algorithm) The matrix FFT convolution algorithm:1. Apply a (length R) FFT on each column.(memory access with C-skips)2.
Multiply each matrix element (index r, c) by exp(±2 π i r c/n).3. Apply a (length C) FFT on each row.(memory access without skips)4. Complex square row (element-wise).5. Apply a (length C) FFT on each row (of the transposed matrix).(memory access is without skips)6. Multiply each matrix element (index r, c) by exp(∓2 π i r c/n).7. Apply a (length R) FFT on each column (of the transposed matrix).(memory access with C-skips)Note that steps 3, 4 and 5 constitute a length-C convolution.[FXT: matrix fft convolution in matrixfft/matrixfftcnvl.cc][FXT: matrix fft convolution0 in matrixfft/matrixfftcnvl.cc][FXT: matrix fft auto convolution in matrixfft/matrixfftcnvla.cc][FXT: matrix fft auto convolution0 in matrixfft/matrixfftcnvla.cc]A simple consideration lets one use the above algorithm for mass storage convolutions, i.e.
convolutionsof data sets that do not fit into the RAM workspace. An important consideration is theMinimization of the number of disk seeksThe number of disk seeks has to be kept minimal because these are slow operations which, if occur toooften, degrade performance unacceptably.√The crucial modification of the use of the MFA is not to choose R and C as close as possible to n asusually done. Instead one chooses R minimal, i.e. the row length C corresponds to the biggest data setthat fits into the RAM memory2 . We now analyze how the number of seeks depends on the choice of Rand C: in what follows it is assumed that the data lies in memory as row0 , row1 , .
. . , rowR−1 , i.e. theway array data is stored in the C language, as opposed to the Fortran language convention. Further letα ≥ 2 be the number of times the data set exceeds the RAM size.1 The naive idea to simply try such an FFT with the virtual memory mechanism will of course, due to the non-localityof FFTs, end in eternal hard disk activity2 more precisely: the amount of RAM where no swapping will occur, some programs plus the operating system have tobe there, too.CHAPTER 2. CONVOLUTIONS47In step 1 and 3 of algorithm 2.5 one reads from disk (row by row, involving R seeks) the number ofcolumns that just fit into RAM, does the (many, short) column-FFTs3 , writes back (again R seeks) andproceeds to the next block; this happens for α of these blocks, giving a total of 4 α R seeks for steps 1and 3.In step 2 one has to read (α times) blocks of one or more rows, which lie in contiguous portions of thedisk, perform the FFT on the rows and write back to disk, leading to a total of 2 α seeks.Thereby one has a number of 2 α + 4 α R seeks during the whole computation, which is minimized by thechoice of maximal C.
This means that one chooses a shape of the matrix so that the rows are as big aspossible subject to the constraint that they have to fit into main memory, which in turn means there areR = α rows, leading to an optimal seek count of K = 2 α + 4 α2 .If one seek takes 10 milliseconds then one has for α = 16 (probably quite a big FFT) a total of K · 10 =1056 · 10 milliseconds or approximately 10 seconds. With a RAM workspace of 64 Megabytes4 the CPUtime alone might be in the order of several minutes.
The overhead for the (linear) read and write wouldbe (throughput of 10MB/sec assumed) 6 · 1024M B/(10M B/sec) ≈ 600sec or approximately 10 minutes.With a multi-threading OS one may want to produce a ‘double buffer’ variant: choose the row length sothat it fits twice into the RAM workspace; then let always one (CPU-intensive) thread do the FFTs inone of the scratch spaces and another (hard disk intensive) thread write back the data from the otherscratch-space and read the next data to be processed.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
