Arndt - Algorithms for Programmers (523138), страница 10
Текст из файла (страница 10)
With not too small main memory (and not tooslow hard disk) and some fine tuning this should allow to keep the CPU busy during much of the harddisk operations.9Using a mass storage convolution as described the calculation of the number 99 ≈ 0.4281247·10369,693,100could be done on a 32 bit machine in 1999.
The computation used two files of size 2GigaBytes each andtook less than eight hours on a system with a AMD K6/2 CPU at 366MHz with 66MHz memory.Cf. [hfloat: examples/run1-pow999.txt]2.3Weighted Fourier transformsLet us define a new kind of transform by slightly modifying the definition of the FT (cf. formula 1.1):cck=Wv [a]n−1X:=(2.13)vx ax z x kvx 6= 0∀xx=0where z := e± 2 π i/n . The sequence c shall be called weighted (discrete) transform of the sequence a withthe weight (sequence) v.
Note the vx that entered: the weighted transform with vx = √1n ∀x is just theusual Fourier transform. The inverse transform isa =ax3 real-complex4 allowing=Wv−1 [c]1n vxn−1X(2.14)ck z −x kk=0FFTs in step 1 and complex-real FFTs in step 3.for 8 million 8 byte floats, so the total FFT size is S = 16 · 64 = 1024 MB or 32 million floatsCHAPTER 2. CONVOLUTIONS48This can be easily seen:Wv−1 [Wv [a]]yn−1 n−11 XXvx ax z x k z −y kn vyx=0=k=01n=n−1X n−1Xvxk=0 x=01ax z x k z −y kvyn−1X11ax δx,y nvxn x=0 vyay==£ −1 ¤(cf. section 1.1). That Wv Wv [a] is also identity is apparent from the definitions.Given an implemented FFT it is trivial to set up a weighted Fourier transform:Code 2.3 (weighted transform) Pseudo code for the discrete weighted Fourier transformprocedure weighted_ft(a[], v[], n, is){for x:=0 to n-1{a[x] := a[x] * v[x]}fft(a[],n,is)}Inverse weighted transform is also easy:Code 2.4 (inverse weighted transform) Pseudo code for the inverse discrete weighted Fourier transformprocedure inverse_weighted_ft(a[], v[], n, is){fft(a[],n,is)for x:=0 to n-1{a[x] := a[x] / v[x]}}is must be negative wrt.
the forward transform.[FXT: weighted fft in weighted/weightedfft.cc][FXT: weighted inverse fft in weighted/weightedfft.cc]Introducing a weighted (cyclic) convolution hv byhv= a ~{v} b=Wv−1(2.15)[Wv [a] Wv [b]](cf. formula 2.5)Then for the special case vx = V x one hashv=h(0) + V n h(1)(2.16)(h(0) and h(1) were defined by formula 2.7). It is not hard to see why: Up to the final division by theweight sequence, the weighted convolution is just the cyclic convolution of the two weighted sequences,which is for the element with index τ equal toXXXax bn+τ −x V n+τ(2.17)(ax V x ) (by V y ) =ax bτ −x V τ +x+y≡τ (mod n)x≤τx>τCHAPTER 2. CONVOLUTIONS49Final division of this element (by V τ ) gives h(0) + V n h(1) as stated.√The cases when V n is some root of unity are particularly interesting: For V n = ±i = ± −1 one getsthe so called right-angle convolution:hv= h(0) ∓ i h(1)(2.18)This gives a nice possibility to directly use complex FFTs for the computation of a linear (acyclic)convolution of two real sequences: for length-n sequences the elements of the linear convolution withindices 0, 1, .
. . , n − 1 are then found in the real part of the result, the elements n, n + 1, . . . , 2 n − 1 are theimaginary part. Choosing V n = −1 leads to the negacyclic convolution (or skew circular convolution):hv=h(0) − h(1)(2.19)Cyclic, negacyclic and right-angle convolution can be understood as a polynomial product modulo z n − 1,z n + 1 and z n ± i, respectively (cf.
[2]).[FXT: weighted complex auto convolution in weighted/weightedconv.cc][FXT: negacyclic complex auto convolution in weighted/weightedconv.cc][FXT: right angle complex auto convolution in weighted/weightedconv.cc]The semi-symbolic table (cf. table 2.1) for the negacyclic convolution is+-|0:1:2:3:4:5:6:7:8:9:10:11:12:13:14:15:012345678910111213141501234567891011121314151234567891011121314150-2345678910111213141501-3456789101112131415012-4567891011121314150123-5678910111213141501234-6789101112131415012345-7891011121314150123456-8910111213141501234567-9101112131415012345678-1011121314150123456789-1112131415012345678910-1213141501234567891011-1314150123456789101112-1415012345678910111213-1501234567891011121314-Here the products that enter with negative sign are indicated with a postfix minus at the correspondingentry.√With right-angle convolution the minuses have to be replaced by i = −1 which means the wrap-around(i.e.
h(1) ) elements go to the imaginary part. With real input one thereby effectively separates h(0) andh(1) .Note that once one has routines for both cyclic and negacyclic convolution the parts h(0) and h(1) can becomputed as sum and difference, respectively. Thereby all expressions of the form α h(0) + β h(1) can betrivially computed.2.4Half cyclic convolution for half the price?The computation of h(0) from formula 2.7 (without computing h(1) ) is called half cyclic convolution.Apparently, one asks for less information than one gets from the acyclic convolution.
One might hope tofind an algorithm that computes h(0) and uses only half the memory compared to the linear convolutionor that needs half the work, possibly both. It may be a surprise that no such algorithm seems to beknown currently5 .5 Ifyou know one, tell me about it!CHAPTER 2. CONVOLUTIONS50Here is a clumsy attempt to find h(0) alone: Use the weighted transform with the weight sequencevx = V x where V n is very small. Then h(1) will in the result be multiplied with a small number andwe hope to make it almost disappear. Indeed, using V n = 1000 for the cyclic self convolution of thesequence {1, 1, 1, 1} (where for the linear self convolution h(0) = {1, 2, 3, 4} and h(1) = {3, 2, 1, 0}) onegets {1.003, 2.002, 3.001, 4.000}. At least for integer sequences one could choose V n (more than two times)bigger than biggest possible value in h(1) and use rounding to nearest integer to isolate h(0) .
Alas, evenfor modest sized arrays numerical overflow and underflow gives spurious results. Careful analysis showsthat this idea leads to an algorithm far worse than simply using linear convolution.2.5Convolution using the MFAWith the weighted convolutions in mind we reformulate the matrix (self-) convolution algorithm (idea 2.1):1.
Apply a FFT on each column.2. On each row apply the weighted convolution with V C = e2 π i r/R = 1r/R where R is the totalnumber of rows, r = 0..R − 1 the index of the row, C the length of each row (or, equivalently thetotal number columns)3. Apply a FFT on each column (of the transposed matrix).First consider2.5.1The case R = 2The cyclic auto convolution of the sequence x can be obtained by two half length convolutions (one cyclic,one negacyclic) of the sequences6 s := x(0/2) + x(1/2) and d := x(0/2) − x(1/2) using the formulax~x =1{s ~ s + d ~− d,2s ~ s − d ~− d}(2.20)The equivalent formula for the cyclic convolution of two sequences x and y isx~y=1{sx ~ sy + dx ~− dy ,2sx ~ sy − dx ~− dy }(2.21)wheresxdxsy:= x(0/2) + x(1/2):= x(0/2) − x(1/2):= y (0/2) + y (1/2)dy:= y (0/2) − y (1/2)For the acyclic (or linear) convolution of sequences one can use the cyclic convolution of the zero paddedsequences zx := {x0 , x1 , .
. . , nn−1 , 0, 0, . . . , 0} (i.e. x with n zeros appended). Using formula 2.20 one getsfor the two sequences x and y (with sx = dx = x, sy = dy = y):x ~ac y=zx ~ zy=1{x ~ y + x ~− y,2x ~ y − x ~− y}(2.22)And for the acyclic auto convolution:x ~ac x6 s,=z~zd lower half plus/minus higher half of x=1{x ~ x + x ~− x,2x ~ x − x ~− x}(2.23)CHAPTER 2. CONVOLUTIONS2.5.251The case R = 3Let ω =12(1 +√3) and defineA :=B :=C :=x(0/3) + x(1/3) + x(2/3)x(0/3) + ω x(1/3) + ω 2 x(2/3)x(0/3) + ω 2 x(1/3) + ω x(2/3)Then, if h := x ~ac x, there isx(0/3)=A ~ A + B ~{ω} B + C ~{ω2 } C(1/3)x=A ~ A + ω (B ~{ω} B) + ω (C ~{ω2 } C)x(2/3)=A ~ A + ω (B ~{ω} B) + ω 2 (C ~{ω2 } C)(2.24)2For real valued data C is the complex conjugate (cc.) of B and (with ω 2 = cc.ω) B ~{ω} B is the cc.
ofC ~{ω2 } C and therefore every B ~{} B-term is the cc. of the C ~{} C-term in the same line. Is there anice and general scheme for real valued convolutions based on the MFA? Read on for the positive answer.2.6Convolution of real valued data using the MFAFor row 0 (which is real after the column FFTs) one needs to compute the (usual) cyclic convolution; forrow R/2 (also real after the column FFTs) a negacyclic convolution is needed7 , the code for that task isgiven on page 68.All other weighted convolutions involve complex computations, but it is easy to see how to reduce thework by 50 percent: As the result must be real the data in row number R − r must, because of thesymmetries of the real and imaginary part of the (inverse) Fourier transform of real data, be the complexconjugate of the data in row r.
Therefore one can use real FFTs (R2CFTs) for all column-transforms forstep 1 and half-complex to real FFTs (C2RFTs) for step 3.Let the computational cost of a cyclic (real) convolution be q, thenFor R even one must perform 1 cyclic (row 0), 1 negacyclic (row R/2) and R/2 − 2 complex (weighted)convolutions (rows 1, 2, . . . , R/2 − 1)For R odd one must perform 1 cyclic (row 0) and (R − 1)/2 complex (weighted) convolutions (rows1, 2, .
. . , (R − 1)/2)Now assume, slightly simplifying, that the cyclic and the negacyclic real convolution involve the samenumber of computations and that the cost of a weighted complex convolution is twice as high. Then inboth cases above the total work is exactly half of that for the complex case, which is about what onewould expect from a real world real valued convolution algorithm.For acyclic convolution one may want to use the right angle convolution (and complex FFTs in the columnpasses).2.7Convolution without transposition using the MFA *Section 7.4 explained the connection between revbin-permutation and transposition.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
