fft.c

用 Hitex 工具软件开发 stm32 的例子

/********************************************************************
 * Project:    Tasking-STM32-Stick
 * File:       fft.c
 *
 * System:     Cortex ARMv7 32 Bit (STM32FRT)
 * Compiler:   Tasking Altuim VX Toolchain v2.01
 *
 * Date:       2007-08-20
 * Author:     [e-mail]
 *
 * Rights:     Hitex Development Tools GmbH
 *             Greschbachstr. 12
 *             D-76229 Karlsruhe
 ********************************************************************
 * Description:
 *
 * This file is part of the Tasking Example chain
 * The code is based on usage of the STmicro library functions
 * This is a small implementation of different features
 * The application runs in ARM mode with high optimization level.
 * Fixed-point Fast Fourier Transform
 *
 *    fft()       perform FFT or inverse FFT
 *    fix_mpy()       perform fixed-point multiplication.
 *    Sinewave[1024]  sinewave normalized to 32767 (= 1.0).
 *
 *    All data are fixed-point short integers, in which
 *    -32768 to +32768 represent -1.0 to +1.0. Integer arithmetic
 *    is used for speed, instead of the more natural floating-point.

 *    For the forward FFT (time -> freq), fixed scaling is
 *    performed to prevent arithmetic overflow, and to map a 0dB
 *    sine/cosine wave (i.e. amplitude = 32767) to two -6dB freq
 *    coefficients; the one in the lower half is reported as 0dB
 *    by fix_loud(). The return value is always 0.
 *
 *    Timing on a Macintosh PowerBook 180.... (using Symantec C6.0)
 *                      fix_fft (1024 points)             8 ticks
 *                      fft (1024 points - Using SANE)  112 Ticks
 *                      fft (1024 points - Using FPU)    11
 *
 ********************************************************************
 * History:
 *
 *    Revision 1.0    2007/08/20      Gn
 *    Initial revision
 ********************************************************************
 * This is a preliminary version.
 *
 * WARRANTY:  HITEX warrants that the media on which the SOFTWARE is
 * furnished is free from defects in materials and workmanship under
 * normal use and service for a period of ninety (90) days. HITEX entire
 * liability and your exclusive remedy shall be the replacement of the
 * SOFTWARE if the media is defective. This Warranty is void if failure
 * of the media resulted from unauthorized modification, accident, abuse,
 * or misapplication.
 *
 * DISCLAIMER:  OTHER THAN THE ABOVE WARRANTY, THE SOFTWARE IS FURNISHED
 * "AS IS" WITHOUT WARRANTY OF ANY KIND. HITEX DISCLAIMS ALL OTHER WARRANTIES,
 * EITHER EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO IMPLIED WARRANTIES
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 * NEITHER HITEX NOR ITS AFFILIATES SHALL BE LIABLE FOR ANY DAMAGES ARISING
 * OUT OF THE USE OF OR INABILITY TO USE THE SOFTWARE, INCLUDING DAMAGES FOR
 * LOSS OF PROFITS, BUSINESS INTERRUPTION, OR ANY SPECIAL, INCIDENTAL, INDIRECT
 * OR CONSEQUENTIAL DAMAGES EVEN IF HITEX HAS BEEN ADVISED OF THE POSSIBILITY
 * OF SUCH DAMAGES.
 ********************************************************************/


/* FIX_MPY() - fixed-point multiplication macro.
   This macro is a statement, not an expression (uses asm).
   Scaling ensures that 32767*32767 = 32767. */

#include "FFTfunc.h"
#include "FFT_tab.h"

short fix_mpy(short a, short b);

/**
   fft() - perform fast Fourier transform.
   fr[n],fi[n] are real,imaginary arrays, INPUT AND RESULT.
*/
void fft(short fr[], short fi[])
{
   int mr,nn,i,j,p,k,istep,m,n;
   short qr,qi,tr,ti,wr,wi;

   n = N_SAMPLE;
   mr = 0;
   nn = n - 1;

   /* decimation in time - re-order data */
   for(m=1; m<=nn; ++m)
   {
      p = n;
      do {
            p >>= 1;
      } while(mr+p > nn);
      mr = (mr & (p-1)) + p;

      if(mr <= m) continue;
      tr = fr[m];
      fr[m] = fr[mr];
      fr[mr] = tr;
   }

   p = 1;
   k = LOG2_N_WAVE-1;
/**
   fixed scaling, for proper normalization -
   there will be log2(n) passes, so this
   results in an overall factor of 1/n,
   distributed to maximize arithmetic accuracy. */
   while(p < n)
   {
      /* it may not be obvious, but the shift will be performed
         on each data point exactly once, during this pass. */
      istep = p << 1;
      for(m=0; m<p; ++m)
      {
         j = m << k;
         /* 0 <= j < N_WAVE/2 */
         wr =  sine_table[j+N_WAVE/4];
         wi = -sine_table[j];
         wr >>= 1;
         wi >>= 1;
         for(i=m; i<n; i+=istep)
         {
            j = i + p;
            tr = fix_mpy(wr,fr[j]) - fix_mpy(wi,fi[j]);
            ti = fix_mpy(wr,fi[j]) + fix_mpy(wi,fr[j]);
            qr = fr[i];
            qi = fi[i];
            qr >>= 1;
            qi >>= 1;
            fr[j] = qr - tr;
            fi[j] = qi - ti;
            fr[i] = qr + tr;
            fi[i] = qi + ti;
         }
      }
      --k;
      p = istep;
   }
   return;
}


/**
   fix_mpy() - fixed-point multiplication
*/
short fix_mpy(short a, short b)
{
   FIX_MPY(a,a,b);
   return a;
}

void Norming (short fr[], short faktor, short dividor)
{
   int n;
   for (n = 0; n < N_SAMPLE; n++)
      NORMING(fr[n]);

   for (n = 0; n < N_SAMPLE; n++)
      fr[n] = (fr[n]* faktor) / dividor;
}