Bits reversal for CMSIS-DSP FFT

genBitsReversal.py

import math
import argparse
from sympy.combinatorics import Permutation

def bits_for_value(value):
  return int(math.log2(value))

def decompose(N, R):
  logN2 = bits_for_value(N)
  logR2 = []

  while (N >= R):
    logR2.append(bits_for_value(R))
    N = N // R

  if (N > 1):
    logR2.append(bits_for_value(N))

  return (logN2, logR2)

def reverse_bits(x, n, bits_list):
  result = 0
  for bits in bits_list:
    mask = (0xffffffff >> (32 - bits))
    result = (result << bits) | (x & mask)
    x = x >> bits
  return result

def create_transpositions(N, R):
  (logN2, logR2) = decompose(N, R)

  indexes = []
  for n in range(N):
    indexes.append(reverse_bits(n, logN2, logR2))

  # Create transpositions table
  tps = []
  for c in Permutation.from_sequence(indexes).cyclic_form:
    for i in range(len(c) - 1):
      tps.append([c[i] * 8, c[-1] * 8])

  return tps

def transpositions_stringify(N, R, tps):
  MAX_LINE_LEN   = 79
  MAX_FFT_IN_U16 = 8192

  index_type = 'uint16_t' if N <= MAX_FFT_IN_U16 else 'uint32_t'
  tps_elements_count = len(tps) * 2

  out = '#define ARMBITREVINDEXTABLE_{}_TABLE_LENGTH {}\n'.format(N, tps_elements_count)
  out += 'const {} armBitRevIndexTable{}[ARMBITREVINDEXTABLE_{}_TABLE_LENGTH] = {{\n'.format(index_type, N, N)

  line = ''
  for tp in tps:
    entry = '{},{}'.format(tp[0], tp[1])

    # Append to line
    exp_line_len = len(line) + len(entry) + len(', ,')

    if (line == ''):
      line = '   ' + entry
    elif (exp_line_len >= MAX_LINE_LEN):
      out += line + ',\n'
      line = '   ' + entry
    else:
      line += ', ' + entry

  out += line + '\n};'
  return out

parser = argparse.ArgumentParser(description='Generate bits reversal tables',
                                 formatter_class=argparse.ArgumentDefaultsHelpFormatter)

parser.add_argument('filename', metavar='out.c', nargs='?', help='output file name')
parser.add_argument('--size',  type=int, default=8192, help='size')
parser.add_argument('--radix', type=int, default=8, choices=[2, 8],
                               help='radix | use 2 for Radix 4 and 4x2 | use 8 for Radix 8, 8x4, 8x2')

args = parser.parse_args()

tps = create_transpositions(args.size, args.radix)
out = transpositions_stringify(args.size, args.radix, tps)

if (args.filename == None):
  print(out)
else:
  f = open(args.filename, 'w')
  f.write(out)
  f.close()

Requirements

pip3 install sympy

Tested using sympy 1.5.1

Usage

python3 genBitsReversal.py --size 8192 --radix 8 out.c

size must be the power of two

• use radix 2 for CMSIS Radix 4 and 4x2
• use radix 8 for CMSIS Radix 8, 8x4 and 8x2

if output file name is not specified, the script writes to standard output

The output doesn't match any table that is generated by ARM in the library.

For Example

fftGenLkup.py --size=16 --radix=8 
#define ARMBITREVINDEXTABLE_16_TABLE_LENGTH 20
const uint16_t armBitRevIndexTable16[ARMBITREVINDEXTABLE_16_TABLE_LENGTH] = {
   8,64, 16,64, 32,64, 24,72, 48,72, 96,72, 40,80, 56,88, 112,88, 104,88
};

https://raw.githubusercontent.com/ARM-software/CMSIS-DSP/2f24a05dee811794b5f32590c21544e811f4f719/Source/CommonTables/arm_common_tables.c

const uint16_t armBitRevIndexTable16[ARMBITREVINDEXTABLE_16_TABLE_LENGTH] ARM_DSP_TABLE_ATTRIBUTE =
{
   /* 8x2, size 20 */
   8,64, 24,72, 16,64, 40,80, 32,64, 56,88, 48,72, 88,104, 72,96, 104,112
};

I just have a very surface level understanding of this - shouldn't the output be the same?