Arndt - Algorithms for Programmers (523138), страница 18
Текст из файла (страница 18)
[FXT: slant inwalsh/slant.cc]The inverse transform is:void inverse_slant(double *f, ulong ldn)// inverse of slant(){ulong n = 1UL<<ldn;ulong ldm=ldn-2;do{ulong m = 1UL<<ldm; // m = n/4, n/2, ..., 4, 2, 1double N = m*2, N2 = N*N;double a = sqrt(3.0*N2/(4.0*N2-1.0));double b = sqrt(1.0-a*a); // == sqrt((N2-1)/(4*N2-1));for (ulong j=m; j<n-1; j+=4*m){ulong t1 = j;ulong t2 = j + m;double f1 = f[t1], f2 = f[t2];f[t1] = b * f2 + a * f1;f[t2] = a * f2 - b * f1;}}while ( ldm-- );walsh_wak(f, ldn);}A sequency ordered version of the transform can be implemented as follows:void slant_seq(double *f, ulong ldn)// sequency ordered slant transform{slant(f, ldn);ulong n = 1UL<<ldn;inverse_gray_permute(f, n);unzip_rev(f, n);revbin_permute(f, n);}This implementation could be optimized by fusing the involved permutations, cf.
[40].The inverse is trivially derived by calling the inverse operations in reversed order:void inverse_slant_seq(double *f, ulong ldn)// inverse of slant_seq(){ulong n = 1UL<<ldn;revbin_permute(f, n);zip_rev(f, n);gray_permute(f, n);inverse_slant(f, ldn);}TBD: figure: slant basis funcsCHAPTER 5. THE WALSH TRANSFORM AND ITS RELATIVES5.899The Reed-Muller transform (RMT)How to make a Reed-Muller transform out of a Walsh transform:‘Replace u+v by u and u-v by u XOR v, done.’0:1:2:3:4:5:6:7:8:9:10:11:12:13:14:15:16:17:18:19:20:21:22:23:24:25:26:27:28:29:30:31:[* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *][ ****************][* ** ** ** ** ** ** ** *][********][* * * ** * * ** * * ** * * *][********][* ** ** ** *][****][* * * * * * * ** * * * * * * *][********][* ** ** ** *][****][* * * ** * * *][****][* ** *][**][* * * * * * * * * * * * * * * *][********][* ** ** ** *][****][* * * ** * * *][****][* ** *][**][* * * * * * * *][****][* ** *][**][* * * *][**][* *][*]Figure 5.8: Basis functions for the Reed-Muller transform (RMT).
Asterisks denote positions where thefunctions are equal to one. Blank entries correspond to zero.There we go:template <typename Type>void word_reed_muller_dif2(Type *f, ulong ldn)// Reed-Muller Transform// Type must have the XOR operator.// decimation in frequency (DIF) algorithm// self-inverse{const ulong n = (1UL<<ldn);for (ulong ldm=ldn; ldm>=1; --ldm){const ulong m = (1UL<<ldm);const ulong mh = (m>>1);for (ulong r=0; r<n; r+=m){ulong t1 = r;ulong t2 = r+mh;for (ulong j=0; j<mh; ++j, ++t1, ++t2){Type u = f[t1];Type v = f[t2];f[t1] = u;f[t2] = u ^ v;}CHAPTER 5. THE WALSH TRANSFORM AND ITS RELATIVES}}100}[FXT: word reed muller dif2 in walsh/reedmuller.h].This routine is almost identical to [FXT: walsh wak dif2 in walsh/walshwak.h].
The only changes arewalshwalshreedmullerreedmullerf[t1]f[t2]f[t1]f[t2]====u + v;u - v;u;u ^ v;The decimation in time algorithm can be obtained from [FXT: walsh wak dit2 in walsh/walshwak.h]by the very same changes.As given, the transforms work word-wise, if the bit-wise transform is wanted usetemplate <typename Type>inline void bit_reed_muller(Type *f, ulong ldn){word_reed_muller_dif2(f, ldn);ulong n = 1UL << ldn;for (ulong k=0; k<n; ++k) f[k] = yellow_code(f[k]);}The yellow_code (see section 8.23) can also be applied before the main loop. In fact, the yellow code isthe Reed-Muller transform (RMT) on a binary word.The other ‘color-transforms’ of section 8.23 lead to variants of the RMT, the blue-code gives anotherself-inverse transform, the red-code and the cyan-code give transforms R and C so thatRRRCCCRC===idR−1 = R R = CidC −1 = C C = RC R = id(5.40)(5.41)(5.42)In fact, all relations given in the referenced section hold.As can be seen from the ‘atomic’ matrices (relations 8.40 .
. . 8.43) the four transforms corresponding tothe ‘color-codes’ are obtained byWalsh:B:Y:R:C:f[t1]f[t1]f[t1]f[t1]f[t1]=====u + v;u ^ v;u;v;u ^ v;f[t2]f[t2]f[t2]f[t2]f[t2]=====u - v;v;u ^ v;u ^ v;u;(Reed-Muller transform)The blue-code equivalent leads to a basis that is obtained by transposing the shown one, the red- andcyan- variants are the mutually inversered11111111111111111.1.1.1.1.1.1.1.11..11..11..11..1...1...1...1...1111....1111....1.1.....1.1.....11......11......1.......1.......11111111........1.1.1.1.........11..11..........1...1...........1111............1.1.............11..............1...............cyan...............1..............11.............1.1............1111...........1...1..........11..11.........1.1.1.1........11111111.......1.......1......11......11.....1.1.....1.1....1111....1111...1...1...1...1..11..11..11..11.1.1.1.1.1.1.1.11111111111111111CHAPTER 5.
THE WALSH TRANSFORM AND ITS RELATIVES101The symbolic powering idea from section 8.23 leads to transforms with bases (using eight element arrays):1........1........1........1........1........1........1........1x=01...1....1...1....1...1....1...1....1........1........1........1x=11.1......1.1......1........1........1.1......1.1......1........1x=21.1.1.1..1.1.1.1..1...1....1...1....1.1......1.1......1........1x=311.......1........11.......1........11.......1........11.......1x=411..11...1...1....11..11...1...1....11.......1........11.......1x=51111.....1.1......11.......1........1111.....1.1......11.......1x=611111111.1.1.1.1..11..11...1...1....1111.....1.1......11.......1x=7[FXT: file demo/bitxtransforms-demo.cc] gives the list for 32-bit words.[FXT: file aux1/wordgray.h] contains some functions that are the Gray code equivalents for GF (2n ):template <typename Type>void word_gray(Type *f, ulong ldn){ulong n = 1UL<<ldn;for (ulong k=0; k<n-1; ++k) f[k] ^= f[k+1];}template <typename Type>void inverse_word_gray(Type *f, ulong ldn){ulong n = 1UL<<ldn;for (ulong s=1; s<n; s*=2){// word_gray ** s:for (ulong k=0, j=k+s; j<n; ++k,++j)}}f[k] ^= f[j];As one might suspect, these are related to the Reed-Muller transform.
Writing Y (‘yellow’) for the RMT,g for the word- Gray code and Sk for the cyclic shift by k words (word zero is moved to position k) onehasY S+1 YY S−1 YY Sk Y===gg −1gk(5.43)(5.44)(5.45)These are exactly the relations 8.32 given on page 196 for the bit-wise transforms.For k >= 0 the operator Sk corresponds to the shift toward element zero (use [FXT: rotate sgn inperm/rotate.h]).The power of the word-wise Gray code is perfectly equivalent to the bit-wise version:template <typename Type>void word_gray_pow(Type *f, ulong ldn, ulong x){ulong n = 1UL<<ldn;x &= (n-1); // modulo nfor (ulong s=1; s<n; s*=2){if ( x & 1){// word_gray ** s:for (ulong k=0, j=k+s; j<n; ++k,++j)}x >>= 1;}}f[k] ^= f[j];With e be the green code operator, then for the ‘blue-variant’:B S+1 BB S−1 B==e−1e(5.46)(5.47)B Sk B=e−k(5.48)CHAPTER 5.
THE WALSH TRANSFORM AND ITS RELATIVES102Further,C Sk RC ek R==ekSk(5.49)(5.50)The transforms as Kronecker products (all operations are modulo two):·B2=·Y2=·R2=·C2=1011110101111110¸log2 (n)Bn =OB2(5.51)Y2(5.52)k=1¸log2 (n)Yn =Ok=1¸log2 (n)Rn =OR2(5.53)C2(5.54)k=1¸log2 (n)Cn =Ok=1A function that computes the k-th base function of the word-wise transform istemplate <typename Type>inline void word_reed_muller_basefunc(Type *f, ulong n, ulong k){for (ulong i=0; i<n; ++i){ulong x = (n-1-i) & (k);f[i] = ( 0==x ? +1 : 0 );}}[FXT: word reed muller basefunc in walsh/reedmuller.h]5.9The arithmetic transformHow to make a arithmetic transform out of a Walsh transform:‘Forward: replace u+v by u and u-v by v-u. Backward: replace u+v by u and u-v by u+v.’On to the code:template <typename Type>void arith_transform_plus(Type *f, ulong ldn)// Arithmetic Transform// Decimation In Frequency (DIF) algorithm{const ulong n = (1UL<<ldn);for (ulong ldm=ldn; ldm>=1; --ldm){const ulong m = (1UL<<ldm);const ulong mh = (m>>1);for (ulong r=0; r<n; r+=m){ulong t1 = r;ulong t2 = r+mh;for (ulong j=0; j<mh; ++j, ++t1, ++t2){Type u = f[t1];Type v = f[t2];f[t1] = u;f[t2] = u + v;CHAPTER 5.
THE WALSH TRANSFORM AND ITS RELATIVES0:1:2:3:4:5:6:7:8:9:10:11:12:13:14:15:16:17:18:19:20:21:22:23:24:25:26:27:28:29:30:31:[ + - - + - + + - - + + - + - - + - + + - + - - + + - - + - + + [++++++++[+ - +- ++ - ++ + - +[++++[+ - - +- + + - + + + - - +[++++[+ - +- ++ [++[+ - - + - + + - + + - + - - +[++++[+ - +- ++ [++[+ - - +- + + [++[+ - +[+[+ - - + - + + - - + + - + - - +[++++[+ - +- ++ [++[+ - - +- + + [++[+ - +[+[+ - - + - + + [++[+ - +[+[+ - - +[+[+ [+103]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]Figure 5.9: Basis functions for the Arithmetic transform (A− , the one with the minus sign). The valuesare ±1, blank entries denote 0.}}}}andtemplate <typename Type>void arith_transform_minus(Type *f, ulong ldn)// Inverse of arith_transform_plus{-- snip -f[t1] = u;f[t2] = v - u;-- snip -}The length-2 transforms can be written as·¸·¸·¸+1 0aa+A2 v ==+1 +1ba+b·¸·¸·¸+1 0aa−A2 v ==−1 +1bb−a(5.55)(5.56)That the one with the minus is called the forward transform is tradition.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
