DSPRelated.com
Forums

Fixed-point FFT on TMS320C6416

Started by phil_levy1 November 6, 2002
Help !!!
I have been trying to do a complex FFT using the DSP_fft() function
provided by the TMS320C64x DSP Library (see TMS320C64x DSP Library
Programmers' Reference - SPRU565A). I have provided the "twiddle
factor" array by running the tw_fft16x16.exe program to produce a *.c
file, also provided in source-code by TI. I have left the "scale"
parameter for this proram at the default value (32765.5).

(Q1)Does this mean Q1.15 ??

My input data represents complex values (stored in the specified
real/imag ordering), and has been quantized to Q1.15 values. Even
when my input data represents a sine wave (created in MATLAB,
quantized to 16-bits in Q1.15 format), the output of the FFT appears
to be garbage.

(Q2) Does anyone have any experience with using (any of) theFFT
functions provided in this TMS320C64x DSP Library ??

Would really, appreciate any help ...:)

Phil



phil_levy1 schrieb:

> Help !!!
> I have been trying to do a complex FFT using the DSP_fft() function
> provided by the TMS320C64x DSP Library (see TMS320C64x DSP Library
> Programmers' Reference - SPRU565A). I have provided the "twiddle
> factor" array by running the tw_fft16x16.exe program to produce a *.c
> file, also provided in source-code by TI. I have left the "scale"
> parameter for this proram at the default value (32765.5).
>
> (Q1)Does this mean Q1.15 ??
>
> My input data represents complex values (stored in the specified
> real/imag ordering), and has been quantized to Q1.15 values. Even
> when my input data represents a sine wave (created in MATLAB,
> quantized to 16-bits in Q1.15 format), the output of the FFT appears
> to be garbage.
>
> (Q2) Does anyone have any experience with using (any of) theFFT
> functions provided in this TMS320C64x DSP Library ??
>
> Would really, appreciate any help ...:)
>
> Phil

Hello Phil,

assembly optimized FFT subroutines require some precautions:
- alignment of input / output buffer(s) in memory
- input and output data maybe presented in a special order (i.e. bit-reversed)
- buffer / FFT size constraints
- twiddle factor size and presentation

thus, I would recommend to check the *.h file or the source code of
the library for calling constraints. Maybe the output is just bit-reversed
i.e. you need to run a bit-reverse function to get in correct order?

Rgds. Phil
--
Dipl.-Inf.
Phil Alder - Field Application Engineer
------------------------------
DSPecialists GmbH
Rotherstr. 22
10245 Berlin, Germany http://www.DSPecialists.de
------------------------------
DSPecialists - Making the impossible work!


Phil-

> I have been trying to do a complex FFT using the DSP_fft() function
> provided by the TMS320C64x DSP Library (see TMS320C64x DSP Library
> Programmers' Reference - SPRU565A). I have provided the "twiddle
> factor" array by running the tw_fft16x16.exe program to produce a *.c
> file, also provided in source-code by TI. I have left the "scale"
> parameter for this proram at the default value (32765.5).
>
> (Q1)Does this mean Q1.15 ??

Yes, seems like it; i.e. 1 sign bit and 15 mantissa bits. The scale value is a
bit
weird, where 32765.5 comes from I'm not sure, I would have expected 32767.

I've attached some C code we've used before on C62xx to get very precise
fixed-point
FFT results. It's the "accurate" version, not the "fast" version, but it should
help
with understanding the optimized routine provided by TI.

Jeff Brower
DSP sw/hw engineer
Signalogic void FFTRadix2(int xy[], int n, short dir)
/*******************************************************************************

* NAME *
* FFTRadix2 *
* USAGE *
* This routine is C-callable and can be called as: *
* *
* void FFTRadix2(short xy[],short n,short w[]) *
* *
* xy[] --- input and output sequences (dim-n) (input/output) *
* n --- FFT size (input) *
* w[] --- FFT coefficients (dim-n/2) (input) *
* *
* ORIGINAL Copyright (C) Texas Instruments 1999 *
* *
* MODIFICATIONS Copyright (C) Signalogic 2001 *
* *
* DESCRIPTION *
* This routine is used to compute FFT of a complex sequence *
* of size n, a power of 2, with "decimation-in-frequency *
* decomposition" method, ie, the output is in bit-reversed order. *
* *
* 4 <= n <= 2048 *
* Each complex value has interleaved 16-bit real and imaginary parts. *
* Input xy and Coefficients w are 16 bit data. *
* sine and cosine data are stored in an integer lookup table as *
* sin(n*pi/2048) and cos(n*pi/2048) coefficients multiplied by 32767. *
* The sin/cos lookup table limits the maximum FFT size to 2048; this *
* can be changed by increasing/decreasing the table size and adjusting *
* the value of MAXFFT. Default value of MAXFFT is 11. *
* Rounding is used in the butterfly calculation to improve output *
* accuracy; can be removed for better performance. *
*
* TO USE TI'S EXACT METHOD, you need to add Bitrev_xtra.c to the *
* project, and then call digitrev_index() and bitrev_cplx() functions *

* prior to calling FFTRadix2(). However, please note that this method *
* requires a *huge* Index array to hold the pre-calculated bit-reverse *
* index values, either on the heap or on the stack (this is why we tend *
* to avoid it.) *
* *
* MEMORY NOTE: *
* For assembly language implementation, align xy and w on different word *
* boundaries to minimize memory bank hits. *
* Both input xy and coefficient w should be aligned on word boundary. *
*******************************************************************************/
{ int n1,n2,ie,ia,i,j,k,l;
short xt,yt;
short c,s;

int* wcosine;
int* wsine;
int skip;

int TR,TI; wcosine = (int*)&COSINE;
wsine = (int*)&SINTAB;

skip = 1 << (MAXFFT - ORDVAL);

n2 = n;
ie = 1;

for (k=n; k>1; k>>=1) {

n1 = n2;
n2 >>= 1;
ia = 0;

for (j=0; j<n2; j++) {

c = (short)-wcosine[skip*ia];
s = (short)wsine[skip*ia];

if (dir < 0) s = -s; // inverse FFT check

ia = ia + ie;

for (i=j; i < n; i += n1) {

l = i + n2;
xt = (short)xy[2*l] - (short)xy[2*i];
xy[2*i] = (short)xy[2*i] + (short)xy[2*l];
yt = (short)xy[2*l+1] - (short)xy[2*i+1];
xy[2*i+1] = (short)xy[2*i+1] + (short)xy[2*l+1];

TR = (c*xt - s*yt);

if (TR < 0) xy[2*l] = (TR - 16384) / 32767;
else xy[2*l] = (TR + 16384) / 32767;

TI = (c*yt + s*xt);

if (TI < 0) xy[2*l+1] = (TI - 16384) / 32767;
else xy[2*l+1] = (TI + 16384) / 32767;
}
}

ie <<= 1;
}
}