Arndt - Algorithms for Programmers (523138), страница 20
Текст из файла (страница 20)
Only the sign of the basisfunctions is shown. At the blank entries the functions are zero.template <typename Type>void transposed_haar_inplace_nn(Type *f, ulong ldn){ulong n = 1UL<<ldn;for (ulong js=n; js>=2; js>>=1){for (ulong j=0, t=js>>1; j<n; j+=js, t+=js){Type x = f[j];Type y = f[t];f[j] = x + y;f[t] = x - y;}}}CHAPTER 6. THE HAAR TRANSFORM113[FXT: transposed haar inplace nn in haar/transposedhaarnn.h]template <typename Type>void inverse_transposed_haar_inplace_nn(Type *f, ulong ldn){ulong n = 1UL<<ldn;for (ulong js=2; js<=n; js<<=1){for (ulong j=0, t=js>>1; j<n; j+=js, t+=js){Type x = f[j] * 0.5;Type y = f[t] * 0.5;f[j] = x + y;f[t] = x - y;;}}}6.40: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:The reversed Haar transform[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[+ + + + + + ++ - + - + - +++++++++++++++++ + + + + + +- + - + - + ++++++++-+ + + + + + ++ - + - + - ++++++++++++++ + + + + + +- + - + - + +++++++-+ + + + ]+ - + - ]+]+- ]]]]- ]]]]]]+]+]++- ]+]+]+]+]+]+]+]+]+]+]+]+]+]+]+]+- ]Figure 6.6: Basis functions for the reversed Haar transform.
Only the sign of the non-zero entries isshown.Let Hni denote the non-normalized in-place Haar transform (haar_inplace_nn), Let Htni denote thetransposed non-normalized in-place Haar transform (transposed_haar_inplace_nn), R the revbin-CHAPTER 6. THE HAAR TRANSFORM114permutation, H̄ the reversed Haar transform and H̄t the transposed reversed Haar transform. ThenH̄H̄tH̄ −1H̄t−1====R Hni RR Htni R−1R HniR−1R Htni R(6.3)(6.4)(6.5)(6.6)Code for the reversed Haar transform:template <typename Type>void haar_rev_nn(Type *f, ulong ldn){//const ulong n = (1UL<<ldn);for (ulong ldm=ldn; ldm>=1; --ldm){const ulong m = (1UL<<ldm);const ulong mh = (m>>1);ulong r = 0;//for (ulong r=0; r<n; r+=m) // almost walsh_wak_dif2(){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 + v;f[t2] = u - v;}}}}[FXT: haar rev nn in haar/haarrevnn.h] Note that this is almost the DIF implementation for theWalsh transform (Wk implemented as walsh_wak_dif2): The only thing that changed is that the linefor (ulong r=0; r<n; r+=m) was replaced by ulong r = 0.The transposed transform is:template <typename Type>void transposed_haar_rev_nn(Type *f, ulong ldn){for (ulong ldm=1; ldm<=ldn; ++ldm){const ulong m = (1UL<<ldm);const ulong mh = (m>>1);ulong r = 0;//for (ulong r=0; r<n; r+=m) // almost walsh_wak_dit2(){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 + v;f[t2] = u - v;}}}}[FXT: transposed haar rev nn in haar/transposedhaarrevnn.h] With the identical change as abovethis is almost the DIT implementation for the Walsh transform (Wk implemented as walsh_wak_dit2).The inverse transforms aretemplate <typename Type>CHAPTER 6.
THE HAAR TRANSFORM115void inverse_haar_rev_nn(Type *f, ulong ldn){for (ulong ldm=1; ldm<=ldn; ++ldm){const ulong m = (1UL<<ldm);const ulong mh = (m>>1);ulong r = 0;//for (ulong r=0; r<n; r+=m) // almost walsh_wak_dit2(){ulong t1 = r;ulong t2 = r + mh;for (ulong j=0; j<mh; ++j, ++t1, ++t2){Type u = f[t1] * 0.5;Type v = f[t2] * 0.5;f[t1] = u + v;f[t2] = u - v;}}}}andtemplate <typename Type>void inverse_transposed_haar_rev_nn(Type *f, ulong ldn){//const ulong n = (1UL<<ldn);for (ulong ldm=ldn; ldm>=1; --ldm){const ulong m = (1UL<<ldm);const ulong mh = (m>>1);ulong r = 0;//for (ulong r=0; r<n; r+=m) // almost walsh_wak_dif2(){ulong t1 = r;ulong t2 = r + mh;for (ulong j=0; j<mh; ++j, ++t1, ++t2){Type u = f[t1] * 0.5;Type v = f[t2] * 0.5;f[t1] = u + v;f[t2] = u - v;}}}}6.56.5.1Relations between Walsh- and Haar- transformsWalsh transforms from Haar transformsA length-n Walsh transform can be obtained from one length-n Haar transform, one transform of lengthnnnn2 , two transforms of length- 4 , four transforms of length- 8 , .
. . and 4 transforms of length-2. Using thereversed Haar transform the implementation is most straightforward: A Walsh transform (Wk , the onewith the Walsh Kronecker base) can be implemented as// algorithm WH1:ulong n = 1UL<<ldn;haar_rev_nn(f, ldn);for (ulong ldk=ldn-1; ldk>0; --ldk){ulong k = 1UL << ldk;for (ulong j=k; j<n; j+=2*k) haar_rev_nn(f+j, ldk);}The idea can be drawn symbolically as in figure 6.5.1.Equivalently, one can compute Wk using the transposed version:CHAPTER 6.
THE HAAR TRANSFORMHaar transforms:H(16)AAAAAAAAaaaaaaaaAAAAaaaaAAaaAa116H(8)BBBBbbbbBBbbBbH(4)CCccCcH(2)DdWalsh(16) =^= 1*H(16) + 1*H(8) + 2*H(4) + 4*H(2)AAAAAAAAaaaaaaaaAAAAaaaaBBBBbbbbAAaaCCccBBbbCCccAaDdCcDdBbDdCcDdFigure 6.7: Symbolic description of how to build a Walsh transform from Haar transforms.Transposed Haar transforms:H(16)H(8)AaAAaaBbAAAAaaaaBBbbAAAAAAAAaaaaaaaaBBBBbbbbH(4)H(2)CcCCccDdWalsh(16) =^= 1*H(16) + 1*H(8) + 2*H(4) + 4*H(2)AaDdCcDdBbDdCcDdAAaaCCccBBbbCCccAAAAaaaaBBBBbbbbAAAAAAAAaaaaaaaaFigure 6.8: Symbolic description of how to build a Walsh transform from Haar transforms, transposedversion.// algorithm WH1T:ulong n = 1UL<<ldn;for (ulong ldk=1; ldk<ldn; ++ldk){ulong k = 1UL << ldk;for (ulong j=k; j<n; j+=2*k) transposed_haar_rev_nn(f+j, ldk);}transposed_haar_rev_nn(f, ldn);The symbolic scheme is obtained by reversing the lines in the non-transposed scheme, this is shown infigure 6.5.1.Moreover, the inverse Walsh transform (Wk−1 =1nWk ) can be computed as// algorithm WH2T:ulong n = 1UL<<ldn;inverse_transposed_haar_rev_nn(f, ldn);for (ulong ldk=ldn-1; ldk>0; --ldk){ulong k = 1UL << ldk;for (ulong j=k; j<n; j+=2*k) inverse_transposed_haar_rev_nn(f+j, ldk);}or as// algorithm WH2:ulong n = 1UL<<ldn;for (ulong ldk=1; ldk<ldn; ++ldk)CHAPTER 6.
THE HAAR TRANSFORM{ulong k = 1UL << ldk;for (ulong j=k; j<n; j+=2*k)}inverse_haar_rev_nn(f, ldn);6.5.2inverse_haar_rev_nn(f+j, ldk);Haar transforms from Walsh transformsWalsh transform:W(16)AaDdCcDdBbDdCcDdAAaaCCccBBbbCCccAAAAaaaaBBBBbbbbAAAAAAAAaaaaaaaaInverse (or transposed) Walsh transforms:W(8):W(4):W(2):BBBBbbbbCCccDdBBbbCCccCcDdBbDdCcDdAaAAaaAAAAaaaaAAAAAAAAaaaaaaaaHaar(16)=^=BBBBbbbbCCccBBbbCCccDdCcDdBbDdCcDdAaDdCcDdBbDdCcDdAAaaCCccBBbbCCccAAAAaaaaBBBBbbbbAAAAAAAAaaaaaaaaW(16) + W(8) + W(4) + W(2)Figure 6.9: Symbolic description of how to build a Haar transform from Walsh transforms.The non-normalized transposed reversed Haar transform can (up to normalization) be obtained via// algorithm HW1: transposed_haar_rev_nn(f, ldn); =^=for (ulong ldk=1; ldk<ldn; ++ldk){ulong k = 1UL << ldk;walsh_wak(f+k, ldk);}walsh_wak(f, ldn);and its inverse as// algorithm HW1I: inverse_transposed_haar_rev_nn(f, ldn); =^=walsh_wak(f, ldn);for (ulong ldk=1; ldk<ldn; ++ldk){ulong k = 1UL << ldk;walsh_wak(f+k, ldk);}The non-normalized transposed Haar transform can (again, up to normalization) be obtained via// algorithm HW2: transposed_haar_nn(f, ldn); =^=for (ulong ldk=1; ldk<ldn; ++ldk){ulong k = 1UL << ldk;117CHAPTER 6.
THE HAAR TRANSFORMwalsh_pal(f+k, ldk);}walsh_pal(f, ldn);and its inverse as// algorithm HW2I: inverse_transposed_haar_nn(f, ldn); =^=walsh_pal(f, ldn); // =^= revbin_permute(f, n); walsh_wak(f, ldn);for (ulong ldk=1; ldk<ldn; ++ldk){ulong k = 1UL << ldk;walsh_pal(f+k, ldk);}The symbolic scheme is given in figure 6.5.2.6.6Integer to integer Haar transformCode 6.1 (integer to integer Haar transform)procedure int_haar(f[], ldn)// real f[0..2**ldn-1] // input, result{n := 2**nreal g[0..n-1] // workspacefor m:=n to 2 div_step 2{mh = m/2k := 0for j=0 to m-1 step 2{x := f[j]y := f[j+1]d := x - ys := y + floor(d/2) // == floor((x+y)/2)g[k]:= sg[mh+k] := dk := k + 1}copy g[0..m-1] to f[0..m-1]m := m/2}}Omit floor() with integer types.
[FXT: haar i2i in haar/haari2i.cc]Code 6.2 (inverse integer to integer Haar transform)procedure inverse_int_haar(f[], ldn)// real f[0..2**ldn-1] // input, result{n := 2**nreal g[0..n-1] // workspacefor m:=2 to n mul_step 2{mh := m/2k := 0for j=0 to m-1 step 2{s := f[k]d := f[mh+k]y := s - floor(d/2)118CHAPTER 6. THE HAAR TRANSFORMx := dg[j]g[j+1]k := k}}+ y // == s+floor((d+1)/2):= x:= y+ 1}copy g[0..m-1] to f[0..m-1]m := m * 2[FXT: inverse haar i2i in haar/haari2i.cc]119Chapter 7Permutations7.10:1:2:3:4:5:6:7:8:9:10:11:12:13:14:15:The revbin permutation[ *[*[*[*[*[*[*[*[*[*[*[*[*[*[*[*]]]]]]]]]]]]]]]]0:1:2:3:4:5:6:7:[ *[*[*[*[*[*[*[*]]]]]]]]0:1:2:3:[ *[*[*[*]]]]Figure 7.1: Permutation matrices of the revbin-permutation for sizes 16, 8 and 4.
The permutation isself-inverse.The procedure revbin_permute(a[], n) used in the DIF and DIT FFT algorithms rearranges the arraya[] in a way that each element ax is swapped with ax̃ , where x̃ is obtained from x by reversing its binarydigits. For example if n = 256 and x = 4310 = 001010112 then x̃ = 110101002 = 21210 . Note that x̃depends on both x and on n.7.1.1A naive versionA first implementation might look likeprocedure revbin_permute(a[], n)// a[0..n-1] input,result{for x:=0 to n-1{r := revbin(x, n)if r>x then swap(a[x], a[r])}}The condition r>x before the swap() statement makes sure that the swapping isn’t undone later whenthe loop variable x has the value of the present r. The function revbin(x, n) shall return the reversed120CHAPTER 7.
PERMUTATIONS0: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:[ *[*[*[*[*[*[*[*[*[*[*[*[*[*[*[*[*[*[*[*[*[*[*[*[*[*[*[*[*[*[*[*121]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]Figure 7.2: Permutation matrix of the revbin-permutation for size 32.bits of x:function revbin(x, n){j := 0ldn := log2(n) // is an integerwhile ldn>0{j := j << 1j := j + (x & 1)x := x >> 1ldn := ldn - 1}return j}This version of the revbin_permute-routine is pretty inefficient (even if revbin() is inlined and ldn isonly computed once). Each execution of revbin() costs proportional ldn operations, giving a total ofproportional n2 log2 (n) operations (neglecting the swaps for the moment).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
