dsp_fft.h64

来自「dm642函数库」· H64 代码 · 共 532 行 · 第 1/3 页

H64
532
字号
;* ======================================================================== *;
;*  TEXAS INSTRUMENTS, INC.                                                 *;
;*                                                                          *;
;*  DSPLIB  DSP Signal Processing Library                                   *;
;*                                                                          *;
;*      Release:        Revision 1.04b                                      *;
;*      CVS Revision:   1.14    Sun Sep 29 03:32:20 2002 (UTC)              *;
;*      Snapshot date:  23-Oct-2003                                         *;
;*                                                                          *;
;*  This library contains proprietary intellectual property of Texas        *;
;*  Instruments, Inc.  The library and its source code are protected by     *;
;*  various copyrights, and portions may also be protected by patents or    *;
;*  other legal protections.                                                *;
;*                                                                          *;
;*  This software is licensed for use with Texas Instruments TMS320         *;
;*  family DSPs.  This license was provided to you prior to installing      *;
;*  the software.  You may review this license by consulting the file       *;
;*  TI_license.PDF which accompanies the files in this library.             *;
;* ------------------------------------------------------------------------ *;
;*          Copyright (C) 2003 Texas Instruments, Incorporated.             *;
;*                          All Rights Reserved.                            *;
;* ======================================================================== *;
;* ======================================================================== *;
;*  Assembler compatibility shim for assembling 4.30 and later code on      *;
;*  tools prior to 4.30.                                                    *;
;* ======================================================================== *;
;* ======================================================================== *;
;*  End of assembler compatibility shim.                                    *;
;* ======================================================================== *;
* ========================================================================= *
*                                                                           *
*   TEXAS INSTRUMENTS, INC.                                                 *
*                                                                           *
*   NAME                                                                    *
*       DSP_fft                                                             *
*                                                                           *
*                                                                           *
*   REVISION DATE                                                           *
*       16-Oct-2000                                                         *
*                                                                           *
*   USAGE                                                                   *
*       This routine is C-callable and can be called as:                    *
*                                                                           *
*       void DSP_fft(const short *w, int nsamp, short *x, short *y);        *
*                                                                           *
*       nsamp = length of DSP_fft in complex samples                        *
*       x     = pointer to complex data input, time domain                  *
*       w     = pointer to complex twiddle factors                          *
*       y     = pointer to complex data output, frequency domain            *
*                                                                           *
*   DESCRIPTION                                                             *
*       This code performs a Radix-4 FFT with digit reversal.  The code     *
*       uses a special ordering of twiddle factors and memory accesses      *
*       to improve performance in the presence of cache.  It operates       *
*       largely in-place, but the final digit-reversed output is written    *
*       out-of-place.                                                       *
*                                                                           *
*       This code requires a special sequence of twiddle factors stored     *
*       in Q.15 fixed-point format.  The following C code illustrates       *
*       one way to generate the desired twiddle-factor array:               *
*                                                                           *
*       #include <math.h>                                                   *
*                                                                           *
*       #ifndef PI                                                          *
*       # define PI (3.14159265358979323846)                                *
*       #endif                                                              *
*                                                                           *
*       short d2s(double d)                                                 *
*       {                                                                   *
*           d = floor(0.5 + d);  /* Explicit rounding to integer */         *
*           if (d >=  32767.0) return  32767;                               *
*           if (d <= -32768.0) return -32768;                               *
*           return (short)d;                                                *
*       }                                                                   *
*                                                                           *
*       void gen_twiddle(short *w, int n)                                   *
*       {                                                                   *
*           double M = 32767.5;                                             *
*           int i, j, k;                                                    *
*                                                                           *
*           for (j = 1, k = 0; j < n >> 2; j = j << 2)                      *
*           {                                                               *
*               for (i = 0; i < n >> 2; i += j << 1)                        *
*               {                                                           *
*                   w[k + 11] = d2s(M * cos(6.0 * PI * (i + j) / n));       *
*                   w[k + 10] = d2s(M * sin(6.0 * PI * (i + j) / n));       *
*                   w[k +  9] = d2s(M * cos(6.0 * PI * (i    ) / n));       *
*                   w[k +  8] = d2s(M * sin(6.0 * PI * (i    ) / n));       *
*                                                                           *
*                   w[k +  7] = d2s(M * cos(4.0 * PI * (i + j) / n));       *
*                   w[k +  6] = d2s(M * sin(4.0 * PI * (i + j) / n));       *
*                   w[k +  5] = d2s(M * cos(4.0 * PI * (i    ) / n));       *
*                   w[k +  4] = d2s(M * sin(4.0 * PI * (i    ) / n));       *
*                                                                           *
*                   w[k +  3] = d2s(M * cos(2.0 * PI * (i + j) / n));       *
*                   w[k +  2] = d2s(M * sin(2.0 * PI * (i + j) / n));       *
*                   w[k +  1] = d2s(M * cos(2.0 * PI * (i    ) / n));       *
*                   w[k +  0] = d2s(M * sin(2.0 * PI * (i    ) / n));       *
*                                                                           *
*                   k += 12;                                                *
*               }                                                           *
*           }                                                               *
*           w[2*n - 1] = w[2*n - 3] = w[2*n - 5] = 32767;                   *
*           w[2*n - 2] = w[2*n - 4] = w[2*n - 6] = 0;                       *
*       }                                                                   *
*                                                                           *
*   ASSUMPTIONS                                                             *
*       n must be a power of 4 and n >= 16 & n < 32768.                     *
*       FFT data x are aligned on a double word boundary, in real/imag      *
*       pairs, FFT twiddle factors w are also aligned on a double word      *
*       boundary in real/imaginary pairs.                                   *
*                                                                           *
*       Input FFT coeffs. are in signed Q.15 format.                        *
*       The memory Configuration is LITTLE ENDIAN.                          *
*       The complex data will be returned in natural order. This code is    *
*       uninteruptable, interupts are disabled on entry to the function and *
*       re-enabled on exit.                                                 *
*                                                                           *
*   MEMORY NOTE                                                             *
*       No bank conflict stalls occur in this code.                         *
*                                                                           *
*   TECHNIQUES                                                              *
*       A special sequence of coefficients. are used (as generated above)   *
*       to produce the DSP_fft. This collapses the inner 2 loops in the     *
*       taditional Burrus and Parks implementation Fortran Code.            *
*                                                                           *
*      The following C code represents an implementation of the Cooley      *
*      Tukey radix 4 DIF FFT. It accepts the inputs in normal order and     *
*      produces the outputs in digit reversed order. The natural C code     *
*      shown in this file on the other hand, accepts the inputs in nor-     *
*      mal order and produces the outputs in normal order.                  *
*                                                                           *
*      Several transformations have been applied to the original Cooley     *
*      Tukey code to produce the natural C code description shown here.     *
*      In order to understand these it would first be educational to        *
*      understand some of the issues involved in the conventional Cooley    *
*      Tukey FFT code.                                                      *
*                                                                           *
*       void radix4(int n, short x[], short wn[])                           *
*       {                                                                   *
*           int    n1,  n2,  ie,   ia1,  ia2, ia3;                          *
*           int    i0,  i1,  i2,    i3,    i, j,     k;                     *
*           short  co1, co2, co3,  si1,  si2, si3;                          *
*           short  xt0, yt0, xt1,  yt1,  xt2, yt2;                          *
*           short  xh0, xh1, xh20, xh21, xl0, xl1,xl20,xl21;                *
*                                                                           *
*           n2 = n;                                                         *
*           ie = 1;                                                         *
*           for (k = n; k > 1; k >>= 2)                                     *
*           {                                                               *
*               n1 = n2;                                                    *
*               n2 >>= 2;                                                   *
*               ia1 = 0;                                                    *
*                                                                           *
*               for (j = 0; j < n2; j++)                                    *
*               {                                                           *
*                    ia2 = ia1 + ia1;                                       *
*                    ia3 = ia2 + ia1;                                       *
*                                                                           *
*                    co1 = wn[2 * ia1    ];                                 *
*                    si1 = wn[2 * ia1 + 1];                                 *
*                    co2 = wn[2 * ia2    ];                                 *
*                    si2 = wn[2 * ia2 + 1];                                 *
*                    co3 = wn[2 * ia3    ];                                 *
*                    si3 = wn[2 * ia3 + 1];                                 *
*                    ia1 = ia1 + ie;                                        *
*                                                                           *
*                    for (i0 = j; i0< n; i0 += n1)                          *
*                    {                                                      *
*                        i1 = i0 + n2;                                      *
*                        i2 = i1 + n2;                                      *
*                        i3 = i2 + n2;                                      *
*                                                                           *
*                                                                           *
*                        xh0  = x[2 * i0    ] + x[2 * i2    ];              *
*                        xh1  = x[2 * i0 + 1] + x[2 * i2 + 1];              *
*                        xl0  = x[2 * i0    ] - x[2 * i2    ];              *
*                        xl1  = x[2 * i0 + 1] - x[2 * i2 + 1];              *

⌨️ 快捷键说明

复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?