# Twiddle Factor calculation for FFT-Radix 4

Started by December 20, 2007
Dear all,

I am running a program for DSK6711 that calculates the FFT-radix 4. As is shown in the attached code, data length is defined as a constant value of 128 samples (NUMDATA)and twiddle arrays (TR,TI), which values are defined in twiddleI.h and twiddleR.h, are 64 samples each one (all files are avalaible in the CD "DSP System Design Using the TMS320C6000"). At this point the program works ok!

My next step is to increase the number of points of FFT but before do it, I want to know how was the twiddle factor generated. To understand it I removed twiddleI.h and twiddleR.h files from my project and replaced code in the main function with the following loop:

for (i=0;i {
W4[i].real = cos(2*PI*i/NUMDATA);
W4[i].imag = sin(2*PI*i/NUMDATA);
}

but now the output values are wrong for the same 128 data record. So I have compared the real and imaginary values generated with loop with values defined in twiddleI/R.h files and they are different and I don't understand why!

So my questions are:

(1) How were generated TR and TI values in header files to make correct the FFT calculation?
(2) How can be possible that TR,TI were 64 samples length when in Chassaing (page 204 of TMS320C6x book) the length of the set up for a FFT radix 4 is N/RADIX length (whichs means in my case 128/42!)?

The source code is listed below.

Any help would be really appreciated! =)

Benjamin Sanchez
Electronic Department
Technical University of Catalonia UPC
Barcelona (SPAIN)

/******************/
/******************/

#include
#include "twiddleR.h"
#include "twiddleI.h"

#define NUMDATA 128 // number of real data samples 128
#define NUMPOINTS 64 // number of point in the DFT,
// NUMDATA/2
#define PI 3.141592653589793 // definition of pi //

typedef struct { //data type for the radix-4 twiddle factors
short imag;
short real;
} COEFF;

typedef struct {
short real;
short imag;
} COMPLEX;
#pragma DATA_ALIGN(x,NUMPOINTS);

COMPLEX x[NUMPOINTS+1];// array of complex DFT data
COEFF W4[NUMPOINTS];
short g[NUMDATA];
COMPLEX A[NUMPOINTS]; // array of complex A coefficients //
COMPLEX B[NUMPOINTS]; // array of complex B coefficients //
COMPLEX G[2*NUMPOINTS]; // array of complex DFT result //
unsigned short IIndex[NUMPOINTS], JIndex[NUMPOINTS];
int count;

int magR[NUMDATA];
int magI[NUMDATA];
int mag[NUMDATA];

...

int data[NUMDATA]={...}; //signal waveform

int output[NUMDATA/2];

void main(void)
{
int i,j,k;
short tr[NUMPOINTS], ti[NUMPOINTS];

//Read Twiddle factors to COMPLEX array and make Q-15;
make_q15(tr, TR, NUMPOINTS); //Data in Header files from Matlab
make_q15(ti, TI, NUMPOINTS);

for(i=0;i {
W4[i].real = tr[i];
W4[i].imag = ti[i];
}

// Initialize A,B arrays //
for(k=0; k {
A[k].imag = (short)(16383.0*(-cos(2*PI/(double)(2*NUMPOINTS)*(double)k)));
A[k].real = (short)(16383.0*(1.0 - sin(2*PI/(double)(2*NUMPOINTS)*(double)k)));
B[k].imag = (short)(16383.0*(cos(2*PI/(double)(2*NUMPOINTS)*(double)k)));
B[k].real = (short)(16383.0*(1.0 + sin(2*PI/(double)(2*NUMPOINTS)*(double)k)));
}

// Initialize tables for FFT digit reversal function //
R4DigitRevIndexTableGen(NUMPOINTS, &count, IIndex, JIndex);

for(i=0;i g[i]a[i];

// Call FFT algorithm //
fft();

for (k=0; k {
magR[k] = (G[k].real*G[k].real) << 1;
magI[k] = (G[k].imag*G[k].imag) << 1;
mag[k]=magR[k]+magI[k];
}
}

void fft()
{
int n;

for (n=0; n {
x[n].imag = g[2*n + 1];
x[n].real = g[2*n];
}

digit_reverse((int *)x, IIndex, JIndex, count);

x[NUMPOINTS].real = x[0].real;
x[NUMPOINTS].imag = x[0].imag;

split1(NUMPOINTS, x, A, B, G);

G[NUMPOINTS].real = x[0].real - x[0].imag;
G[NUMPOINTS].imag = 0;

for (n=1; n {
G[2*NUMPOINTS-n].real = G[n].real;
G[2*NUMPOINTS-n].imag = -G[n].imag;
}

}

void radix4(int n, short x[], short w[])
{
int n1, n2, ie, ia1, ia2, ia3, i0, i1, i2, i3, j, k;
short t, r1, r2, s1, s2, co1, co2, co3, si1, si2, si3;
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 = w[ia1 * 2 + 1];
si1 = w[ia1 * 2];
co2 = w[ia2 * 2 + 1];
si2 = w[ia2 * 2];
co3 = w[ia3 * 2 + 1];
si3 = w[ia3 * 2];
ia1 = ia1 + ie;
for (i0 = j; i0 < n; i0 += n1)
{
i1 = i0 + n2;
i2 = i1 + n2;
i3 = i2 + n2;
r1 = x[2 * i0] + x[2 * i2];
r2 = x[2 * i0] - x[2 * i2];
t = x[2 * i1] + x[2 * i3];
x[2 * i0] = r1 + t;
r1 = r1 - t;
s1 = x[2 * i0 + 1] + x[2 * i2 + 1];
s2 = x[2 * i0 + 1] - x[2 * i2 + 1];
t = x[2 * i1 + 1] + x[2 * i3 + 1];
x[2 * i0 + 1] = s1 + t;
s1 = s1 - t;
x[2 * i2] = (r1 * co2 + s1 * si2) >> 15;
x[2 * i2 + 1] = (s1 * co2-r1 * si2)>>15;
t = x[2 * i1 + 1] - x[2 * i3 + 1];
r1 = r2 + t;
r2 = r2 - t;
t = x[2 * i1] - x[2 * i3];
s1 = s2 - t;
s2 = s2 + t;
x[2 * i1] = (r1 * co1 + s1 * si1) >>15;
x[2 * i1 + 1] = (s1 * co1-r1 * si1)>>15;
x[2 * i3] = (r2 * co3 + s2 * si3) >>15;
x[2 * i3 + 1] = (s2 * co3-r2 * si3)>>15;
}
}
ie <<= 2;
}
}

void digit_reverse(int *yx, unsigned short *JIndex, unsigned short *IIndex, int count)
{
int i;
unsigned short I, J;
int YXI, YXJ;
for (i = 0; i {
I = IIndex[i];
J = JIndex[i];
YXI = yx[I];
YXJ = yx[J];
yx[J] = YXI;
yx[I] = YXJ;
}
}

void R4DigitRevIndexTableGen(int n, int *count, unsigned short *IIndex, unsigned short
*JIndex)
{
int j, n1, k, i;
j = 1;
n1 = n - 1;
*count = 0;
for(i=1; i<=n1; i++)
{
if(i < j)
{
IIndex[*count] = (unsigned short)(i-1);
JIndex[*count] = (unsigned short)(j-1);
*count = *count + 1;
}
k = n >> 2;
while(k*3 < j)
{
j = j - k*3;
k = k >> 2;
}
j = j + k;
}
}
void split1(int N, COMPLEX *X, COMPLEX *A, COMPLEX *B, COMPLEX *G)
{
int k;
int Tr, Ti;

for (k=0; k {
Tr = (int)X[k].real * (int)A[k].real - (int)X[k].imag * (int)A[k].imag + (int)X[N-k].real * (int)B[k].real + (int)X[N-k].imag * (int)B[k].imag;

G[k].real = (short)(Tr>>15);

Ti = (int)X[k].imag * (int)A[k].real + (int)X[k].real * (int)A[k].imag + (int)X[N-k].real * (int)B[k].imag - (int)X[N-k].imag * (int)B[k].real;

G[k].imag = (short)(Ti>>15);
}
}

void make_q15(short out[], float in[], int N)
{
int i;

for(i=0;i {
out[i]=0x7fff*in[i]; //Convert to Q-15, good approximate
}
}
/******************/
twiddleR.h
/******************/

float TR[64] ={0.9999,
0.9952,
0.9808,
0.9569,
0.9239,
0.8819,
0.8315,
0.7730,
0.7071,
0.6344,
0.5556,
0.4714,
0.3827,
0.2903,
0.1951,
0.0980,
0.0000,
-0.0980,
-0.1951,
-0.2903,
-0.3827,
-0.4714,
-0.5556,
-0.6344,
-0.7071,
-0.7730,
-0.8315,
-0.8819,
-0.9239,
-0.9569,
-0.9808,
-0.9952,
-1.0000,
-0.9952,
-0.9808,
-0.9569,
-0.9239,
-0.8819,
-0.8315,
-0.7730,
-0.7071,
-0.6344,
-0.5556,
-0.4714,
-0.3827,
-0.2903,
-0.1951,
-0.0980,
0.0000,
0.0980,
0.1951,
0.2903,
0.3827,
0.4714,
0.5556,
0.6344,
0.7071,
0.7730,
0.8315,
0.8819,
0.9239,
0.9569,
0.9808,
0.9952};

/******************/
twiddleI.h
/******************/

float TI[64] = {0.0000,
0.0980,
0.1951,
0.2903,
0.3827,
0.4714,
0.5556,
0.6344,
0.7071,
0.7730,
0.8315,
0.8819,
0.9239,
0.9569,
0.9808,
0.9952,
0.9999,
0.9952,
0.9808,
0.9569,
0.9239,
0.8819,
0.8315,
0.7730,
0.7071,
0.6344,
0.5556,
0.4714,
0.3827,
0.2903,
0.1951,
0.0980,
0.0000,
-0.0980,
-0.1951,
-0.2903,
-0.3827,
-0.4714,
-0.5556,
-0.6344,
-0.7071,
-0.7730,
-0.8315,
-0.8819,
-0.9239,
-0.9569,
-0.9808,
-0.9952,
-1.0000,
-0.9952,
-0.9808,
-0.9569,
-0.9239,
-0.8819,
-0.8315,
-0.7730,
-0.7071,
-0.6344,
-0.5556,
-0.4714,
-0.3827,
-0.2903,
-0.1951,
-0.0980};