Arndt - Algorithms for Programmers (523138), страница 19
Текст из файла (страница 19)
Similar to the Fourier transformwe avoid the forward- backward- naming scheme and put:template <typename Type>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:[ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +[++++++++++++++++[+ ++ ++ ++ ++ ++ ++ ++ +[++++++++[+ + + ++ + + ++ + + ++ + + +[++++++++[+ ++ ++ ++ +[++++[+ + + + + + + ++ + + + + + + +[++++++++[+ ++ ++ ++ +[++++[+ + + ++ + + +[++++[+ ++ +[++[+ + + + + + + + + + + + + + + +[++++++++[+ ++ ++ ++ +[++++[+ + + ++ + + +[++++[+ ++ +[++[+ + + + + + + +[++++[+ ++ +[++[+ + + +[++[+ +[+104]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]Figure 5.10: Basis functions for the inverse Arithmetic transform (A+ , the one without minus sign).
Thevalues are ±1, blank entries denote 0.inline void arith_transform(Type *f, ulong ldn, int is){if ( is>0 ) arith_transform_plus(f, ldn);else arith_transform_minus(f, ldn);}In Kronecker product notation the arithmetic transform and its inverse can be written as·A+2=·A−2=+1+10+1+1 0−1 +1¸log2 (n)A+n=OA+2(5.57)A−2(5.58)k=1¸log2 (n)A−n =Ok=1Chapter 6The Haar transform0: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:[+[+[+[[+[[[[+[[[[[[[[+[[[[[[[[[[[[[[[+ + + + + + + + + + + + + + + + ++ + + + + + + + + + + + + + + - + + + + + + + - - - - - - - + ++ + + - - - + + + + - - - + ++ + + + + + + + + + + + + +]- - - - - - - - - - - - - -]]+ + + + + + - - - - - - - -]]]+ + - - - ]+ + + + - - - -]+ - ]+ + - ]+ + - ]+ + - ]+ + - ]+ + - ]+ + - ]+ + - -]]+ ]+ ]+ ]+ ]+ ]+ ]+ ]+ ]+ ]+ ]+ ]+ ]+ ]+ ]+ -]1/sqrt(32)1/sqrt(32)1/sqrt(16)1/sqrt(16)1/sqrt(8)1/sqrt(8)1/sqrt(8)1/sqrt(8)1/sqrt(4)1/sqrt(4)1/sqrt(4)1/sqrt(4)1/sqrt(4)1/sqrt(4)1/sqrt(4)1/sqrt(4)1/sqrt(2)1/sqrt(2)1/sqrt(2)1/sqrt(2)1/sqrt(2)1/sqrt(2)1/sqrt(2)1/sqrt(2)1/sqrt(2)1/sqrt(2)1/sqrt(2)1/sqrt(2)1/sqrt(2)1/sqrt(2)1/sqrt(2)1/sqrt(2)Figure 6.1: Basis functions for the Haar transform.
Only the sign of the basis functions is shown. Theabsolute value of the non-zero entries in each row is given at the right. The norm of each row is one. Atthe blank entries the functions are zero.Code for the Haar transform:template <typename Type>void haar(Type *f, ulong ldn, Type *ws=0){ulong n = (1UL<<ldn);Type s2 = sqrt(0.5);Type v = 1.0;Type *g = ws;105CHAPTER 6. THE HAAR TRANSFORM}106if ( !ws ) g = new Type[n];for (ulong m=n; m>1; m>>=1){v *= s2;ulong mh = (m>>1);for (ulong j=0, k=0; j<m; j+=2, k++){Type x = f[j];Type y = f[j+1];g[k]= x + y;g[mh+k] = (x - y) * v;}copy(g, f, m);}f[0] *= v; // v == 1.0/sqrt(n);if ( !ws ) delete [] g;The above routine uses a temporary workspace that can be supplied by the caller.
The computationalcost is only ∼ n. [FXT: haar in haar/haar.h]Code for the inverse Haar transform:template <typename Type>void inverse_haar(Type *f, ulong ldn, Type *ws=0){ulong n = (1UL<<ldn);Type s2 = sqrt(2.0);Type v = 1.0/sqrt(n);Type *g = ws;if ( !ws ) g = new Type[n];f[0] *= v;for (ulong m=2; m<=n; m<<=1){ulong mh = (m>>1);for (ulong j=0, k=0; j<m; j+=2, k++){Type x = f[k];Type y = f[mh+k] * v;g[j]= x + y;g[j+1] = x - y;}copy(g, f, m);v *= s2;}if ( !ws ) delete [] g;}[FXT: inverse haar in haar/haar.h]That the given routines use a temporary storage may be seen as a disadvantage. A rather simplereordering of the basis functions, however, allows for to an in-place algorithm.
This leads to the in-placeHaar transform.6.1In-place Haar transformCode for the in-place version of the Haar transform:template <typename Type>void haar_inplace(Type *f, ulong ldn){ulong n = 1UL<<ldn;Type s2 = sqrt(0.5);Type v = 1.0;for (ulong js=2; js<=n; js<<=1){v *= s2;for (ulong j=0, t=js>>1; j<n;j+=js, t+=js)CHAPTER 6. THE HAAR TRANSFORM0: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:[+[+[+[[+[[[[+[[[[[[[[+[[[[[[[[[[[[[[[+ + + + + + + + + + + + + ++ - + + + + - - - + + + - + + + + + + + + - - - - - - + + + - + + + + + - - + + + ++ + + + + + + + + + + + + ++ + + + + + + + + + + + + + + + +]]]]]]]]]]]]]]]]+ - - - - - - - - - - - - - - - -]+ ]+ + - ]+ ]+ + + + - - - ]+ ]+ + - ]+ ]+ + + + + + + + - - - - - - - -]+ ]+ + - ]+ ]+ + + + - - - -]+ ]+ + - -]+ -]1071/sqrt(32)1/sqrt(4)1/sqrt(4)1/sqrt(4)1/sqrt(8)1/sqrt(4)1/sqrt(4)1/sqrt(4)1/sqrt(16)1/sqrt(4)1/sqrt(4)1/sqrt(4)1/sqrt(8)1/sqrt(4)1/sqrt(4)1/sqrt(4)1/sqrt(32)1/sqrt(4)1/sqrt(4)1/sqrt(4)1/sqrt(8)1/sqrt(4)1/sqrt(4)1/sqrt(4)1/sqrt(16)1/sqrt(4)1/sqrt(4)1/sqrt(4)1/sqrt(8)1/sqrt(4)1/sqrt(4)1/sqrt(4)Figure 6.2: Haar basis functions, in-place order.
Only the sign of the basis functions is shown. Theabsolute value of the non-zero entries in each row is given at the right. The norm of each row is one. Atthe blank entries the functions are zero.{}Type x = f[j];Type y = f[t];f[j] = x + y;f[t] = (x - y) * v;}}f[0] *= v; // v==1.0/sqrt(n);[FXT: haar inplace in haar/haar.h]. . . and its inverse:template <typename Type>void inverse_haar_inplace(Type *f, ulong ldn){ulong n = 1UL<<ldn;Type s2 = sqrt(2.0);Type v = 1.0/sqrt(n);f[0] *= v;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] * v;f[j] = x + y;f[t] = x - y;}CHAPTER 6. THE HAAR TRANSFORM0: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:[+[+[+[[+[[[[+[[[[[[[[+[[[[[[[[[[[[[[[108+ + + + + + + + + + + + + + + + ++ + + + + + + + + + + + + + + - + + + + + + + - - - - - - - + ++ + + - - - + ++ + + + - - - -+ + + + + + + + + + + + + +]- - - - - - - - - - - - - -]]+ + + + + + - - - - - - - -]]+ + - - - ]]+ + + + - - - -]+ - ]+ + - ]+ + - ]+ + - ]+ + - ]+ + - ]+ + - ]+ + - -]]+ ]+ ]+ ]+ ]+ ]+ ]+ ]+ ]+ ]+ ]+ ]+ ]+ ]+ ]+ -]1/sqrt(32)1/sqrt(32)1/sqrt(16)1/sqrt(16)1/sqrt(8)1/sqrt(8)1/sqrt(8)1/sqrt(8)1/sqrt(4)1/sqrt(4)1/sqrt(4)1/sqrt(4)1/sqrt(4)1/sqrt(4)1/sqrt(4)1/sqrt(4)1/sqrt(2)1/sqrt(2)1/sqrt(2)1/sqrt(2)1/sqrt(2)1/sqrt(2)1/sqrt(2)1/sqrt(2)1/sqrt(2)1/sqrt(2)1/sqrt(2)1/sqrt(2)1/sqrt(2)1/sqrt(2)1/sqrt(2)1/sqrt(2)Figure 6.3: Basis functions of the in-place order Haar transform followed by a revbin permutation.
Notethat the ordering is such that basis functions that are identical up to a shift appear consecutively.}}v *= s2;[FXT: inverse haar inplace in haar/haar.h]The in-place Haar transform Hi is related to the ‘usual’ Haar transform H by a permutation PH via therelationsHH−1==P H · Hi−1Hi−1 · PHPH can be programmed astemplate <typename Type>void haar_permute(Type *f, ulong ldn){revbin_permute(f, 1UL<<ldn);for (ulong ldm=1; ldm<=ldn-1; ++ldm){ulong m = (1UL<<ldm); // m=2, 4, 8, ..., n/2revbin_permute(f+m, m);}}while its inverse istemplate <typename Type>(6.1)(6.2)CHAPTER 6. THE HAAR TRANSFORM109void inverse_haar_permute(Type *f, ulong ldn){for (ulong ldm=1; ldm<=ldn-1; ++ldm){ulong m = (1UL<<ldm); // m=2, 4, 8, ..., n/2revbin_permute(f+m, m);}revbin_permute(f, 1UL<<ldn);}(cf. [FXT: file perm/haarpermute.h])Then, as given above, haar is equivalent toinplace_haar();haar_permute();and inverse_haar is equivalent toinverse_haar_permute();inverse_inplace_haar();6.2Non-normalized Haar transformsVersions of the Haar transform without normalization are given in [FXT: file haar/haarnn.h].
The basisfunctions are the same as for the normalized versions, only the absolute value of the non-zero entries aredifferent.template <typename Type>void haar_nn(Type *f, ulong ldn, Type *ws=0){ulong n = (1UL<<ldn);Type *g = ws;if ( !ws ) g = new Type[n];for (ulong m=n; m>1; m>>=1){ulong mh = (m>>1);for (ulong j=0, k=0; j<m; j+=2, k++){Type x = f[j];Type y = f[j+1];g[k]= x + y;g[mh+k] = x - y;}copy(g, f, m);}if ( !ws ) delete [] g;}[FXT: haar nn in haar/haarnn.h]template <typename Type>void inverse_haar_nn(Type *f, ulong ldn, Type *ws=0){ulong n = (1UL<<ldn);Type s2 = 2.0;Type v = 1.0/n;Type *g = ws;if ( !ws ) g = new Type[n];f[0] *= v;for (ulong m=2; m<=n; m<<=1){ulong mh = (m>>1);for (ulong j=0, k=0; j<m; j+=2, k++){CHAPTER 6. THE HAAR TRANSFORMType x = f[k];Type y = f[mh+k] * v;g[j]= x + y;g[j+1] = x - y;}copy(g, f, m);v *= s2;}}if ( !ws )delete [] g;template <typename Type>void haar_inplace_nn(Type *f, ulong ldn)//// transform wrt.
to non-normalized haar base// in-place operation//// the sequence//haar_inplace_nn(); haar_permute();// is equivalent to//haar_nn()//{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];Type y = f[t];f[j] = x + y;f[t] = x - y;}}}[FXT: haar inplace nn in haar/haarnn.h]template <typename Type>void inverse_haar_inplace_nn(Type *f, ulong ldn)//// inverse transform wrt. to haar base// in-place operation//// the sequence//inverse_haar_permute();//inverse_haar_inplace();// is equivalent to//inverse_haar()//{ulong n = 1UL<<ldn;Type s2 = 2.0;Type v = 1.0/n;f[0] *= v;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] * v;f[j] = x + y;f[t] = x - y;}v *= s2;}}110CHAPTER 6.
THE HAAR TRANSFORM0: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:[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[++++++++++++++++++++++++++++++++++++++++++++++++-++++++++-++++-++++++++-++++++++-+-++++++++-++++++++-++++-111]]+]]+]]+]]+]]+]]+]]+]]+]]+]]+]]+]]+]]+]]+]]+ ]- ]Figure 6.4: Basis functions for the transposed Haar transform. Only the sign of the basis functions isshown.
At the blank entries the functions are zero.6.3Transposed Haar transformstemplate <typename Type>void transposed_haar_nn(Type *f, ulong ldn, Type *ws=0){ulong n = (1UL<<ldn);Type *g = ws;if ( !ws ) g = new Type[n];for (ulong m=2; m<=n; m<<=1){ulong mh = (m>>1);for (ulong j=0, k=0; j<m; j+=2, k++){Type x = f[k];Type y = f[mh+k];g[j]= x + y;g[j+1] = x - y;}copy(g, f, m);}if ( !ws ) delete [] g;}[FXT: transposed haar nn in haar/transposedhaarnn.h]template <typename Type>void inverse_transposed_haar_nn(Type *f, ulong ldn, Type *ws=0){ulong n = (1UL<<ldn);Type *g = ws;CHAPTER 6. THE HAAR TRANSFORM}if ( !ws ) g = new Type[n];for (ulong m=n; m>1; m>>=1){ulong mh = (m>>1);for (ulong j=0, k=0; j<m;{Type x = f[j]* 0.5;Type y = f[j+1] * 0.5;g[k]= x + y;g[mh+k] = x - y;}copy(g, f, m);}if ( !ws ) delete [] g;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:[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[112j+=2, k++)+ + ++++]+ - ++++]+- + +++]+- - +++]+- + +++]+- - +++]+- + ++]+- - ++]+- + +++]+- - +++]+- + ++]+- - ++]+- + ++]+- - ++]+- + +]+- - +]+- + +++]+- - +++]+- + ++]+- - ++]+- + ++]+- - ++]+- + +]+- - +]+- + ++]+- - ++]+- + +]+- - +]+- + +]+- - +]+- + ]+- - ]Figure 6.5: Basis functions for the transposed in-place Haar transform.