FFT: ビット列反転インデックス

FFT をシングル バッファ（入力も出力も同じ配列）で処理するとき欲しくなるのは、ビット列反転インデックス。

インデックスから、メモリ アドレスを計算するときは、index のビット幅で反転して欲しいので
　　address = base + (bit_reverse(index, width) << N)
だけど、これをプロセッサ側で実装してもらったら
　　address = base + (bit_reverse(index) >> K)
になりそう。

使った事はないのでよくわからないが、ビット列反転アドレッシング モードを持っている DSP があるらしい。

FFT 以外では、RFC1951 (ZLIB 圧縮) のハフマン符号のところで欲しくなります。

---

32ビット幅のビット列反転

まずは 32ビット幅のビット列反転を実装してみる。

brev_std.cpp

// ビット並び反転の１ステップ
inline unsigned int bit_reverse_sw(unsigned int n, unsigned int m, int s)
{
  return ((n >> s) & m) | ((n & m) << s);
}

// ビット並びを反転させる
unsigned int bit_reverse(unsigned int n)
{
  n = bit_reverse_sw(n, 0x0000ffff, 16);
  n = bit_reverse_sw(n, 0x00ff00ff, 8);
  n = bit_reverse_sw(n, 0x0f0f0f0f, 4);
  n = bit_reverse_sw(n, 0x33333333, 2);
  n = bit_reverse_sw(n, 0x55555555, 1);
  return n;
}

x86 でコンパイルしたアセンブリ出力

$ g++ -S -O3 brev_std.cpp -o brev_x86.s
Apple clang version 11.0.0 (clang-1100.0.33.17)
Target: x86_64-apple-darwin19.3.0

brev_x86.s

        .section        __TEXT,__text,regular,pure_instructions
        .build_version macos, 10, 15    sdk_version 10, 15
        .globl  __Z11bit_reversej       ## -- Begin function _Z11bit_reversej
        .p2align        4, 0x90
__Z11bit_reversej:                      ## @_Z11bit_reversej
        .cfi_startproc
## %bb.0:
        pushq   %rbp
        .cfi_def_cfa_offset 16
        .cfi_offset %rbp, -16
        movq    %rsp, %rbp
        .cfi_def_cfa_register %rbp
                                        ## kill: def $edi killed $edi def $rdi
        bswapl  %edi
        movl    %edi, %eax
        shrl    $4, %eax
        andl    $252645135, %eax        ## imm = 0xF0F0F0F
        shll    $4, %edi
        andl    $-252645136, %edi       ## imm = 0xF0F0F0F0
        orl     %eax, %edi
        movl    %edi, %eax
        shrl    $2, %eax
        andl    $858993459, %eax        ## imm = 0x33333333
        andl    $858993459, %edi        ## imm = 0x33333333
        leal    (%rax,%rdi,4), %eax
        movl    %eax, %ecx
        shrl    %ecx
        andl    $1431655765, %ecx       ## imm = 0x55555555
        andl    $1431655765, %eax       ## imm = 0x55555555
        leal    (%rcx,%rax,2), %eax
        popq    %rbp
        retq
        .cfi_endproc
                                        ## -- End function

.subsections_via_symbols

では、良いコードがでているようです。

ARM でコンパイルしたアセンブリ出力

$ g++ -S -O3 brev_std.cpp -o brev_arm.s
gcc version 8.3.0 (Raspbian 8.3.0-6+rpi1)
Target: arm-linux-gnueabihf

brev_arm.s

        .type   _Z11bit_reversej, %function
_Z11bit_reversej:
        .fnstart
.LFB1:
        @ args = 0, pretend = 0, frame = 0
        @ frame_needed = 0, uses_anonymous_args = 0
        @ link register save eliminated.
        rev     r0, r0
        ldr     r2, .L3
        ldr     r3, .L3+4
        and     r3, r3, r0, lsr #4
        and     r0, r2, r0, lsl #4
        orr     r0, r3, r0
        ldr     r2, .L3+8
        ldr     r3, .L3+12
        and     r3, r3, r0, lsr #2
        and     r0, r2, r0, lsl #2
        orr     r3, r3, r0
        ldr     r2, .L3+16
        ldr     r0, .L3+20
        and     r0, r0, r3, lsl #1
        and     r3, r2, r3, lsr #1
        orr     r0, r0, r3
        bx      lr
.L4:
        .align  2
.L3:
        .word   -252645136
        .word   252645135
        .word   -858993460
        .word   858993459
        .word   1431655765
        .word   -1431655766
        .cantunwind
        .fnend

は、ちょっと気になります。

インライン アセンブラで書き直す

brev_arm.cpp

// ビット並びを反転させる
unsigned int bit_reverse(unsigned int n)
{
  unsigned int h, l, m;

  asm("rev\t%0, %1" : "=r" (n) : "r" (n));

  m = 0x0f0f0f0f;
  asm("and\t%0, %2, %1, lsr #4": "=r" (l) : "r" (n), "r" (m));
  asm("and\t%0, %1, %2"        : "=r" (h) : "r" (n), "r" (m));
  asm("orr\t%0, %1, %2, lsl #4": "=r" (n) : "r" (l), "r" (h));

  m = 0x33333333;
  asm("and\t%0, %2, %1, lsr #2": "=r" (l) : "r" (n), "r" (m));
  asm("and\t%0, %1, %2"        : "=r" (h) : "r" (n), "r" (m));
  asm("orr\t%0, %1, %2, lsl #2": "=r" (n) : "r" (l), "r" (h));

  m = 0x55555555;
  asm("and\t%0, %2, %1, lsr #1": "=r" (l) : "r" (n), "r" (m));
  asm("and\t%0, %1, %2"        : "=r" (h) : "r" (n), "r" (m));
  asm("orr\t%0, %1, %2, lsl #1": "=r" (n) : "r" (l), "r" (h));

  return n;
}

アセンブリ出力は（余計な部分は削除してあります）

brev_arm.s

_Z11bit_reversej:
        .fnstart
.LFB0:
        @ args = 0, pretend = 0, frame = 0
        @ frame_needed = 0, uses_anonymous_args = 0
        @ link register save eliminated.
        ldr     r2, .L3
        ldr     r3, .L3+4
        ldr     ip, .L3+8
        rev     r0, r0
        and     r1, r2, r0, lsr #4
        and     r0, r0, r2
        orr     r1, r1, r0, lsl #4
        and     r2, r3, r1, lsr #2
        and     r1, r1, r3
        orr     r3, r2, r1, lsl #2
        and     r0, ip, r3, lsr #1
        and     r3, r3, ip
        orr     r0, r0, r3, lsl #1
        bx      lr
.L4:
        .align  2
.L3:
        .word   252645135
        .word   858993459
        .word   1431655765
        .cantunwind
        .fnend

ldr 命令と .word がともに３つ減りました。 ldr の位置が気になりますが…

---

配列の大きさからビット数を得る

整数版 log2を使います。

配列の大きさ N は $2^n$ でないと困りますが

((N > 0) && (N == (N & -N))) /* 真で 2^n */

で判定できます。

---

ビット列反転インデックス

動作テストのためのサンプル プログラム

sample.cpp

#include "log2u.h"
#include "brev_std.cpp"

#include <iostream>
#include <iomanip>

using namespace std;

void dump_bin(unsigned int n, int l)
{
  for (int i = 0; i < l; i++)
    {
      int b = l - i - 1;
      cout << ((n & (1 << b)) ? "1" : "0");
    }
}

int main(int argc, char **argv)
{
  if (argc != 2)
    {
      cout << "Usage: " << argv[0] << " N" << endl;
      return 1;
    }

  char *arg = argv[1];
  int base = 10;

  if (arg[0] == '0')
    {
      switch (arg[1])
        {
        case 'X': case 'x': base = 16; arg += 2; break;
        case 'B': case 'b': base =  2; arg += 2; break;
        }
    }

  unsigned int N = strtoul(arg, nullptr, base);

  if (!((N > 0) && (N == (N & -N))))
    {
      cout << "error: " << N << " != 2^n" << endl;
      return 2;
    }

  int W = log2u(N);
  int K = 32 - W;
  int hw = (W + 3) >> 2;

  for (unsigned int i = 0; i < N; i++)
    {
      unsigned int ri = bit_reverse(i) >> K; // ビット列反転インデックス

      cout << "0x" << hex << setw(hw) << setfill('0') << i << " --> ";
      dump_bin(ri, W);
      cout  << " : 0x" << hex << setw(hw) << setfill('0') << ri << endl;
    }

  return 0;
}

で動作確認する。

$ ./sample 8
0x0 --> 000 : 0x0
0x1 --> 100 : 0x4
0x2 --> 010 : 0x2
0x3 --> 110 : 0x6
0x4 --> 001 : 0x1
0x5 --> 101 : 0x5
0x6 --> 011 : 0x3
0x7 --> 111 : 0x7

単純増加のループならば bit_reverse は不要になるので

-  for (unsigned int i = 0; i < N; i++)
-    {
-      unsigned int ri = bit_reverse(i) >> K; // ビット列反転インデックス
-
-      cout << "0x" << hex << setw(hw) << setfill('0') << i << " --> ";
-      dump_bin(ri, W);
-      cout  << " : 0x" << hex << setw(hw) << setfill('0') << ri << endl;
-    }
-
+  unsigned int bm = N - 1;
+  unsigned int b0 = N >> 1;
+  unsigned int ri = 0; // ビット列反転インデックス
+  for (unsigned int i = 0; i < N; i++)
+    {
+      cout << "0x" << hex << setw(hw) << setfill('0') << i << " --> ";
+      dump_bin(ri, W);
+      cout  << " : 0x" << hex << setw(hw) << setfill('0') << ri << endl;
+
+      // ビット列反転インデックスを更新
+#if 1
+      unsigned int bn = b0;
+      do
+        {
+          ri ^= bn;
+          if ((ri & bn) != 0)
+            break;
+          bn >>= 1;
+        }
+      while (bn != 0);
+#else
+      unsigned int ii = i + 1;
+      int bs = log2u(ii & -ii);
+      ri = (ri & (bm >> bs)) ^ (b0 >> bs);
+#endif
+    }

とできます。

変更部分のアセンブリ出力：GCC (x86)

;#if 1 の場合、実行される命令数は不定

        movl    -52(%rbp), %eax         ## 4-byte Reload
        .p2align        4, 0x90
LBB2_28:                                ##   Parent Loop BB2_25 Depth=1
                                        ## =>  This Inner Loop Header: Depth=2
        xorl    %eax, %r12d
        testl   %eax, %r12d
        jne     LBB2_30
## %bb.29:                              ##   in Loop: Header=BB2_28 Depth=2
        shrl    %eax
        jne     LBB2_28
LBB2_30:                                ##   in Loop: Header=BB2_25 Depth=1

;#else の場合、実行される命令数は固定

        incl    %r14d
        movl    %r14d, %eax
        movl    -64(%rbp), %edx         ## 4-byte Reload
        andl    %edx, %eax
        movl    $-1, %ecx
        ## InlineAsm Start
        bsrl    %eax, %ecx
        ## InlineAsm End
        movl    -56(%rbp), %eax         ## 4-byte Reload
        shrl    %cl, %eax
        andl    %r12d, %eax
        movl    -52(%rbp), %r12d        ## 4-byte Reload
                                        ## kill: def $cl killed $cl killed $ecx
        shrl    %cl, %r12d
        xorl    %eax, %r12d

;#endif