testfft.c

来自「realview22.rar」· C语言 代码 · 共 527 行

/*
 * $Copyright: 
 * ----------------------------------------------------------------
 * This confidential and proprietary software may be used only as
 * authorised by a licensing agreement from ARM Limited
 *   (C) COPYRIGHT 1999,2000,2002 ARM Limited
 *       ALL RIGHTS RESERVED
 * The entire notice above must be reproduced on all authorised
 * copies and copies may only be made to the extent permitted
 * by a licensing agreement from ARM Limited.
 * ----------------------------------------------------------------
 * File:     testfft.c,v
 * Revision: 1.45
 * ----------------------------------------------------------------
 * $
 *
 * Test of ARM assembler FFT implementations
 *
 */

#include <stdio.h>
#include <math.h>
#include <time.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include "testfft.h"

#define VERBOSECOMP 0
#define TIMING 1
#define FAIL_THRESHOLD 24

#define MAXTESTFNS 32
static tFFTFunction testFunctions[MAXTESTFNS];
static int nTestFunctions;

#define NEWTEST(na, qs, inS, outS, fl) { \
  extern fFFT na; \
  test->fn = &na; \
  test->name = #na; \
  test->qshift = qs; \
  test->inSize = inS; \
  test->outSize = outS; \
  test->flags = fl ; \
  test++; nTestFunctions++; \
}

static void InitialiseTests(void) {
  tFFTFunction *test;
  
  nTestFunctions=0;
  test = &testFunctions[0];
  
  NEWTEST(FFT_4EFQ14R, 14, 2, 2, FFT_forward | FFT_evenlog | FFT_postshift | FFT_real);
  NEWTEST(FFT_4EIQ30R, 30, 4, 4, FFT_inverse | FFT_evenlog | FFT_postshift | FFT_real);
  NEWTEST(FFT_4OIQ14R, 14, 2, 2, FFT_inverse | FFT_oddlog  | FFT_postshift | FFT_real);
  NEWTEST(FFT_4OFQ30R, 28, 4, 4, FFT_forward | FFT_oddlog  | FFT_postshift | FFT_real);
  
  NEWTEST(FFT_4EFQ14, 14, 2, 2, FFT_forward | FFT_evenlog | FFT_postshift);
  NEWTEST(FFT_4EIQ30, 29, 4, 4, FFT_inverse | FFT_evenlog | FFT_postshift);
  NEWTEST(FFT_4OIQ14, 14, 2, 2, FFT_inverse | FFT_oddlog  | FFT_postshift);
  NEWTEST(FFT_4OFQ30, 28, 4, 4, FFT_forward | FFT_oddlog  | FFT_postshift);
  
//  NEWTEST(FFT_4EFQ30, 29, 4, 4, FFT_evenlog | FFT_postshift);
//  NEWTEST(FFT_4EIQ14, 14, 2, 2, FFT_inverse | FFT_evenlog | FFT_postshift);
//  NEWTEST(FFT_4OIQ14, 14, 2, 2, FFT_inverse | FFT_oddlog | FFT_postshift);
}

/* global value */
static real pi;

static int log2(int N) {
  int k=0;
  
  if (N>>16) { k+=16; N>>=16; }
  if (N>>8) { k+=8; N>>=8; }
  if (N>>4) { k+=4; N>>=4; }
  if (N>>2) { k+=2; N>>=2; }
  if (N>>1) { k+=1; N>>=1; }

  return k;
}

/* Perform bit reversal on an array of objects of given size
 * Currently only supports N a power of 2
 */

static void BitReverse(char *out, char *in, unsigned N, size_t size) {
  unsigned i, ibrev;
    
  for (i=0, ibrev=0; i<N; i++) {
    unsigned b;
      
    /* move next element */
    memcpy(out+ibrev*size, in+i*size, size);
    
    /* bit reverse increment of ibrev */
    for (b = N>>1; ibrev & b; b>>=1) ibrev ^= b;
    ibrev ^= b;
    
  }
}

static void FFT_stage_fp(
  comp *data,
  comp *table,
  unsigned nBlocks,
  unsigned entriesPerBlock
)
{
  int i,j;

  for (i=0; i<entriesPerBlock/2; i++) {
    real c,s;
    comp *dl,*dr;

    s = table->y;
    c = table->x;
    table += nBlocks;
    dl = data+i;
    dr = dl+(entriesPerBlock/2);

    for (j=nBlocks; j!=0; j--) {
      real x,y,z;

      x = dr->x;
      y = dr->y;
      z = x * c - y * s; // real part of complex product
      y = x * s + y * c; // imag part of complex product
      x = dl->x;
      dl->x = x+z;
      dr->x = x-z;
      x = dl->y;
      dl->y = x+y;
      dr->y = x-y;

      dl += entriesPerBlock;
      dr += entriesPerBlock;
    }
  }
}

/* Create cos/sin table for FFT
 * This can be optimised a lot since N is a power of 2
 *
 * Given and t on input it creates a table of
 * 1 w w^2 w^3 ...... w^(N-1)
 *
 * where w = exp(it) = cos(t) + i*sin(t)
 */

static void FFT_table_fp(comp *table, real t, unsigned N) {
  unsigned i,p;
  real s, c;
  
  printf("Creating cos/sin table\n");
  c = cos(t);
  s = sin(t);
  
  table[0].x = (real)1;
  table[0].y = (real)0;

  for (p=1; p<N; p<<=1) {
    
    /* extend table to be 2*p long */
    for (i=0; i<p; i++) {
      real x,y;
      
      x = table[i].x;
      y = table[i].y;
      table[p+i].x = x * c - y * s;
      table[p+i].y = x * s + y * c;
    }
    
    /* square w */
    s = (real)2*c*s;
    c = (real)2*c*c - (real)1;
  }

  printf("Finished\n");
}

int FFT_fp(comp *input, comp *output, int logN, int direction) {
  real w;
  unsigned nBlocks;
  unsigned entriesPerBlock;
  comp *table;
  unsigned N;
  
  N = 1<<logN;
  
  /* do bit reversal */
  BitReverse((char*)output, (char*)input, N, sizeof(comp));

  /* do FFT stages */
  if (direction) {
    w = -pi;
  } else {
    w = pi;
  }
  nBlocks = N/2;
  entriesPerBlock = 2;

  /* create table */
  table = (comp *)malloc(N/2*sizeof(comp));
  if (!table) {
    printf("Couldn't claim table\n");
    return -1;
  }
  FFT_table_fp(table, w/nBlocks, nBlocks);

  do {
    printf("nblock = %d\n", nBlocks);
    FFT_stage_fp(output, table, nBlocks, entriesPerBlock);
    w = w/2.0;
    nBlocks >>= 1;
    entriesPerBlock  <<= 1;
  } while (nBlocks!=0);

  free(table);
  return 0;
}


static int compare_output(
  tFFTFunction *fft,
  void *out,
  real *ref,
  int N,
  int qshift,
  int reconstruct
)
{
  int i, v, d, r;
  real Q = (real)((unsigned int)(1<<qshift));
  unsigned int err=0;
  int size = fft->outSize;
  int rflag = fft->flags & FFT_real;
  int conjugate = 0;
  
  /* do we need to conjugate the output? */
  if (fft->flags & FFT_inverse) {
    conjugate ^= 1;
  }
  if (reconstruct) {
    conjugate ^= 1;
  }

  /* print out reference results first */
  printf("  Ref: ");
  for (i=0; i<10*2 && i<N; i++) {
    v = (int)(ref[i]*Q);
    if (i&conjugate) v=-v;
    if (i==1 && rflag) v=(int)(ref[N]*Q);
    if (i&1) printf("%+di", v);
    else     printf(" %d", v);
  }
  printf("\n");

  /* now compare actual results */  
  printf("  Act: ");
  for (i=0; i<N; i++) {
    if (size==4) {
      v = (int) (*((int *)out+i));
    } else {
      v = (int) (*((short *)out+i));
    }
    r = (int)(ref[i]*Q);
    
    /* check if we need to conjugate the output */
    if (i & conjugate) {
      r = -r;
    }
    
    /* for a real FFT output we place the half way result at offset 1 */
    if (i==1 && rflag) {
      r = (int)(ref[N]*Q);
    }
    
#if VERBOSECOMP
    printf(" ref[%d]=%d out[%d]=%d\n", i, r, i, v);
#else
    if (i<10*2) {
      if (i&1) printf("%+di", v);
      else     printf(" %d", v);
    }
#endif
    d = v - r;
    if (d<0) d=-d;
    if (d>err) {
      err=d;
#if VERBOSECOMP
      printf("New max error of %d: [%d] %d %d\n", d, i, r, v);
#endif
    }
  }
  printf("\n  Rounding error = %d (%lf)\n\n", err, (real)err/Q);
  if (err>=FAIL_THRESHOLD) return 1;
  return 0;
}

  
static void comp_input(void *out, comp *in, tFFTFunction *transform, int N) {
  real Q = (real)((unsigned int)(1<<transform->qshift));
  int inverse=0;
  int size = transform->inSize;
  int i;
  
  if (transform->flags & FFT_inverse) inverse=1;
  for (i=0; i<2*N; i++) {
    real v = ((real*)in)[i] * Q;
    
    /* conjugate for an inverse transform */
    if (i & inverse) v=-v;
    
    if (size==4) {
      *((int *)out+i) = (int)v;
    } else {
      *((short *)out+i) = (short)v;
    }
  }
}

static void real_input(void *out, comp *in, int size, int N, int qshift) {
  int i;
  real Q = (real)((unsigned int)(1<<qshift));
  
  for (i=0; i<N; i++) {
    real v = in[i].x * Q;
    
    if (size==4) {
      *((int *)out+i) = (int)v;
    } else {
      *((short *)out+i) = (short)v;
    }
  }
}

static void shift_input(void *out, void *in, int size, int N, int shift) {
  int i;
  
  for (i=0; i<N; i++) {
    if (size==4) {
      *((int *)out+i) = *((int *)in+i) >> shift;
    } else {
      *((short *)out+i) = (short)(*((short *)in+i) >> shift);
    }
  }
}

static int DoFFT(tFFTFunction *fft, void *src, void *dest, int NC, int qshift) {
  int result;
#if TIMING
  int start;
#endif
  
  /* test prescale feature if present */
  if (fft->flags & FFT_prescale) {
    shift_input(src, src, fft->inSize, 2*NC, 1);
  }
  
  /* bit reverse the input if the FFT requires this */
  if (fft->flags & FFT_reversed) {
    BitReverse((char*)dest, (char*)src, NC, 2*fft->inSize);
    src=dest;
  }
  
#if TIMING
  start = clock();
#endif
  result = fft->fn(src, dest, NC, 1);
#if TIMING
  start = clock()-start;
  printf("  Time taken = %d\n", start);
#endif

  if (result) {
    printf("FFT failed: %d - ", result);
    switch (result) {
      case 1: printf("Table too small\n"); break;
    }
    printf("\n");
    return 0xDEAD;
  }
  
  if (fft->flags & FFT_outinbuf) {
    shift_input(dest, src, fft->outSize, 2*NC, 0);
  }
  
  qshift -= log2(NC);
  if (fft->flags & FFT_real) qshift--;
  
  return qshift;
}

/* returns 1 if the FFT failed */

static int TestFFT(
  tFFTTest *FFTTest,            /* FFT test data */
  tFFTFunction *transform,      /* FFT transform function to test */
  tFFTFunction *recover         /* transform to recover the data if available */
)
{
  void *src, *destFor;
  void *destInv;
  int qshift;
  int N=FFTTest->N;
  int NC=N;                     /* number of complex points */
  int result;
  
  if (!transform) return 0;
  
  if (transform->flags & FFT_real) {
    /* FFT can only handle real data and produces
     * half the number of outputs
     */
    
    if (!(FFTTest->flags & FFT_real)) return 0;
    NC = N/2;
  }
  
  if (NC&0x55555555) {
    /* even parity NC */
    if (!(transform->flags & FFT_evenlog)) return 0;
  } else {
    /* odd parity NC */
    if (!(transform->flags & FFT_oddlog)) return 0;
  }
  
  src     = malloc(2*NC*transform->inSize);
  destFor = malloc(2*NC*transform->outSize);
  
  if (!src || !destFor) {
    printf("Couldn't claim FFT buffers\n");
    return 0;
  }
  
  printf("Performing %s on %d-point test data: %s\n",
    transform->name, N, FFTTest->name);
  if (NC==N) {
    comp_input(src, FFTTest->in, transform, N);
  } else {
    printf("Seeding real data\n");
    real_input(src, FFTTest->in, transform->inSize, N, transform->qshift);
  }
  qshift = DoFFT(transform, src, destFor, NC, transform->qshift);
  if (qshift==0xDEAD) return 1;
  result = compare_output(transform, destFor, (real *)FFTTest->out, 2*NC, qshift, 0);
  
  if (recover && 0) {
    
    printf("Recovering data using %s\n", recover->name);
    assert(transform->outSize == recover->inSize);
    destInv = malloc(2*N*recover->outSize);
    if (!destInv) {
      printf("Couldn't claim FFT inverse buffer\n");
      return 0;
    }
    
    qshift = DoFFT(recover, destFor, destInv, N, qshift);
    qshift += log2(N);
    compare_output(recover, destInv, (real *)FFTTest->in, 2*N, qshift, 1);
    free(destInv);
  }
  
  free(src);
  free(destFor);

  return result;
}

static int RunFFTTest(tFFTTest *FFTTest) {
  int i, failed=0;
  
  for (i=0; i<nTestFunctions; i++) {
    failed += TestFFT(FFTTest, &testFunctions[i], NULL);
  }
//  TestFFT(FFTTest, &TEST_4EFQ14, &TEST_4EIQ14);
//  TestFFT(FFTTest, &TEST_4OFQ14, &TEST_4OIQ14);
//  TestFFT(FFTTest, &TEST_4EIQ14, &TEST_4EIQ14);
//  TestFFT(FFTTest, &TEST_4OIQ14, &TEST_4OIQ14);
//  TestFFT(FFTTest, &TEST_4EFQ30, NULL);
//  TestFFT(FFTTest, &TEST_4OIQ14R, NULL);
//  TestFFT(FFTTest, &TEST_4OIQ30R, NULL);
//  TestFFT(FFTTest, &TEST_4EIQ14R, NULL);
  return failed;
}

static int TestFFTSize(const tFFTTests *FFTTests) {
  int N,i,failed=0;
  
  N = FFTTests->N;
  printf("Testing %d point FFT's\n", N);
  for (i=0; i<FFTTests->n; i++) {
    failed+=RunFFTTest(&FFTTests->FFTTests[i]);
  }
  return failed;
}

#define TestSize(n) { \
  extern const tFFTTests FFTTests_##n; \
  failed += TestFFTSize(&FFTTests_##n); \
}

int main() {
  int failed=0;
  
  InitialiseTests();
  TestSize(8);
  TestSize(16);
  TestSize(32);
  TestSize(64);
  TestSize(128);
  TestSize(256);
  TestSize(512);
  TestSize(1024);
  TestSize(2048);
  TestSize(4096);
  
  if (failed) {
    printf("There were %d failures\n", failed);
  } else {
    printf("Test passed\n");
  }
  printf("Finished\n");
}