Arndt - Algorithms for Programmers (523138), страница 4
Текст из файла (страница 4)
.LL L"¶#∞ µ1a0 X2πkx2πkx√+ak cos+ bk sinLLL 2k=1(1.33)(1.34)(1.35)withck=a0212 (ak12 (ak(k = 0)− ibk ) (k > 0)+ ibk ) (k < 0)(1.36)The discrete Fourier transformThe discrete Fourier transform (DFT) of a sequence f of length n with elements fx is defined byck:=n−11 X√fx eσ 2 π i x k/nn x=0(1.37)fx=n−11 X√ck eσ 2 π i x k/nn(1.38)Back-transform isk=01.41.4.1Radix 2 FFT algorithmsA little bit of notationAlways assume a is a length-n sequence (n a power of two) in what follows:Let a(even) , a(odd) denote the (length-n/2) subsequences of those elements of a that have even or oddindices, respectively.Let a(lef t) denote the subsequence of those elements of a that have indices 0 .
. . n/2 − 1.Similarly, a(right) for indices n/2 . . . n − 1.Let S k a denote the sequence with elements ax e± k 2 π i x/n where n is the length of the sequence a andthe sign is that of the transform. The symbol S shall suggest a shift operator. In the next two sectionsonly S 1/2 will appear. S 0 is the identity operator.CHAPTER 1. THE FOURIER TRANSFORM1.4.214Decimation in time (DIT) FFTThe following observation is the key to the decimation in time (DIT) FFT2 algorithm:For n even the k-th element of the Fourier transform isn−1Xn/2−1ax z x kX=x=0n/2−1a2 x z 2 x k +x=0XX(1.39)x=0n/2−1=a2 x+1 z (2 x+1) kn/2−1a2 x z 2 x k + z kx=0Xa2 x+1 z 2 x k(1.40)x=0where z = e±i 2 π/n and k ∈ {0, 1, .
. . , n − 1}.The last identity tells us how to compute the k-th element of the length-n Fourier transform from thelength-n/2 Fourier transforms of the even and odd indexed subsequences.To actually rewrite the length-n FT in terms of length-n/2 FTs one has to distinguish the cases 0 ≤k < n/2 and n/2 ≤ k < n, therefore we rewrite k ∈ {0, 1, 2, . . . , n − 1} as k = j + δ n2 where j ∈{0, 1, . .
. , n/2 − 1}, δ ∈ {0, 1}.n−1Xn/2−1ax zx (j+δn2)=x=0Xn/2−1a(even)xz2 x (j+δn2)+zj+δn2x=0Xa(odd)z 2 x (j+δx(1.41)x=0 n/2−1n/2−1XX(even) 2 x jjaz+za(odd)z2 x jxx=n2)x=0n/2−1x=0n/2−1x=0x=0XXa(even)z2 x j − zja(odd)z2 x jxxforδ=0(1.42)forδ=1Noting that z 2 is just the root of unity that appears in a length-n/2 FT one can rewrite the last twoequations as theIdea 1.1 (FFT radix 2 DIT step) Radix 2 decimation in time step for the FFT:hihin/2(lef t)F [a]= F a(even) + S 1/2 F a(odd)hihin/2(right)F [a]= F a(even) − S 1/2 F a(odd)(1.43)(1.44)(Here it is silently assumed that ’+’ or ’−’ between two sequences denotes element-wise addition orsubtraction.)The length-n transform has been replaced by two transforms of length n/2.
If n is a power of 2 thisscheme can be applied recursively until length-one transforms (identity operation) are reached. Therebythe operation count is improved to proportional n · log2 (n): There are log2 (n) splitting steps, the workin each step is proportional to n.Code 1.1 (recursive radix 2 DIT FFT) Pseudo code for a recursive procedure of the (radix 2) DITFFT algorithm, is must be +1 (forward transform) or -1 (backward transform):procedure rec_fft_dit2(a[], n, x[], is)// complex a[0..n-1] input// complex x[0..n-1] result{complex b[0..n/2-1], c[0..n/2-1]// workspacecomplex s[0..n/2-1], t[0..n/2-1]// workspace2 alsocalled Cooley-Tukey FFT.CHAPTER 1. THE FOURIER TRANSFORM}15if n == 1 then // end of recursion{x[0] := a[0]return}nh := n/2for k:=0 to nh-1 // copy to workspace{s[k] := a[2*k]// even indexed elementst[k] := a[2*k+1] // odd indexed elements}// recursion: call two half-length FFTs:rec_fft_dit2(s[],nh,b[],is)rec_fft_dit2(t[],nh,c[],is)fourier_shift(c[],nh,is*1/2)for k:=0 to nh-1 // copy back from workspace{x[k]:= b[k] + c[k];x[k+nh] := b[k] - c[k];}The data length n must be a √power of 2.
The result is in x[]. Note that normalization (i.e. multiplicationof each element of x[] by 1/ n) is not included here.[FXT: recursive fft dit2 in slow/recfft2.cc] The procedure uses the subroutineCode 1.2 (Fourier shift) For each element in c[0..n-1] replace c[k] by c[k] times ev 2 π i k/n .
Used withv = ±1/2 for the Fourier transform.procedure fourier_shift(c[], n, v){for k:=0 to n-1{c[k] := c[k] * exp(v*2.0*PI*I*k/n)}}cf. [FXT: fourier shift in fft/fouriershift.cc]The recursive FFT-procedure involves n log2 (n) function calls, which can be avoided by rewriting it ina non-recursive way. One can even do all operations in-place, no temporary workspace is needed atall. The price is the necessity of an additional data reordering: The procedure revbin_permute(a[],n)rearranges the array a[] in a way that each element ax is swapped with ax̃ , where x̃ is obtained from xby reversing its binary digits.
This is discussed in section 7.1.Code 1.3 (radix 2 DIT FFT, localized) Pseudo code for a non-recursive procedure of the (radix 2)DIT algorithm, is must be -1 or +1:procedure fft_dit2_localized(a[], ldn, is)// complex a[0..2**ldn-1] input, result{n := 2**ldn // length of a[] is a power of 2revbin_permute(a[],n)for ldm:=1 to ldn // log_2(n) iterations{m := 2**ldmmh := m/2for r:=0 to n-m step m // n/m iterations{for j:=0 to mh-1 // m/2 iterations{e := exp(is*2*PI*I*j/m) // log_2(n)*n/m*m/2 = log_2(n)*n/2 computationsCHAPTER 1.
THE FOURIER TRANSFORM}}}}16u := a[r+j]v := a[r+j+mh] * ea[r+j]:= u + va[r+j+mh] := u - v[FXT: fft localized dit2 in fft/fftdit2.cc]This version of a non-recursive FFT procedure already avoids the calling overhead and it works in-place.It works as given, but is a bit wasteful. The (expensive!) computation e := exp(is*2*PI*I*j/m) isdone n/2 · log2 (n) times. To reduce the number of trigonometric computations, one can simply swap thetwo inner loops, leading to the first ‘real world’ FFT procedure presented here:Code 1.4 (radix 2 DIT FFT) Pseudo code for a non-recursive procedure of the (radix 2) DIT algorithm, is must be -1 or +1:procedure fft_dit2(a[], ldn, is)// complex a[0..2**ldn-1] input, result{n := 2**ldnrevbin_permute(a[],n)for ldm:=1 to ldn // log_2(n) iterations{m := 2**ldmmh := m/2for j:=0 to mh-1 // m/2 iterations{e := exp(is*2*PI*I*j/m) // 1 + 2 + ...
+ n/8 + n/4 + n/2 = n-1 computationsfor r:=0 to n-m step m{u := a[r+j]v := a[r+j+mh] * ea[r+j]:= u + va[r+j+mh] := u - v}}}}[FXT: fft dit2 in fft/fftdit2.cc]Swapping the two inner loops reduces the number of trigonometric (exp()) computations to n but leadsto a feature that many FFT implementations share: Memory access is highly nonlocal. For each recursionstage (value of ldm) the array is traversed mh times with n/m accesses in strides of mh. As mh is a powerof 2 this can (on computers that use memory cache) have a very negative performance impact for largevalues of n. On a computer where the CPU clock (366MHz, AMD K6/2) is 5.5 times faster than thememory clock (66MHz, EDO-RAM) I found that indeed for small n the localized FFT is slower by afactor of about 0.66, but for large n the same ratio is in favor of the ‘naive’ procedure!It is a good idea to extract the ldm==1 stage of the outermost loop, this avoids complex multiplicationswith the trivial factors 1 + 0 i: Replacefor ldm:=1 to ldn{byfor r:=0 to n-1 step 2{{a[r], a[r+1]} := {a[r]+a[r+1], a[r]-a[r+1]}}for ldm:=2 to ldn{CHAPTER 1.
THE FOURIER TRANSFORM1.4.317Decimation in frequency (DIF) FFTThe simple splitting of the Fourier sum into a left and right half (for n even) leads to the decimation infrequency (DIF) FFT3 :n−1Xn/2−1ax z x kX=x=0ax z x k +x=0ax z x k(1.45)ax+n/2 z (x+n/2) k(1.46)x=n/2n/2−1X=nXn/2−1ax zxk+x=0Xx=0n/2−1X=t)(a(lef+ z k n/2 a(right)) zx kxx(1.47)x=0(where z = e±i 2 π/n and k ∈ {0, 1, . . . , n − 1})Here one has to distinguish the cases k even or odd, therefore we rewrite k ∈ {0, 1, 2, . .
. , n − 1} ask = 2 j + δ where j ∈ {0, 2, . . . , n2 − 1}, δ ∈ {0, 1}.n−1Xn/2−1ax z x (2 j+δ)=x=0Xt)(a(lef+ z (2 j+δ) n/2 a(right)) z x (2 j+δ)xx(1.48)x=0 n/2−1Xt)) z2 x j+ a(right)(a(lefxx=for δ = 0x=0n/2−1Xt)) z2 x j− a(right)z x (a(lefxx(1.49)for δ = 1x=0z (2 j+δ) n/2 = e±π i δ is equal to plus/minus 1 for δ = 0/1 (k even/odd), respectively.The last two equations are, more compactly written, theIdea 1.2 (radix 2 DIF step) Radix 2 decimation in frequency step for the FFT:(even)n/2(odd)n/2F [a]F [a]==hiF a(lef t) + a(right)h³´iF S 1/2 a(lef t) − a(right)(1.50)(1.51)Code 1.5 (recursive radix 2 DIF FFT) Pseudo code for a recursive procedure of the (radix 2) decimation in frequency FFT algorithm, is must be +1 (forward transform) or -1 (backward transform):procedure rec_fft_dif2(a[], n, x[], is)// complex a[0..n-1] input// complex x[0..n-1] result{complex b[0..n/2-1], c[0..n/2-1]// workspacecomplex s[0..n/2-1], t[0..n/2-1]// workspaceif n == 1 then{x[0] := a[0]return}nh := n/2for k:=0 to nh-1{3 alsocalled Sande-Tukey FFT, cf.
[18].CHAPTER 1. THE FOURIER TRANSFORMs[k] := a[k]t[k] := a[k+nh]}18// ’left’ elements// ’right’ elements}for k:=0 to nh-1{{s[k], t[k]} := {(s[k]+t[k]), (s[k]-t[k])}}fourier_shift(t[],nh,is*0.5)rec_fft_dif2(s[],nh,b[],is)rec_fft_dif2(t[],nh,c[],is)j := 0for k:=0 to nh-1{x[j]:= b[k]x[j+1] := c[k]j := j+2}The data length n must be a power of 2.
The result is in x[].[FXT: recursive fft dif2 in slow/recfft2.cc]The non-recursive procedure looks like this:Code 1.6 (radix 2 DIF FFT) Pseudo code for a non-recursive procedure of the (radix 2) DIF algorithm, is must be -1 or +1:procedure fft_dif2(a[],ldn,is)// complex a[0..2**ldn-1] input, result{n := 2**ldnfor ldm:=ldn to 1 step -1{m := 2**ldmmh := m/2for j:=0 to mh-1{e := exp(is*2*PI*I*j/m)for r:=0 to n-1 step m{u := a[r+j]v := a[r+j+mh]a[r+j]:= (u + v)a[r+j+mh] := (u - v) * e}}}revbin_permute(a[],n)}cf. [FXT: fft dif2 in fft/fftdif2.cc]In DIF FFTs the revbin_permute()-procedure is called after the main loop, in the DIT code it wascalled before the main loop.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
