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

THE HARTLEY TRANSFORM (HT)}}}}revbin_permute(a[], n)[FXT: fht localized dif2 in fht/fhtdif2.cc]The ‘usual’ DIF2 FHT algorithm then isvoid fht_dif2(double *f, ulong ldn)// decimation in frequency radix 2 fht{const ulong n = (1UL<<ldn);for (ulong ldm=ldn; ldm>=1; --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;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]);double fj = f[tj];double fk = f[tk];f[tj] = fj * c + fk * s;f[tk] = fj * s - fk * c;}}}revbin_permute(f, n);}[FXT: fht dif2 in fht/fhtdif2.cc]TBD: higher radix FHT60CHAPTER 3.

THE HARTLEY TRANSFORM (HT)3.361Complex FT by HTThe relations between the HT and the FT can be read off directly from their definitions and theirsymmetry relations. Let σ be the sign of the exponent in the FT, then the HT of a complex sequenced ∈ C is:³´´1 ³H [d] + H [d] + σ i H [d] − H [d]F [d] =(3.11)2Written out for the real and imaginary part d = a + i b (a, b ∈ R):<F [a + i b]==F [a + i b]=³´´1 ³H [a] + H [a] − σ H [b] − H [b]2³´´1 ³H [b] + H [b] + σ H [a] − H [a]2(3.12)(3.13)Alternatively, one can recast the relations (using the symmetry relations 3.5 and 3.6) as<F [a + i b] ==F [a + i b] =1H [aS − σ bA ]21H [bS + σ aA ]2(3.14)(3.15)Both formulations lead to the very sameCode 3.6 (complex FT by HT conversion)fht_fft_conversion(a[],b[],n,is)// preprocessing to use two length-n FHTs// to compute a length-n complex FFT// or// postprocessing to use two length-n FHTs// to compute a length-n complex FFT//// self-inverse{for k:=1 to n/2-1{t := n-kas := a[k] + a[t]aa := a[k] - a[t]bs := b[k] + b[t]ba := b[k] - b[t]aa := is * aaba := is * baa[k] := 1/2 * (as - ba)a[t] := 1/2 * (as + ba)b[k] := 1/2 * (bs + aa)b[t] := 1/2 * (bs - aa)}}[FXT: fht fft conversion in fht/fhtfft.cc] [FXT: fht fft conversion in fht/fhtcfft.cc]Now we have two options to compute a complex FT by two HTs:Code 3.7 (complex FT by HT, version 1) Pseudo code for the complex Fourier transform that usesthe Hartley transform, is must be -1 or +1:fft_by_fht1(a[],b[],n,is)// real a[0..n-1] input,result (real part)// real b[0..n-1] input,result (imaginary part)CHAPTER 3.

THE HARTLEY TRANSFORM (HT){62fht(a[], n)fht(b[], n)fht_fft_conversion(a[], b[], n, is)}andCode 3.8 (complex FT by HT, version 2) Pseudo code for the complex Fourier transform that usesthe Hartley transform, is must be -1 or +1:fft_by_fht2(a[],b[],n,is)// real a[0..n-1] input,result (real part)// real b[0..n-1] input,result (imaginary part){fht_fft_conversion(a[], b[], n, is)fht(a[], n)fht(b[], n)}Note that the real and imaginary parts of the FT are computed independently by this procedure.For convolutions it would be sensible to use procedure 3.7 for the forward and 3.8 for the backwardtransform. The complex squarings are then combined with the pre- and post-processing steps, therebyinterleaving the most nonlocal memory accesses with several arithmetic operations.[FXT: fht fft in fht/fhtcfft.cc]3.4Complex FT by complex HT and vice versaA complex valued HT is simply two HTs (one of the real, one of the imaginary part).

So we can useboth of 3.7 or 3.8 and there is nothing new. Really? If one writes a type complex version of both theconversion and the FHT the routine 3.7 will look likefft_by_fht1(c[], n, is)// complex c[0..n-1] input,result{fht(c[], n)fht_fft_conversion(c[], n, is)}(the 3.8 equivalent is hopefully obvious)This may not make you scream but here is the message: it makes sense to do so. It is pretty easy toderive a complex FHT from the real (i.e.

usual) version1 and with a well optimized FHT you get an evenbetter optimized FFT. Note that this trivial rewrite virtually gets you a length-n FHT with the bookkeeping and trig-computation overhead of a length-n/2 FHT.[FXT: fht dit core in fht/cfhtdit.cc][FXT: fht dif core in fht/cfhtdif.cc][FXT: fht fft conversion in fht/fhtcfft.cc][FXT: fht fft in fht/fhtcfft.cc]Vice versa: Let T be the operator corresponding to the fht_fft_conversion, T is its own inverse:T = T −1 , or, equivalently T · T = 1. We have seen thatF =H·T1 infact this is done automatically in FXTandF =T ·H(3.16)CHAPTER 3. THE HARTLEY TRANSFORM (HT)63Therefore triviallyH=T ·FandH=F ·T(3.17)Hence we have eitherfht_by_fft(c[], n, is)// complex c[0..n-1] input,result{fft(c[], n)fht_fft_conversion(c[], n, is)}or the same thing with swapped lines.

Of course the same ideas also work for separate real- and imaginaryparts.3.5Real FT by HT and vice versaTo express the real and imaginary part of a Fourier transform of a purely real sequence a ∈ R by itsHartley transform use relations 3.12 and 3.13 and set b = 0:<F [a] ==F [a] =1(H [a] + H [a])21(H [a] − H [a])2(3.18)(3.19)The pseudo code is straight forward:Code 3.9 (real to complex FFT via FHT)procedure real_complex_fft_by_fht(a[], n)// real a[0..n-1] input,result{fht(a[], n)for i:=1 to n/2-1{t := n - iu := a[i]v := a[t]a[i] := 1/2 * (u+v)a[t] := 1/2 * (u-v)}}At the end of this procedure the ordering of the output data c ∈ C isa[0]a[1]a[2]===<c0<c1<c2a[n/2]a[n/2 + 1]a[n/2 + 2]a[n/2 + 3]...====...<cn/2=cn/2−1=cn/2−2=cn/2−3a[n − 1]==c1[FXT: fht real complex fft in realfft/realfftbyfht.cc]The inverse procedure is:(3.20)CHAPTER 3.

THE HARTLEY TRANSFORM (HT)64Code 3.10 (complex to real FFT via FHT)procedure complex_real_fft_by_fht(a[], n)// real a[0..n-1] input,result{for i:=1 to n/2-1{t := n - iu := a[i]v := a[t]a[i] := u+va[t] := u-v}fht(a[], n)}[FXT: fht complex real fft in realfft/realfftbyfht.cc]Vice versa: same line of thought as for complex versions. Let Trc be the operator corresponding to the post-processing in real_complex_fft_by_fht, and Tcr correspond to the preprocessing incomplex_real_fft_by_fht. That isFcr = H · TcrandFrc = Trc · H(3.21)−1−1It should be no surprise that Trc · Tcr = 1, or, equivalently Trc = Tcrand Tcr = Trc. ThereforeH = Tcr · FrcandH = Fcr · Trc(3.22)The corresponding code should be obvious. Watch out for real/complex FFTs that use a different orderingthan 3.20.3.6Discrete cosine transform (DCT) by HTThe discrete cosine transform wrt.

the basisu(k) =ν(k) · cosπ k (i + 1/2)n(3.23)√(where ν(k) = 1 for k = 0, ν(k) = 2 else) can be computed from the FHT using an auxiliary routinenamed cos_rot. TBD: give cosrot’s action mathematicallyprocedure cos_rot(x[], y[], n)// real x[0..n-1] input// real y[0..n-1] result{nh := n/2x[0] := y[0]x[nh] := y[nh]phi := PI/2/nfor (ulong k:=1; k<nh; k++){c := cos(phi*k)s := sin(phi*k)cps := (c+s)*sqrt(1/2)cms := (c-s)*sqrt(1/2)x[k]:= cms*y[k] + cps*y[n-k]x[n-k] := cps*y[k] - cms*y[n-k]}}which is its own inverse.

(cf. [FXT: cos rot in dctdst/cosrot.cc]) ThenCHAPTER 3. THE HARTLEY TRANSFORM (HT)65Code 3.11 (DCT via FHT) Pseudo code for the computation of the DCT via FHT:procedure dcth(x[], ldn)// real x[0..n-1] input,result{n := 2**nreal y[0..n-1] // workspaceunzip_rev(x, y, n)fht(y[],ldn)cos_rot(y[], x[], n)}whereprocedure unzip_rev(a[], b[], n)// real a[0..n-1] input// real b[0..n-1] result{nh := n/2for k:=0 to nh-1{k2 := 2*kb[k]:= a[k2]b[nh+k] := a[n-1-k2]}}(see section 7.6, page 131).The inverse routine isCode 3.12 (IDCT via FHT) Pseudo code for the computation of the IDCT via FHT:procedure idcth(x[], ldn)// real x[0..n-1] input,result{n := 2**nreal y[0..n-1] // workspacecos_rot(x[], y[], n);fht(y[],ldn)zip_rev(y[], x[], n)}whereprocedure zip_rev(a[], b[], n)// real a[0..n-1] input// real b[0..n-1] result{nh := n/2for k:=0 to nh-1{k2 := 2*kb[k]:= a[k2]b[nh+k] := a[n-1-k2]}}The implementation of both the forward and the backward transform (cf.

[FXT: dcth and idcth indctdst/dcth.cc]) avoids the temporary array y[] if no scratch space is supplied.Cf. [32], [33].An alternative variant for the computation of the DCT that also uses the FHT is given in [FXT:dcth zapata in dctdst/dctzapata.cc]. The algorithm is described in [34].3.7Discrete sine transform (DST) by DCTTBD: definition dst, idstCHAPTER 3.