Fast Walsh-Hadamard Transform——快速沃爾什變換(二)


上次的博客有點模糊的說...我把思路和算法實現說一說吧...

 

思路

關於快速沃爾什變換,為了方便起見,我們采用線性變換(非線性變換不會搞)。

那么,就會有一個變化前各數值在變換后各處的系數,即前一篇博文中的$f(i,j)$,表示線性變換中第$i$項到第$j$項的系數。

$$DWT(A)_i = \sum_{j=0}^{n-1} A_j * f(i,j)$$

那么,我們既然要求$\oplus$卷積在變換后等價於乘積,就有

$$DWT(A)_i * DWT(B)_i = DWT(C)_i$$

其中$C$是$A,B$的$\oplus$卷積。

那么,將線性變換式及卷積定義代入,並令對應項系數相等,就可得到結論:

$$f(i, j) * f(i, k) = f(i, j \oplus k)$$

那么,通過構造合適的$f(i,j)$,就可得出變換。

$f$函數的構造

我們以異或為例(因為這時最麻煩的)。

既然是位運算,那么我們就從位運算入手。

考慮位運算,我們只需要考慮一位的情況,對於多位運算需要將每位的結果相乘(注意是相乘,這很重要!)。

通過定義,我們得到

$$f(i, 0) * f(i, j) = f(i, 0 \oplus j) = f(i, j)$$

$$f(i, 1) * f(i, 1) = f(i, 0)$$

又由於我們不能使所有系數都相同,那么只能取

$$\begin{cases}
f(i, 0) &= 1\\
f(0, j) &= 1\\
f(1, 1) &= -1
\end{cases}$$

那么此時

$$DWT(a_0, a_1) = (a_0 + a_1, a_0 - a_1)$$

可以驗證其滿足規則。

算法實現:

要實現$FWT$,我們仿照快速傅里葉變換,首先將$A$分成兩半$A_0, A_1$(按二進制最高位,即直接取前一半和后一半)遞歸調用$FWT$,然后,

合並時,根據下面的式子($i_0,i_1$分別表示$i$的最高位和剩余位):

$$\begin{aligned}
DWT(A)_i &= \sum_{j=0}^{n-1} f(i, j)A_j\\
&= \sum_{j=0}^{n/2-1}f(i,j)A_j + \sum_{j=n/2}^{n-1}f(i,j)A_j\\
&= \sum_{j=0}^{n/2-1}f(i_0,j_0)f(i_1,j_1)A_j + \sum_{j=n/2}^{n-1}f(i_0,j_0)f(i_1,j_1)A_j\\
&= \sum_{j=0}^{n/2-1}f(i_0,0)f(i_1,j_1)A_j + \sum_{j=n/2}^{n-1}f(i_0,1)f(i_1,j_1)A_j\\
&= f(i_0, 0)\sum_{j=0}^{n/2-1}f(i_1,j_1)A_j + f(i_1)\sum_{j=n/2}^{n-1}f(i_0,1)f(i_1,j_1)A_j\\
&= f(i_0, 0)*DWT(A_0)_{i_1} + f(i_0, 1)*DWT(A_1)_{i_1}\end{aligned}$$

其中由第2行到第3行是因為我們對多位的$f$直接定義為各位的$f$相乘,即

$$f(i,j)=f(i_0,j_0)*f(i_1,j_1)$$

那么,算法實現就簡單多了。

FWT(A, len)
1 divide A into (A0, A1) 
2 FWT(A0, len/2)
3 FWT(A1, len/2)
4 for i=0 to len/2
5   A[i]=f00*A0[i]+f01*A1[i]
6   A[i+len/2]=f10*A0[i]+f11*A1[i]

關於逆變換嘛。。。你只需要告訴自己:我只是要處理剛被$FWT$變換過的數組,我只需要把它倒過來。

那么

IFWT(A, len)
1 divide A into (A0, A1) 
for i=0 to len/2
3   solve the equation set, where the unknown numbers are A0[i] and A1[i]:
4     f00*A0[i]+f01*A1[i]=A[i]
5     f10*A0[i]+f11*A1[i]=A[i+len/2]
6 IFWT(A0, len/2)
7 IFWT(A1, len/2)

(當然解二元一次方程組只需要手算出來就好啦)

把三種運算(即與,或,異或)的f給出來:

$$\begin{array}{c|cccc}
Ops & f_{00} & f_{01} & f_{10} & f_{11} \\
\hline
And & 1 & 1 & 0 & 1 \\
Or & 1 & 0 & 1 & 1 \\
Xor & 1 & 1 & 1 & -1 
\end{array}
$$

 

附三種運算變換代碼:

 1 void FWT_And(int *A, int len) {
 2   if (len == 1) return;
 3   int len2 = len >> 1;
 4   FWT_And(A, len2);
 5   FWT_And(A + len2, len2);
 6   for (int i = 0; i < len2; ++i)
 7     A[i] += A[i + len2];
 8 }
 9 void IFWT_And(int *A, int len) {
10   if (len == 1) return;
11   int len2 = len >> 1;
12   for (int i = 0; i < len2; ++i)
13     A[i] -= A[i + len2];
14   IFWT_And(A, len2);
15   IFWT_And(A + len2, len2);
16 }
17 
18 void FWT_Or(int *A, int len) {
19   if (len == 1) return;
20   int len2 = len >> 1;
21   FWT_Or(A, len2);
22   FWT_Or(A + len2, len2);
23   for (int i = 0; i < len2; ++i)
24     A[i + len2] += A[i];
25 }
26 void IFWT_Or(int *A, int len) {
27   if (len == 1) return;
28   int len2 = len >> 1;
29   for (int i = 0; i < len2; ++i)
30     A[i + len2] -= A[i];
31   IFWT_Or(A, len2);
32   IFWT_Or(A + len2, len2);
33 }
34 
35 void FWT_Xor(int *A, int len) {
36   if (len == 1) return;
37   int len2 = len >> 1;
38   FWT_Xor(A, len2);
39   FWT_Xor(A + len2, len2);
40   for (int i = 0; i < len2; ++i) {
41     int x = A[i], y = A[i + len2];
42     A[i] = x + y;
43     A[i + len2] = x - y;
44   }
45 }
46 void IFWT_Xor(int *A, int len) {
47   if (len == 1) return;
48   int len2 = len >> 1;
49   for (int i = 0; i < len2; ++i) {
50     int x = A[i], y = A[i + len2];
51     A[i] = (x + y) >> 1;
52     A[i + len2] = (x - y) >> 1;
53   }
54   IFWT_Xor(A, len2);
55   IFWT_Xor(A + len2, len2);
56 }

 


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM