📄 dsp_fft32x32.h64
字号:
;* ======================================================================== *;
;* TEXAS INSTRUMENTS, INC. *;
;* *;
;* DSPLIB DSP Signal Processing Library *;
;* *;
;* Release: Revision 1.04b *;
;* CVS Revision: 1.7 Sun Sep 29 03:32:21 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_fft32x32: Double Precision FFT *
* *
* USAGE *
* This routine is C-callable and can be called as: *
* *
* void DSP_fft32x32(const int * ptr_w, int npoints, *
* int * ptr_x, int * ptr_y ) ; *
* *
* ptr_w = input twiddle factors *
* npoints = number of points *
* ptr_x = transformed data reversed *
* ptr_y = linear transformed data *
* *
* (See the C compiler reference guide.) *
* *
* DESCRIPTION *
* The following code performs a mixed radix FFT for "npoints" which *
* is either a multiple of 4 or 2. It uses logN4 - 1 stages of radix4 *
* transform and performs either a radix2 or radix4 transform on the *
* last stage depending on "npoints". If "npoints" is a multiple of 4, *
* then this last stage is also a radix4 transform, otherwise it is a *
* radix2 transform. This program is available as a C compilable file *
* to automatically generate the twiddle factors "twiddle_split.c" *
* *
* Generate special vector of twiddle factors *
* *
* for (j=1, k=0; j < npoints>>2; j = j <<2 ) *
* { *
* for (i=0; i < npoints>>2; i += j) *
* { *
* theta1 = 2*PI*i/npoints; *
* x_t = M*cos(theta1); *
* y_t = M*sin(theta1); *
* ptr_w[k+1] = (int) x_t; *
* if (x_t >= M) ptr_w[k+1] = 0x7fffffff; *
* ptr_w[k+0] = (int) y_t; *
* if (y_t >= M) ptr_w[k+0] = 0x7fffffff; *
* *
* theta2 = 4*PI*i/npoints; *
* x_t = M*cos(theta2); *
* y_t = M*sin(theta2); *
* ptr_w[k+3] = (int) x_t; *
* *
* if (x_t >= M) ptr_w[k+3] = 0x7fffffff; *
* ptr_w[k+2] = (int) y_t; *
* if (y_t >= M) ptr_w[k+2] = 0x7fffffff; *
* *
* theta3 = 6*PI*i/npoints; *
* x_t = M*cos(theta3); *
* y_t = M*sin(theta3); *
* ptr_w[k+5] = (int) x_t; *
* if (x_t >= M) ptr_w[k+5] = 0x7fffffff; *
* ptr_w[k+4] = (int) y_t; *
* if (y_t >= M) ptr_w[k+4] = 0x7fffffff; *
* k += 6; *
* } *
* } *
* *
* *
* ASSUMPTIONS *
* This code works for both "npoints" a multiple of 2 or 4. *
* The arrays 'x[]', 'y[]', and 'w[]' all must be aligned on a *
* double-word boundary for the "optimized" implementations. *
* The input and output data are complex, with the real/imaginary *
* components stored in adjacent locations in the array. The real *
* components are stored at even array indices, and the imaginary *
* components are stored at odd array indices. The input, twiddle *
* factors are in 32 bit precision. The 32 by 32 multiplies are *
* done with a 1.5 bit loss in accuracy. This comes about because *
* the contribution of the low sixteen bits to the 32 bit result *
* is not computed. In addition the contribution of the low * high *
* term is shifted by 16 as opposed to 15, for a loss 0f 0.5 bits *
* after rounding. To illustrate real part of complex multiply of: *
* (X + jY) ( C + jS) = *
* *
* _mpyhir(si10 , yt1_0) + _mpyhir(co10 , xt1_0) + *
* (((MPYLUHS(si10,yt1_0) + MPYLUHS(co10, xt1_0) *
* + 0x8000) >> 16) << 1) *
* *
* The intrinsic C version of this code performs this function as: *
* *
* _mpyhir(si10 , yt1_0) + _mpyhir(co10 , xt1_0) + *
* (_dotprsu2(yt1_0xt1_0, si10co10) << 1); *
* *
* *
* where the functions _mpyhir, MPYLUHS are as follows: *
* *
* #define _mpyhir(x,y) \ *
* (((int)((short)(x>>16)*(unsigned short)(y&0x0000FFFF)+0x4000) >> 15) *
* + \ ((int)((short)(x >> 16) * (short)((y) >> 16)) << 1)) *
* *
* #define MPYLUHS(x,y) \ *
* ( (int) ((unsigned short)(x & 0x0000FFFF) * (short) (y >> 16)) ) *
* *
* *
* TECHNIQUES *
* 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; *
* *
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -