Arndt - Algorithms for Programmers (523138), страница 17
Текст из файла (страница 17)
If one wants the basis functionsordered by their sequency one can use the procedureconst ulong n = (1UL<<ldn);walsh_wak_dif2(f, ldn);revbin_permute(f, n);inverse_gray_permute(f, n);That isWw=G−1 R Wk = Wk R GA function that computes the k-th base function of the transform istemplate <typename Type>void walsh_wal_basefunc(Type *f, ulong n, ulong k){k = revbin(k, ld(n)+1);k = gray_code(k);//// =^=//k = revbin(k, ld(n));//k = green_code(k);for (ulong i=0; i<n; ++i){ulong x = i & k;x = parity(x);f[i] = ( 0==x ? +1 : -1 );}}[FXT: walsh wal basefunc in walsh/walshbasefunc.h]A version of the transform that avoids the Gray-permutation is based on1 Notethat the sequency of a signal with frequency f usually is 2 f .(5.31)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:[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[********************************************************************************************************************************* * * * * * ** * * * * * **** ** ** * * ** * * ** * ** * *** * * *** * * ** ** ** ** ** ** ** ** ** ** *** **** *** *** ********* **** ****** **** ************** * * * * * * * * * ** *** * * * * * * ** * * * * ** ** * * * ** * * ** * * *** * * ** * ** ** * * ** * * *** ** ** * * ** ** ** ** * ** ** ** ** *** ** ** *** ** *** ** ***** *** *** *** *** *** ******* *** **** **** ***** ********* ********* ************* * * * * * * ]]* * * * * * * ]]* * * * ]* * *]* * * * ]* * *]** * ]* * * *]** * ]* * * *]* ** * ]** *]* ** * ]** *]* ** ]* ** *]* ** ]* ** *]* *** ]** *]* *** ]** *]*** ]** **]*** ]** **]**** ]***]**** ]***]92( 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.4: Basis functions for the sequency-ordered Walsh transform (Walsh-Kacmarz basis).
Asterisksdenote the value +1, blank entries denote −1.template <typename Type>void walsh_wal_dif2_core(Type *f, ulong ldn)// decimation in frequency (DIF) algorithm// gray_permute is absorbed//// walsh_wal(f, ldn)//=^=// revbin_permute(f, n); walsh_wal_dif2_core(f, ldn);{const ulong n = (1UL<<ldn);for (ulong ldm=ldn; ldm>=2; --ldm){const ulong m = (1UL<<ldm);const ulong mh = (m>>1);const ulong m4 = (mh>>1);for (ulong r=0; r<n; r+=m){ulong j;for (j=0; j<m4; ++j){ulong t1 = r+j;ulong t2 = t1+mh;double u = f[t1];double v = f[t2];f[t1] = u + v;f[t2] = u - v;}for ( ; j<mh; ++j){ulong t1 = r+j;ulong t2 = t1+mh;CHAPTER 5.
THE WALSH TRANSFORM AND ITS RELATIVES}//}}93double u = f[t1];double v = f[t2];f[t1] = u + v;f[t2] = v - u; // reversed}if ( ldn ){// ulong ldm=1;const ulong m = 2; //(1UL<<ldm);const ulong mh = 1; //(m>>1);for (ulong r=0; r<n; r+=m){ulong j = 0;for (ulong j=0; j<mh; ++j){ulong t1 = r+j;ulong t2 = t1+mh;double u = f[t1];double v = f[t2];f[t1] = u + v;f[t2] = u - v;}}}The transform still needs the revbin-permutation:template <typename Type>inline void walsh_wal(Type *f, ulong ldn){revbin_permute(f, (1UL<<ldn));walsh_wal_dif2_core(f, ldn);// =^=//walsh_wal_dit2_core(f, ldn);//revbin_permute(f, (1UL<<ldn));}A decimation in time (DIT) version of the core-routine is also given in [FXT: file walsh/walshwal.h].The procedure gray_permute() is introduced in section 7.8.For the variantwalsh_gray(f, ldn);grs_negate(f, n);revbin_permute(f, n);is equivalent to walsh wal(f, ldn) and might be faster for large arrays.
We haveWw=R Q Wg = Wg−1 R Q(5.32)An alternative ordering of the base functions (first even sequencies ascending then odd sequencies descending [FXT: walsh wal rev in walsh/walshwalrev.h]) can be obtained by either ofrevbin_permute(f, 1UL<<ldn);gray_permute(f, n);walsh_wak(f,ldn);walsh_wak(f,ldn);inverse_gray_permute(f, n);revbin_permute(f, 1UL<<ldn);zip_rev(f, n);walsh_wal(f, ldn);walsh_wal(f, ldn);unzip_rev(f, n);walsh_wak(f, ldn);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:[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[***************************************************************** * * * * ** * * * ***** * ** ** ** *** ** ** **** *** **** ***************** *** **** ** *** ** ** *** ** ** * ***** * * * ** * * * * ** * * * * * * * * * * * * * * *** * * * * * * ** * ** * * *** * * ** * ** * * ** * * *** ** * * ** *** ** ** ** ** ** ** ** ** ** ** *** *** *** *** ** ****** **** *** **** ******* *************************** *** *********** **** **** *** **** *** ** ** *** *** *** ** *** ** ** *** ** ** ** * ** * * ** ** *** ** * * ** * ** * * *** * * ** * * * ** * * * * * * ** * * * * * ** * * * ** * * * ** ** **** ** ** ** ** *** ****************** *** ****** *** **** *** ** * * ** * * ** * ** * *********************************94]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]]( 0)( 2)( 4)( 6)( 8)(10)(12)(14)(16)(18)(20)(22)(24)(26)(28)(30)(31)(29)(27)(25)(23)(21)(19)(17)(15)(13)(11)( 9)( 7)( 5)( 3)( 1)Figure 5.5: Basis functions for the reversed sequency ordered Walsh transform.
Asterisks denote thevalue +1, blank entries denote −1.inverse_gray_permute(f, n);revbin_permute(f, n);revbin_permute(f, 1UL<<ldn);walsh_gray(f, ldn);grs_negate(f, n);That is,W¯w===Wk G R = R G−1 WkWw Z̄ = Z̄ −1 WwQ Wg R(5.33)(5.34)(5.35)However, an implementation that is more efficient uses the core-routines that have the Gray-permutation‘absorbed’:template <typename Type>inline void walsh_wal_rev(Type *f, ulong ldn){revbin_permute(f, (1UL<<ldn));walsh_wal_dit2_core(f, ldn);// =^=//walsh_wal_dif2_core(f, ldn);//revbin_permute(f, n);}[FXT: walsh wal rev in walsh/walshwalrev.h].CHAPTER 5.
THE WALSH TRANSFORM AND ITS RELATIVES95This implementation uses the fact thatW¯w= R Ww R(5.36)We can do still better:revbin_permute(f, n);walsh_gray(f, ldn);grs_negate(f, n);gives the same transform.Similar relations as for the transform with Walsh-Paley basis (5.29 and 5.30) hold for Ww :Ww==G Ww G = G−1 Ww G−1Z̄ Ww Z̄ = Z̄ −1 Ww Z̄ −1(5.37)(5.38)The k-th base function of the transform can be computed astemplate <typename Type>void walsh_wal_rev_basefunc(Type *f, ulong n, ulong k){k = revbin(k, ld(n));k = gray_code(k);// =^=//k = green_code(k);//k = revbin(k, ld(n));for (ulong i=0; i<n; ++i){ulong x = i & k;x = parity(x);f[i] = ( 0==x ? +1 : -1 );}}[FXT: walsh wal rev basefunc in walsh/walshbasefunc.h]The next variant of the Walsh transform has the interesting feature that the basis functions for a length-ntransform have only sequencies n/2 and n/2 − 1 at the even and odd indices, respectively. The transformis self-inverse (the basis is shown in figure 5.6) and can be obtained viatemplate <typename Type>void walsh_q1(Type *f, ulong ldn){ulong n = 1UL << ldn;grs_negate(f, n);walsh_gray(f, ldn);revbin_permute(f, n);}[FXT: walsh q1 in walsh/walshq.h]A different transform with sequency n/2 for the first half of the basis, sequency n/2 − 1 for the secondhalf ([FXT: walsh q2 in walsh/walshq.h], basis shown in figure 5.7) is computed bytemplate <typename Type>void walsh_q2(Type *f, ulong ldn){ulong n = 1UL << ldn;revbin_permute(f, n);grs_negate(f, n);walsh_gray(f, ldn);// =^=//grs_negate(f, n);//revbin_permute(f, n);//walsh_gray(f, ldn);}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:[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[* * ** * ** * ** ** ** *** * ** * ** * ** ********* **** ** * ********* * ******* * * *********** **** * *** * ***** * ****** * * ** * * ***** * ********* * ** * *** * ****** *** * *** * **** ***** * ** * ** ** * *** * ** * ** ** * ****** ** ***** **** **** ** * ******* ***** * *** ** ** * *** * *** * *** *** * ****** ********* * ** * ** * * **** * * ******* * * ******* * ** ** ************* * ** * ** * ** *** * ******* * * *********** **** * ** * ****** ***** * ** * ** * *** * ** * ** ** * ** * ** ******* ****** *** * ** * *** ****** * ** * *** ********* * * ***** * ***** * **** * *** * ***** *** *** *** * ** ** * ** * ** ** ***** ]]* ]* ]* ]]]]* ]]* ]* ]]* ]* ]* ]* ]]* ]* ]* ]]]]]* ]]]* ]]]]96(16)(15)(16)(15)(16)(15)(16)(15)(16)(15)(16)(15)(16)(15)(16)(15)(16)(15)(16)(15)(16)(15)(16)(15)(16)(15)(16)(15)(16)(15)(16)(15)Figure 5.6: Basis functions for a self-inverse Walsh transform that has sequencies n/2 and n/2 − 1 only.Asterisks denote the value +1, blank entries denote −1.One has:Wq2= R Wq1 RThe base functions of the transforms can be computed astemplate <typename Type>void walsh_q1_basefunc(Type *f, ulong n, ulong k){ulong qk = (grs_negative_q(k) ? 1 : 0);k = gray_code(k);k = revbin(k, ld(n));for (ulong i=0; i<n; ++i){ulong x = i & k;x = parity(x);ulong qi = (grs_negative_q(i) ? 1 : 0);x ^= (qk ^ qi);f[i] = ( 0==x ? +1 : -1 );}}andtemplate <typename Type>void walsh_q2_basefunc(Type *f, ulong n, ulong k){ulong qk = (grs_negative_q(k) ? 1 : 0);k = revbin(k, ld(n));k = gray_code(k);(5.39)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:[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[* * ** * ** * ** ** *** * **** ** *** *************************** * * ******* ***** ** ** ** *** ********* * ** *** * *** * ***** * ** ** * * ******* * ***** * *** *** *** * ** * *** * **** ****** ******* * *** **** * *** *** * *** * *** * ** * * *** ** * ***** * ******* **** * ***** * ******** * ** * ** ***** *** ** ******* *** * ** * ** ******* * ****** *** *** *** **** *** *** * ** * ** * *** ***** * ********* *** * ****************** ** ** ** *** * ** ******* * ***** *** * ** **** ***** * *** **** ****** ** * ** * *** ******* * *** * *** * ***** *** * ** * ** ** ***** * ** ** *** *** * ******* * * ****** * ***97* ] (16)* ] (16)* ] (16)] (16)* ] (16)* ] (16)] (16)* ] (16)* ] (16)* ] (16)* ] (16)] (16)] (16)] (16)* ] (16)] (16)] (15)] (15)] (15)* ] (15)] (15)] (15)* ] (15)] (15)* ] (15)* ] (15)* ] (15)] (15)] (15)] (15)* ] (15)] (15)Figure 5.7: Basis functions for a self-inverse Walsh transform (second form) that has sequencies n/2 andn/2 − 1 only.
Asterisks denote the value +1, blank entries denote −1.}for (ulong i=0; i<n; ++i){ulong x = i & k;x = parity(x);ulong qi = (grs_negative_q(i) ? 1 : 0);x ^= (qk ^ qi);f[i] = ( 0==x ? +1 : -1 );}[FXT: walsh q1 basefunc and walsh q2 basefunc in walsh/walshbasefunc.h]5.7The slant transformThe slant transform can be implemented using a Walsh Transform and just a little pre/post-processing:void slant(double *f, ulong ldn)// slant transform{walsh_wak(f, ldn);ulong n = 1UL<<ldn;for (ulong ldm=0; ldm<ldn-1; ++ldm){ulong m = 1UL<<ldm; // m = 1, 2, 4, 8, ..., n/4double 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));CHAPTER 5. THE WALSH TRANSFORM AND ITS RELATIVES}}98for (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] = a * f1 - b * f2;f[t2] = b * f1 + a * f2;}The ldm-loop executes ldn−1 times, the inner loop is executed is n/2 − 1 times. That is, apart fromthe Walsh transform only an amount of work linear with the array size has to be done.