dsp_fft16x16r.h64

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

H64
689
字号
;* ======================================================================== *;
;*  TEXAS INSTRUMENTS, INC.                                                 *;
;*                                                                          *;
;*  DSPLIB  DSP Signal Processing Library                                   *;
;*                                                                          *;
;*      Release:        Revision 1.04b                                      *;
;*      CVS Revision:   1.9     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_fft16x16r : 16x16 FFT with rounding and scaling                *
*                                                                            *
*                                                                            *
*   REVISION DATE                                                            *
*       23-May-2002                                                          *
*                                                                            *
*   USAGE                                                                    *
*         This routine is C-callable and can be called as:                   *
*                                                                            *
*         void DSP_fft16x16r                                                 *
*         (                                                                  *
*             int                       n,                                   *
*             short *restrict           ptr_x,                               *
*             const short *restrict     ptr_w,                               *
*             unsigned char *restrict   brev,                                *
*             short *restrict           ptr_y,                               *
*             int                       radix,                               *
*             int                       offset,                              *
*             int                       nmax                                 *
*         )                                                                  *
*                                                                            *
*         N      = length of fft in complex samples, power of 2 <=16384      *
*         ptr_x  = pointer to complex data input                             *
*         ptr_w  = pointer to complex twiddle factor (see below)             *
*         brev   = pointer to bit reverse table containing 64 entries        *
*         n_min  = smallest fft butterfly used in computation                *
*                  used for decomposing fft into subffts, see notes          *
*         offset = index in complex samples of sub-fft from start of main ff *
*         n_max  = size of main fft in complex samples                       *
*                                                                            *
*   DESCRIPTION                                                              *
*        The benchmark performs a mixed radix forwards fft using             *
*        a special sequence of coefficients generated in the following       *
*        way:                                                                *
*                                                                            *
*          -* generate vector of twiddle factors for optimized algorithm *-  *
*         void tw_gen(short * w, int N)                                      *
*         {                                                                  *
*           int j, k;                                                        *
*           double x_t, y_t, theta1, theta2, theta3;                         *
*           const double PI = 3.141592654, M = 32767.0;                      *
*                                        -* M is 16383 for scale by 4 *-     *
*                                                                            *
*           for (j=1, k=0; j <= N>>2; j = j<<2)                              *
*           {                                                                *
*               for (i=0; i < N>>2; i+=j)                                    *
*               {                                                            *
*                   theta1 = 2*PI*i/N;                                       *
*                   x_t = M*cos(theta1);                                     *
*                   y_t = M*sin(theta1);                                     *
*                   w[k]   =  (short)x_t;                                    *
*                   w[k+1] =  (short)y_t;                                    *
*                                                                            *
*                   theta2 = 4*PI*i/N;                                       *
*                   x_t = M*cos(theta2);                                     *
*                   y_t = M*sin(theta2);                                     *
*                   w[k+2] =  (short)x_t;                                    *
*                   w[k+3] =  (short)y_t;                                    *
*                                                                            *
*                   theta3 = 6*PI*i/N;                                       *
*                   x_t = M*cos(theta3);                                     *
*                   y_t = M*sin(theta3);                                     *
*                   w[k+4] =  (short)x_t;                                    *
*                   w[k+5] =  (short)y_t;                                    *
*                   k+=6;                                                    *
*               }                                                            *
*           }                                                                *
*         }                                                                  *
*        This redundent set of twiddle factors is size 2*N short samples.    *
*        As pointed out later dividing these twiddle factors by 2 will give  *
*        an effective divide by 4 at each stage to guarentee no overflow.    *
*        The function is accurate to about 68dB of signal to noise ratio     *
*        to the DFT function below:                                          *
*                                                                            *
*         void dft(int n, short x[], short y[])                              *
*         {                                                                  *
*            int k,i, index;                                                 *
*            const double PI = 3  4159654;                                   *
*            short * p_x;                                                    *
*            double arg, fx_0, fx_1, fy_0, fy_1, co, si;                     *
*                                                                            *
*            for(k = 0; k<n; k++)                                            *
*            {                                                               *
*              p_x = x;                                                      *
*              fy_0 = 0;                                                     *
*              fy_1 = 0;                                                     *
*              for(i=0; i<n; i++)                                            *
*              {                                                             *
*                fx_0 = (double)p_x[0];                                      *
*                fx_1 = (double)p_x[1];                                      *
*                p_x += 2;                                                   *
*                index = (i*k) % n;                                          *
*                arg = 2*PI*index/n;                                         *
*                co = cos(arg);                                              *
*                si = -sin(arg);                                             *
*                fy_0 += ((fx_0 * co) - (fx_1 * si));                        *
*                fy_1 += ((fx_1 * co) + (fx_0 * si));                        *
*              }                                                             *
*              y[2*k] = (short)2*fy_0/sqrt(N);                               *
*              y[2*k+1] = (short)2*fy_1/sqrt(N);                             *
*            }                                                               *
*         }                                                                  *
*        Scaling takes place at each stage except the last one.              *
*        This is a divide by 2 to prevent overflow. All shifts are rounded t *
*        reduce truncation noise power by 3dB.                               *
*        The function takes the table and input data and calculates the fft  *
*        producing the frequency domain data in the Y array.                 *
*        As the fft allows every input point to effect every output point in *
*        a cache based system such as the c6211, this causes cache thrashing *
*        This is mitigated by allowing the main fft of size N to be divided  *
*        into several steps, allowing as much data reuse as possible.        *
*                                                                            *
*        For example the following function:                                 *
*                                                                            *
*        DSP_fft16x16r  (1024, &x_asm[0],&w[0],y_asm,brev,4,  0,1024);       *
*                                                                            *
*        is equvalent to:                                                    *
*                                                                            *
*        DSP_fft16x16r  (1024,&x_asm[2*0],  &w[0]  ,y_asm,brev,256, 0,1024); *
*        DSP_fft16x16r  (256, &x_asm[2*0],  &w[2*768],y_asm,brev,4, 0,1024); *
*        DSP_fft16x16r (256, &x_asm[2*256],&w[2*768],y_asm,brev,4, 256,1024); *
*        DSP_fft16x16r (256, &x_asm[2*512],&w[2*768],y_asm,brev,4, 512,1024); *
*        DSP_fft16x16r (256, &x_asm[2*768],&w[2*768],y_asm,brev,4, 768,1024); *
*                                                                            *
*        Notice how the 1st fft function is called on the entire 1K data set *
*        it covers the 1st pass of the fft until the butterfly size is 256.  *
*        The following 4 ffts do 256 pt ffts 25% of the size. These continue *
*        down to the end when the buttefly is of size 4. The use an index to *
*        the main twiddle factor array of 0.75*2*N. This is because the      *
*        twiddle factor array is composed of successively decimated versions *
*        of the main array.                                                  *
*                                                                            *
*        N not equal to a power of 4 can be used, i.e. 512. In this case to  *
*        decompose the fft the following would be needed :                   *
*                                                                            *
*        DSP_fft16x16r   (512, &x_asm[0],&w[0],y_asm,brev,2,  0,512);        *
*                                                                            *
*        is equvalent to:                                                    *
*                                                                            *

⌨️ 快捷键说明

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