gentests.c

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

/*
 * $Copyright: 
 * ----------------------------------------------------------------
 * This confidential and proprietary software may be used only as
 * authorised by a licensing agreement from ARM Limited
 *   (C) COPYRIGHT 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:     gentests.c,v
 * Revision: 1.9
 * ----------------------------------------------------------------
 * $
 *
 * Create the test vectors required for the FFT
 * Ideally compile this file natively (eg using 'cl' or 'gcc')
 *  as it takes a long time to generate the tables
 *
 */

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

#define MAXTESTS 50
tFFTTest LFFTTests[MAXTESTS];
static int nLFFTTests=0;

void generate_test(FILE *fp, int N, int type, int flags);
void parabola(tFFTTest *FFTTest, comp s, comp b, comp c);
void sinewave(tFFTTest *FFTTest, comp s, int k1, int k2);
void normalise(tFFTTest *FFTTest);
void output(FILE *fp, comp *in, int N, int nTest, char *name);
void output_tests(FILE *fp, int N);

int main(int argc, char **argv) {
  int N;
  char path[256];
  FILE *fp;
  
  /* set defaults */
  N=1024;

  /* read command line arguments */
  if (--argc>0) {
    N = atoi(*(++argv)); 
  } else {
    printf("Please enter the FFT size to test:\n");
    scanf("%d", &N);
    if (N==0) N=2048;
  }
  printf("FFT size to test = %d\n", N);
  
  /* formulate filename */
  sprintf(path, "t%06d.c", N);
        
  printf("Output file = %s\n", path);
  fp = fopen(path, "wt");
  if (fp == NULL) {
    printf ("Couldn't open file %s for output.\n", path);
    exit (1);
  }
        
  fprintf(fp, "/*\n");
  fprintf(fp, " * Copyright (C) ARM Limited 1998-2002. All rights reserved.\n");
  fprintf(fp, " *\n");
  fprintf(fp, " * t%06d.c\n", N);
  fprintf(fp, " *\n");
  fprintf(fp, " * Tests for FFT's of size N=%d\n", N);
  fprintf(fp, " */\n");
  fprintf(fp, "\n");
  fprintf(fp, "#include \"testfft.h\"\n\n");
  
  /* generate real parabola and sinewave test patterns */
  generate_test(fp, N, 0, FFT_real);
  generate_test(fp, N, 1, FFT_real);
  
  /* generate complex parabola and sinewave test patterns */
  generate_test(fp, N, 0, 0);
  generate_test(fp, N, 1, 0);

  output_tests(fp, N);
  printf ("Tests created\n");
  return 0;
}

/* N = number of points in FFT to test (need not be a power of 2)
 * type = the type of shape to test
 *        0 => parabola input shape (accuracy test)
 *        1 => sinwave input shape (maximum abs value test)
 * flags = direction and type
 *        FFT_real
 *        FFT_inverse => test FFT in the inverse direction
 * The range of the test input is scaled to absolute value of 1
 * on each scalar.
 */

void generate_test(FILE *fp, int N, int type, int flags) {
  tFFTTest *FFTTest;
  int nTest;
  comp s;
  
  assert(nLFFTTests<MAXTESTS);
  nTest = nLFFTTests++;
  FFTTest = &LFFTTests[nTest];
  
  printf("Creating test N=%d type=%d flags=%d\n", N, type, flags);
  
  FFTTest->N = N;
  FFTTest->flags = flags;
  FFTTest->in  = (comp*)malloc(N*sizeof(comp));
  FFTTest->out = (comp*)malloc(N*sizeof(comp));
  
  assert(FFTTest->in);
  assert(FFTTest->out);
  
  s.x = 1.0;
  s.y = 0.3;
  if (flags & FFT_real) s.y=0;
  
  if (type==0) {
    /* parabola */
    comp b, c;
    
    FFTTest->name="parabola-complex";
    b.x = (real)N*(real)4/(real)3;
    b.y = b.x * 0.3;
    c.x = (real)N*(real)4/(real)9;
    c.y = c.x * 0.7;
    if (flags & FFT_real) {
      FFTTest->name="parabola-real";
      b.y=0;
      c.y=0;
    }
    parabola(FFTTest, s, b, c);
  }
  
  if (type==1) {
    /* sinewave */
    if (flags & FFT_real) {
      FFTTest->name="sinewave-real";
      sinewave(FFTTest, s, 3, 3);
    } else {
      FFTTest->name="sinewave-complex";
      sinewave(FFTTest, s, 3, N-3);
    }
  }
  
  normalise(FFTTest);
  output(fp, FFTTest->in, N, nTest, "in");
  output(fp, FFTTest->out, N, nTest, "out");
  free(FFTTest->in);
  free(FFTTest->out);
}

/* Take the FFT of the parabola:
 *
 * x[k] = s*(k*k + b*k + c)
 *
 * if x[k]=1 then X=FFT(x) has
 *  X[0] = N
 *  X[k] = 0 (k!=0)
 *
 * if x[k]=k then X=FFT(x) has
 *  X[0] = N*(N-1)/2
 *  X[k] = N/(w^k-1) (k!=0)
 *
 * if x[k]=k*k then X=FFT(x) has
 *  X[0] = N*(N-1)*(2N-1)/6
 *  X[k] = N/(w^k-1) * (N-2-2/(w^k-1)) (k!=0)
 *
 * if x[k]= s*(k*k + b*k + c) then X=FFT(x) has
 *  X[k] = (-s*N/2) 2/(1-w^k) * (N-2+b + 2/(1-w^k))
 *
 */

void parabola(tFFTTest *FFTTest, comp s, comp b, comp c) {
  real x,y,n1,n2,t;
  int k;
  int N=FFTTest->N;
  comp *in =FFTTest->in;
  comp *out=FFTTest->out;
  
  /* t = - 2*pi/N */
  t = -(real)8*atan(1.0)/(real)N;
  
  /* write the input parabola */
  for (k=0; k<N; k++) {
    x = k*k + b.x*k + c.x;
    y = 0   + b.y*k + c.y;
    in[k].x = x*s.x - y*s.y;
    in[k].y = x*s.y + y*s.x;
  }
  
  /* write the output transformed data
   * NB if w = cos(t) + i*sin(t)  then
   * 2/(1-w^k) = 1 + i*[sin(kt)/(1-cos(kt))]
   */
  
  n1 = (real)(N-1)/(real)2;
  n2 = (real)(2*N-1)/(real)3;
  x = ((n2 + b.x)*n1 + c.x)*(real)N;
  y = (( 0 + b.y)*n1 + c.y)*(real)N;
  out[0].x = x*s.x - y*s.y;
  out[0].y = x*s.y + y*s.x;
  
  for (k=1; k<N; k++) {
    real x1,y1,wi;
    real ang = (real)k*t;
    
    wi = sin(ang)/((real)1-cos(ang));
    x1 = (real)(N-1) + b.x;
    y1 = b.y + wi;
    x = x1 - y1*wi;
    y = y1 + x1*wi;
    wi = (-(real)(N/2));
    out[k].x = wi*(x*s.x - y*s.y);
    out[k].y = wi*(x*s.y + y*s.x);
  }
}

/* Take FFT of the sinewave:
 *
 * x[k] = s*(exp(i*k*k1*t) + exp(-i*k*k2*t))
 *
 * where t=2*pi/N is the angle to FFT by. Then,
 *
 * X[k] = s * N * ( (k==k1) + (k==N-k2) )
 */

void sinewave(tFFTTest *FFTTest, comp s, int k1, int k2) {
  real x,y,t;
  int k;
  int N=FFTTest->N;
  comp *in =FFTTest->in;
  comp *out=FFTTest->out;
  
  /* t = 2*pi/N */
  t = (real)8*atan(1.0)/(real)N;
  
  /* range k1 and k2 bins */
  while (k1<0) k1+=N;
  while (k1>=N) k1-=N;
  while (k2<=0) k2+=N;
  while (k2>N) k2-=N;
  
  /* write the input sinewave */
  for (k=0; k<N; k++) {
    real ang = k*t;
    
    x = cos(ang*k1) + cos(ang*k2);
    y = sin(ang*k1) - sin(ang*k2);
    in[k].x = x*s.x - y*s.y;
    in[k].y = x*s.y + y*s.x;
  }
  
  /* write the output transformed data */
  for (k=0; k<N; k++) {
    out[k].x = 0;
    out[k].y = 0;
  }
  out[k1].x += s.x*N;
  out[k1].y += s.y*N;
  out[N-k2].x += s.x*N;
  out[N-k2].y += s.y*N;
}

/* Normalise the input to the range -1...+1 */

void normalise(tFFTTest *FFTTest) {
  real max=0;
  real *in  = (real *)FFTTest->in;
  real *out = (real *)FFTTest->out;
  int N = 2*FFTTest->N;
  int i;
  
  for (i=0; i<N; i++) {
    real r=in[i];
    if (r<0) r=-r;
    if (r>max) max=r;
  }
  max = 1/max;
  for (i=0; i<N; i++) {
    in[i] *= max;
    out[i] *= max;
  }
}

/* output the FFTTest arrays to the C file */

void output(FILE *fp, comp *in, int N, int nTest, char *name) {
  int k;
  
  fprintf(fp, "static const comp ffttest_%d_%d_%s[%d] = {\n",
    N, nTest, name, N
  );
  for (k=0; k<N-1; k++) {
    fprintf(fp, "  { %.18lf, %.18lf },\n", in[k].x, in[k].y);
  }
  fprintf(fp, "  { %.18lf, %.18lf }\n};\n\n", in[k].x, in[k].y);
}

/* write the test structure */

void output_tests(FILE *fp, int N) {
  int i;
  
  fprintf(fp, "static const tFFTTest FFTTest_%d[%d]={\n", N, nLFFTTests);
  for (i=0; i<nLFFTTests; ) {
    tFFTTest *FFTTest = &LFFTTests[i];

    fprintf(fp, "  { %d, \"%s\", (comp*)ffttest_%d_%d_in, "
      "(comp*)ffttest_%d_%d_out, %d }",
      FFTTest->N, FFTTest->name, N, i, N, i, FFTTest->flags);
    i++;
    if (i<nLFFTTests) fprintf(fp, ",");
    fprintf(fp, "\n");
  }
  fprintf(fp, "};\n\n");
  fprintf(fp, "extern const tFFTTests FFTTests_%d;\n", N);
  fprintf(fp, "const tFFTTests FFTTests_%d = {\n", N);
  fprintf(fp, " %d, %d, (tFFTTest*)FFTTest_%d\n", nLFFTTests, N, N);
  fprintf(fp, "};\n");
}