Arndt - Algorithms for Programmers, страница 11

Equipped with thatknowledge an algorithm for convolution using the MFA that uses revbin_permute instead of transposeis almost straight forward:rows=87 Forcolumns=4R odd there is no such row and no negacyclic convolution is needed.CHAPTER 2. CONVOLUTIONS52input data (symbolic format: R00C):0:01231:10001001100210032:20002001200220033:30003001300230034:40004001400240035:50005001500250036:60006001600260037:7000700170027003FULL REVBIN_PERMUTE for transposition:0:040002000600010001:240022002600210022:140012001600110013:34003200360031003500050025001500330003002300130037000700270017003DIT FFTs on revbin_permuted rows (in revbin_permuted sequence), i.e.

unrevbin_permute rows:(apply weight after each FFT)0:010002000300040005000600070001:210022002300240025002600270022:110012001300140015001600170013:31003200330034003500360037003FULL REVBIN_PERMUTE for transposition:0:01231:40004001400240032:20002001200220033:60006001600260034:10001001100210035:50005001500250036:30003001300230037:7000700170027003CONVOLUTIONS on rows (do not care revbin_permuted sequence), no reordering.FULL REVBIN_PERMUTE for transposition:0:010002000300040001:210022002300240022:110012001300140013:31003200330034003500050025001500360006002600160037000700270017003(apply inverse weight before each FFT)DIF FFTs on rows (in revbin_permuted sequence), i.e.

revbin_permute rows:0:040002000600010005000300070001:240022002600210025002300270022:140012001600110015001300170013:34003200360031003500330037003FULL REVBIN_PERMUTE for transposition:0:01231:10001001100210032:20002001200220033:30003001300230034:40004001400240035:50005001500250036:60006001600260037:7000700170027003CHAPTER 2. CONVOLUTIONS53As shown works for sizes that are a power of two, generalizes for sizes a power of some prime.add text2.8TBD:The z-transform (ZT)In this section we will learn a technique to compute the FT by a (linear) convolution.

In fact, thetransform computed is the z-transform, a more general transform that in a special case is identical to theFT.2.8.1Definition of the ZTThe z-transform (ZT) Z [a] = â of a (length n) sequence a with elements ax is defined asâkn−1X:=ax z k x(2.25)x=0The z-transform is a linear transformation, its most important property is the convolution property(formula 2.3): Convolution in original space corresponds to ordinary (element-wise) multiplication inz-space. (See [13] and [17].)Note that the special case z = e±2 π i/n is the discrete Fourier transform.2.8.2Computation of the ZT via convolutionIn the definition of the (discrete) z-transform we rewrite8 the product x k asxkfˆk =n−1X=fx z x k¢1 ¡ 2x + k 2 − (k − x)22= zk2/2x=0n−1X³2fx z x/2´(2.26)2z −(k−x)/2(2.27)x=0This leads to the followingIdea 2.2 (chirp z-transform) Algorithm for the chirp z-transform:1.

Multiply f element-wise with z x2/2.2. Convolve (acyclically) the resulting sequence with the sequence z −xis required here.3. Multiply element-wise with the sequence z k2/22/2, zero padding of the sequences.The above algorithm constitutes a ‘fast’ (∼ n log(n)) algorithm for the ZT because fast convolution ispossible via FFT.8 cf.[2]CHAPTER 2. CONVOLUTIONS2.8.354Arbitrary length FFT by ZTWe first note that the length n of the input sequence a for the fast z-transform is not limited to highlycomposite values (especially n prime is allowed): For values of n where a FFT is not feasible pad thesequence with zeros up to a length L with L >= 2 n and a length L FFT becomes feasible (e.g. L is apower of 2).Second remember that the FT is the special case z = e±2 π i/n of the ZT: With the chirp ZT algorithmone also has an (arbitrary length) FFT algorithmThe transform takes a few times more than an optimal transform (by direct FFT) would take.

The worstcase (if only FFTs for n a power of 2 are available) is n = 2p + 1: One must perform 3 FFTs of length2p+2 ≈ 4 n for the computation of the convolution. So the total work amounts to about 12 times thework a FFT of length n = 2p would cost. It is of course possible to lower this ‘worst case factor’ to 6 byusing highly composite L slightly greater than 2 n.[FXT: fft arblen in chirp/fftarblen.cc]TBD: show shortcuts for n even/odd2.8.4Fractional Fourier transform by ZTThe z-transform with z = eα 2 π i/n and α 6= 1 is called the fractional Fourier transform (FRFT).

Uses ofthe FRFT are e.g. the computation of the DFT for data sets that have only few nonzero elements and thedetection of frequencies that are not integer multiples of the lowest frequency of the DFT. A thoroughdiscussion can be found in [51].[FXT: fft fract in chirp/fftfract.cc]Chapter 3The Hartley transform (HT)3.1Definition of the HTThe Hartley transform (HT) is defined like the Fourier transform with ‘cos + sin’ instead of ‘cos +i · sin’.The (discrete) Hartley transform of a is defined ascck=:=H [a]µ¶12πkx2πkx√ax cos+ sinnnn x=0n−1X(3.1)(3.2)It has the obvious property that real input produces real output,H [a]∈ Rfora∈R(3.3)It also is its own inverse:H [H [a]]= a(3.4)The symmetries of the HT are simply:H [aS ] = H [aS ] = H [aS ]H [aA ] = H [aA ] = −H [aA ](3.5)(3.6)i.e.

symmetry is, like for the Fourier transform, conserved.The n log(n)-implementations of the HT are called fast Hartley transforms (FHT).3.23.2.1Radix 2 FHT algorithmsDecimation in time (DIT) FHTFor a sequence a of length n let X 1/2 a denote the sequence with elements ax cos π x/n + ax sin π x/n(this is the ‘shift operator’ for the Hartley transform).Idea 3.1 (FHT radix 2 DIT step) Radix 2 decimation in time step for the FHT:iihhn/2(lef t)H [a]= H a(even) + X 1/2 H a(odd)iihhn/2(right)H [a]= H a(even) − X 1/2 H a(odd)55(3.7)(3.8)CHAPTER 3. THE HARTLEY TRANSFORM (HT)56Code 3.1 (recursive radix 2 DIT FHT) Pseudo code for a recursive procedure of the (radix 2) DITFHT algorithm:procedure rec_fht_dit2(a[], n, x[])// real a[0..n-1] input// real x[0..n-1] result{real b[0..n/2-1], c[0..n/2-1]// workspacereal s[0..n/2-1], t[0..n/2-1]// workspaceif n == 1 then{x[0] := a[0]return}nh := n/2;for k:=0 to nh-1{s[k] := a[2*k]// even indexed elementst[k] := a[2*k+1] // odd indexed elements}rec_fht_dit2(s[], nh, b[])rec_fht_dit2(t[], nh, c[])hartley_shift(c[], nh, 1/2)for k:=0 to nh-1{x[k]:= b[k] + c[k];x[k+nh] := b[k] - c[k];}}[FXT: recursive fht dit2 in slow/recfht2.cc]The procedure hartley_shift replaces element ck of the input sequence c by ck cos(π k/n) +cn−k sin(π k/n).

Here is the pseudo code:Code 3.2 (Hartley shift) procedure hartley_shift_05(c[], n)// real c[0..n-1] input, result{nh := n/2j := n-1for k:=1 to nh-1{c := cos( PI*k/n )s := sin( PI*k/n ){c[k], c[j]} := {c[k]*c+c[j]*s, c[k]*s-c[j]*c}j := j-1}}[FXT: hartley shift 05 in fht/hartleyshift.cc]Code 3.3 (radix 2 DIT FHT, localized) Pseudo code for a non-recursive procedure of the (radix 2)DIT FHT algorithm:procedure fht_dit2(a[], ldn)// real a[0..n-1] input,result{n := 2**ldn // length of a[] is a power of 2revbin_permute(a[], n)for ldm:=1 to ldn{m := 2**ldmmh := m/2m4 := m/4for r:=0 to n-m step m{for j:=1 to m4-1 // hartley_shift(a+r+mh,mh,1/2)CHAPTER 3. THE HARTLEY TRANSFORM (HT){}}}k := mh - ju := a[r+mh+j]v := a[r+mh+k]c := cos(j*PI/mh)s := sin(j*PI/mh){u, v} := {u*c+v*s, u*s-v*c}a[r+mh+j] := ua[r+mh+k] := v}for j:=0 to mh-1{u := a[r+j]v := a[r+j+mh]a[r+j]:= u + va[r+j+mh] := u - v}The derivation of the ‘usual’ DIT2 FHT algorithm starts by fusing the shift with the sum/diff step:void fht_localized_dit2(double *f, ulong ldn){const ulong n = 1UL<<ldn;revbin_permute(f, n);for (ulong ldm=1; ldm<=ldn; ++ldm){const ulong m = (1UL<<ldm);const ulong mh = (m>>1);const ulong m4 = (mh>>1);const double phi0 = M_PI/mh;for (ulong r=0; r<n; r+=m){{ // j == 0:ulong t1 = r;ulong t2 = t1 + mh;sumdiff(f[t1], f[t2]);}if ( m4 ){ulong t1 = r + m4;ulong t2 = t1 + mh;sumdiff(f[t1], f[t2]);}for (ulong j=1, k=mh-1; j<k; ++j,--k){double s, c;SinCos(phi0*j, &s, &c);ulong tj = r + mh + j;ulong tk = r + mh + k;double fj = f[tj];double fk = f[tk];f[tj] = fj * c + fk * s;f[tk] = fj * s - fk * c;ulong t1 = r + j;ulong t2 = tj; // == t1 + mh;sumdiff(f[t1], f[t2]);t1 = r + k;t2 = tk; // == t1 + mh;sumdiff(f[t1], f[t2]);}}}}57CHAPTER 3.

THE HARTLEY TRANSFORM (HT)58[FXT: fht localized dit2 in fht/fhtdit2.cc] Swapping the innermost loops then yields (considerationsas for DIT FFT, page 16, hold)void fht_dit2(double *f, ulong ldn)// decimation in time radix 2 fht{const ulong n = 1UL<<ldn;revbin_permute(f, n);for (ulong ldm=1; ldm<=ldn; ++ldm){const ulong m = (1UL<<ldm);const ulong mh = (m>>1);const ulong m4 = (mh>>1);const double phi0 = M_PI/mh;for (ulong r=0; r<n; r+=m){{ // j == 0:ulong t1 = r;ulong t2 = t1 + mh;sumdiff(f[t1], f[t2]);}if ( m4 ){ulong t1 = r + m4;ulong t2 = t1 + mh;sumdiff(f[t1], f[t2]);}}for (ulong j=1, k=mh-1; j<k; ++j,--k){double s, c;SinCos(phi0*j, &s, &c);for (ulong r=0; r<n; r+=m){ulong tj = r + mh + j;ulong tk = r + mh + k;double fj = f[tj];double fk = f[tk];f[tj] = fj * c + fk * s;f[tk] = fj * s - fk * c;ulong t1 = r + j;ulong t2 = tj; // == t1 + mh;sumdiff(f[t1], f[t2]);t1 = r + k;t2 = tk; // == t1 + mh;sumdiff(f[t1], f[t2]);}}}}[FXT: fht dit2 in fht/fhtdit2.cc]3.2.2Decimation in frequency (DIF) FHTIdea 3.2 (FHT radix 2 DIF step) Radix 2 decimation in frequency step for the FHT:ihn/2(even)H [a]= H a(lef t) + a(right)´i³hn/2(odd)H [a]= H X 1/2 a(lef t) − a(right)(3.9)(3.10)Code 3.4 (recursive radix 2 DIF FHT) Pseudo code for a recursive procedure of the (radix 2) DIFFHT algorithm:CHAPTER 3.

THE HARTLEY TRANSFORM (HT)procedure rec_fht_dif2(a[], n, x[])// real a[0..n-1] input// real x[0..n-1] result{real b[0..n/2-1], c[0..n/2-1]real s[0..n/2-1], t[0..n/2-1]if n == 1 then{x[0] := a[0]return}nh := n/2;for k:=0 to nh-1{s[k] := a[k]// ’left’t[k] := a[k+nh] // ’right’}for k:=0 to nh-1{{s[k], t[k]} := {s[k]+t[k],}hartley_shift(t[], nh, 1/2)rec_fht_dif2(s[], nh, b[])rec_fht_dif2(t[], nh, c[])j := 0for k:=0 to nh-1{x[j]:= b[k]x[j+1] := c[k]j := j+2}}59// workspace// workspaceelementselementss[k]-t[k]}[FXT: recursive fht dif2 in slow/recfht2.cc]Code 3.5 (radix 2 DIF FHT, localized) Pseudo code for a non-recursive procedure of the (radix 2)DIF FHT algorithm:procedure fht_dif2(a[], ldn)// real a[0..n-1] input,result{n := 2**ldn // length of a[] is a power of 2for ldm:=ldn to 1 step -1{m := 2**ldmmh := m/2m4 := m/4for r:=0 to n-m step m{for j:=0 to mh-1{u := a[r+j]v := a[r+j+mh]a[r+j]:= u + va[r+j+mh] := u - v}for j:=1 to m4-1{k := mh - ju := a[r+mh+j]v := a[r+mh+k]c := cos(j*PI/mh)s := sin(j*PI/mh){u, v} := {u*c+v*s, u*s-v*c}a[r+mh+j] := ua[r+mh+k] := vCHAPTER 3.