benchmarks/vec-fft/fft-data-gen

   1 #!/usr/bin/env python
   2
   3 import sys
   4 import numpy
   5
   6 def array_print(name, data):
   7         if type_scalar == numpy.float16:
   8                 a = ['{:#04x}'.format(h.astype(int)) for h in data.view(dtype=numpy.uint16)]
   9         else:
  10                 a = [repr(x) for x in data]
  11         print 'fftval_t {}[{}] = {{ {} }};'.format(name, len(a), ', '.join(a))
  12
  13 def array_permute(data, radix = 2):
  14         logradix = int(numpy.log2(radix))
  15         term_mask = int(radix - 1)
  16 #       num_term = int(numpy.log2(len(data)) / logradix)
  17         for i in xrange(0, len(data)):
  18                 # Obtain permuted address
  19                 i_left = i
  20                 permuted = 0
  21                 cur_fft_size = radix
  22                 while cur_fft_size <= len(data):
  23                         permuted = (permuted << logradix) | (i_left & term_mask)
  24                         i_left >>= logradix
  25                         cur_fft_size <<= logradix
  26                 # Permute only once and when addresses are different
  27                 if i < permuted:
  28                         tmp = data[i]
  29                         data[i] = data[permuted]
  30                         data[permuted] = tmp
  31
  32 fft_size = 1024
  33 type_scalar = numpy.float16
  34 permute = True
  35
  36 if len(sys.argv) > 1:
  37         type_scalar = {
  38                 '16' : numpy.float16,
  39                 '32' : numpy.float32,
  40                 '64' : numpy.float64 }.get(sys.argv[1].strip(), None)
  41         if type_scalar == None:
  42                 sys.exit('Invalid datatype')
  43
  44 if len(sys.argv) > 2:
  45         fft_size = sys.argv[2].strip()
  46         if not fft_size.isdigit():
  47                 sys.exit('Invalid FFT size')
  48         fft_size = int(fft_size)
  49
  50 numpy.random.seed(seed=0)
  51
  52 print '#include "fft_const.h"'
  53
  54 # Generate input data
  55 input_real = numpy.random.uniform(size=fft_size).astype(type_scalar)
  56 input_imag = numpy.random.uniform(size=fft_size).astype(type_scalar)
  57
  58 # Compute reference FFT output
  59 type_cmplx = numpy.complex128 if (type_scalar == numpy.float64) else numpy.complex64
  60 input_cmplx = input_real.astype(type_cmplx)
  61 input_cmplx.imag = input_imag
  62 output_cmplx = numpy.fft.fft(input_cmplx)
  63
  64 if permute:
  65         array_permute(input_real)
  66         array_permute(input_imag)
  67
  68 array_print('input_data_real', input_real)
  69 array_print('input_data_imag', input_imag)
  70
  71 output_real = output_cmplx.real.astype(type_scalar)
  72 output_imag = output_cmplx.imag.astype(type_scalar)
  73
  74 array_print('output_data_real', output_real)
  75 array_print('output_data_imag', output_imag)
  76
  77 # Generate twiddle factors (TFs)
  78 # Negate sine since twiddle factor angles are generated *clockwise* from 0
  79 rad = [(2.0 * numpy.pi * (float(i) / fft_size)) for i in xrange(fft_size)]
  80 tf_real = numpy.cos(rad).astype(type_scalar)
  81 tf_imag = -numpy.sin(rad).astype(type_scalar)
  82
  83 array_print('tf_real', tf_real)
  84 array_print('tf_imag', tf_imag)