8 point radix-2 DIF FFT (Sande- Tukey) algorithm

Status
Not open for further replies.
T

tarjina

Junior Member level 3
Joined
Jun 7, 2012
Messages
30
Helped
0
Reputation
0
Reaction score
0
Trophy points
1,286
Activity points
1,536
Hello,
I have a project submission on DIF FFT in C code. I have tried so many things but can not get the code to work.
Code: 
#include "fft.h"
int fix_fft(fixed fr[], fixed fi[], int m, int inverse)
{
	int mr,nn,i,j,l,k,istep, n, scale, shift,twid;

	fixed qr,qi;		//even input
	fixed tr,ti;		//odd input
	fixed wr,wi;		//twiddle factor
	fixed qr0,qi0,a,b,c,d;
        n = 1<<m;

	if(n > N_WAVE) return -1;

	mr = 0;
	nn = n - 1;
	scale = 0;

	l = n>>1;
	k = LOG2_N_WAVE-m;
	while(l > 0)
	{
		if(inverse)
		{
shift = 0;
			for(i=0; i<n; ++i)
			{
				j = fr[i];
				if(j < 0) j = -j;

				m = fi[i];
				if(m < 0) m = -m;

				if(j > 16383 || m > 16383)
				{
					shift = 1;
					++scale;
					break;
				}
			}
		}
		else {shift =1;}
        		istep = l << 1;		
		for(m=0; m<l; ++m)
		{
			j = m << k;
			wr =  Sinewave[j+N_WAVE/4];
			wi = -Sinewave[j];
			if(inverse) wi = -wi;
			if(shift)
			{
				wr >>= 1;
				wi >>= 1;
			}
			for(i=m; i<n; i+=istep)
			{
				 
				j = i + l;
				qr = fr[i];
				qi = fi[i];
				tr = qr - fr[j];
				ti = qi - fi[j];
				qr0 = qr + fr[j];
				qi0 = qi + fi[j];

				a = fix_mpy(wr,tr); 
				b = fix_mpy(wi,ti);
				c = fix_mpy(wr,ti) ;
				d = fix_mpy(wi,tr);

				if(shift)
				{
					qr0 >>= 1;
					qi0 >>= 1;
				}

				fr[j] = a-b;
				fi[j] = c+d;
				fr[i] = qr0;
				fi[i] = qi0;
			}
		}
		++k;
		l >>= 1;
	}
	/* decimation in time - re-order data */
	for(m=1; m<=nn; ++m) {
		l = n;
		do{
			l >>= 1;
		}while(mr+l > nn);
		mr = (mr & (l-1)) + l;

		if(mr <= m) continue;
		tr = fr[m];
		fr[m] = fr[mr];
		fr[mr] = tr;

		ti = fi[m];
		fi[m] = fi[mr];
		fi[mr] = ti;
	}

	return scale;
}
fixed fix_mpy(fixed a, fixed b)
{
	FIX_MPY(a,a,b);
	return a;
}





fixed Sinewave[1024] = {
		0,    201,    402,    603,    804,   1005,   1206,   1406,
		1607,   1808,   2009,   2209,   2410,   2610,   2811,   3011,
		3211,   3411,   3611,   3811,   4011,   4210,   4409,   4608,
		4807,   5006,   5205,   5403,   5601,   5799,   5997,   6195,
		6392,   6589,   6786,   6982,   7179,   7375,   7571,   7766,
		7961,   8156,   8351,   8545,   8739,   8932,   9126,   9319,
		9511,   9703,   9895,  10087,  10278,  10469,  10659,  10849,
		11038,  11227,  11416,  11604,  11792,  11980,  12166,  12353,
		12539,  12724,  12909,  13094,  13278,  13462,  13645,  13827,
		14009,  14191,  14372,  14552,  14732,  14911,  15090,  15268,
		15446,  15623,  15799,  15975,  16150,  16325,  16499,  16672,
		16845,  17017,  17189,  17360,  17530,  17699,  17868,  18036,
		18204,  18371,  18537,  18702,  18867,  19031,  19194,  19357,
		19519,  19680,  19840,  20000,  20159,  20317,  20474,  20631,
		20787,  20942,  21096,  21249,  21402,  21554,  21705,  21855,
		22004,  22153,  22301,  22448,  22594,  22739,  22883,  23027,
		23169,  23311,  23452,  23592,  23731,  23869,  24006,  24143,
		24278,  24413,  24546,  24679,  24811,  24942,  25072,  25201,
		25329,  25456,  25582,  25707,  25831,  25954,  26077,  26198,
		26318,  26437,  26556,  26673,  26789,  26905,  27019,  27132,
		27244,  27355,  27466,  27575,  27683,  27790,  27896,  28001,
		28105,  28208,  28309,  28410,  28510,  28608,  28706,  28802,
		28897,  28992,  29085,  29177,  29268,  29358,  29446,  29534,
		29621,  29706,  29790,  29873,  29955,  30036,  30116,  30195,
		30272,  30349,  30424,  30498,  30571,  30643,  30713,  30783,
		30851,  30918,  30984,  31049,
		31113,  31175,  31236,  31297,
		31356,  31413,  31470,  31525,  31580,  31633,  31684,  31735,
		31785,  31833,  31880,  31926,  31970,  32014,  32056,  32097,
		32137,  32176,  32213,  32249,  32284,  32318,  32350,  32382,
		32412,  32441,  32468,  32495,  32520,  32544,  32567,  32588,
		32609,  32628,  32646,  32662,  32678,  32692,  32705,  32717,
		32727,  32736,  32744,  32751,  32757,  32761,  32764,  32766,
		32767,  32766,  32764,  32761,  32757,  32751,  32744,  32736,
		32727,  32717,  32705,  32692,  32678,  32662,  32646,  32628,
		32609,  32588,  32567,  32544,  32520,  32495,  32468,  32441,
		32412,  32382,  32350,  32318,  32284,  32249,  32213,  32176,
		32137,  32097,  32056,  32014,  31970,  31926,  31880,  31833,
		31785,  31735,  31684,  31633,  31580,  31525,  31470,  31413,
		31356,  31297,  31236,  31175,  31113,  31049,  30984,  30918,
		30851,  30783,  30713,  30643,  30571,  30498,  30424,  30349,
		30272,  30195,  30116,  30036,  29955,  29873,  29790,  29706,
		29621,  29534,  29446,  29358,  29268,  29177,  29085,  28992,
		28897,  28802,  28706,  28608,  28510,  28410,  28309,  28208,
		28105,  28001,  27896,  27790,  27683,  27575,  27466,  27355,
		27244,  27132,  27019,  26905,  26789,  26673,  26556,  26437,
		26318,  26198,  26077,  25954,  25831,  25707,  25582,  25456,
		25329,  25201,  25072,  24942,  24811,  24679,  24546,  24413,
		24278,  24143,  24006,  23869,  23731,  23592,  23452,  23311,
		23169,  23027,  22883,  22739,  22594,  22448,  22301,  22153,
		22004,  21855,  21705,  21554,  21402,  21249,  21096,  20942,
		20787,  20631,  20474,  20317,  20159,  20000,  19840,  19680,
		19519,  19357,  19194,  19031,  18867,  18702,  18537,  18371,
		18204,  18036,  17868,  17699,  17530,  17360,  17189,  17017,
		16845,  16672,  16499,  16325,  16150,  15975,  15799,  15623,
		15446,  15268,  15090,  14911,  14732,  14552,  14372,  14191,
		14009,  13827,  13645,  13462,  13278,  13094,  12909,  12724,
		12539,  12353,  12166,  11980,  11792,  11604,  11416,  11227,
		11038,  10849,  10659,  10469,  10278,  10087,   9895,   9703,
		9511,   9319,   9126,   8932,   8739,   8545,   8351,   8156,
		7961,   7766,   7571,   7375,   7179,   6982,   6786,   6589,
		6392,   6195,   5997,   5799,   5601,   5403,   5205,   5006,
		4807,   4608,   4409,   4210,   4011,   3811,   3611,   3411,
		3211,   3011,   2811,   2610,   2410,   2209,   2009,   1808,
		1607,   1406,   1206,   1005,    804,    603,    402,    201,
		0,   -201,   -402,   -603,   -804,  -1005,  -1206,  -1406,
		-1607,  -1808,  -2009,  -2209,  -2410,  -2610,  -2811,  -3011,
		-3211,  -3411,  -3611,  -3811,  -4011,  -4210,  -4409,  -4608,
		-4807,  -5006,  -5205,  -5403,  -5601,  -5799,  -5997,  -6195,
		-6392,  -6589,  -6786,  -6982,  -7179,  -7375,  -7571,  -7766,
		-7961,  -8156,  -8351,  -8545,  -8739,  -8932,  -9126,  -9319,
		-9511,  -9703,  -9895, -10087, -10278, -10469, -10659, -10849,
		-11038, -11227, -11416, -11604, -11792, -11980, -12166, -12353,
		-12539, -12724, -12909, -13094, -13278, -13462, -13645, -13827,
		-14009, -14191, -14372, -14552, -14732, -14911, -15090, -15268,
		-15446, -15623, -15799, -15975, -16150, -16325, -16499, -16672,
		-16845, -17017, -17189, -17360, -17530, -17699, -17868, -18036,
		-18204, -18371, -18537, -18702, -18867, -19031, -19194, -19357,
		-19519, -19680, -19840, -20000, -20159, -20317, -20474, -20631,
		-20787, -20942, -21096, -21249, -21402, -21554, -21705, -21855,
		-22004, -22153, -22301, -22448, -22594, -22739, -22883, -23027,
		-23169, -23311, -23452, -23592, -23731, -23869, -24006, -24143,
		-24278, -24413, -24546, -24679, -24811, -24942, -25072, -25201,
		-25329, -25456, -25582, -25707, -25831, -25954, -26077, -26198,
		-26318, -26437, -26556, -26673, -26789, -26905, -27019, -27132,
		-27244, -27355, -27466, -27575, -27683, -27790, -27896, -28001,
		-28105, -28208, -28309, -28410, -28510, -28608, -28706, -28802,
		-28897, -28992, -29085, -29177, -29268, -29358, -29446, -29534,
		-29621, -29706, -29790, -29873, -29955, -30036, -30116, -30195,
		-30272, -30349, -30424, -30498, -30571, -30643, -30713, -30783,
		-30851, -30918, -30984, -31049, -31113, -31175, -31236, -31297,
		-31356, -31413, -31470, -31525, -31580, -31633, -31684, -31735,
		-31785, -31833, -31880, -31926, -31970, -32014, -32056, -32097,
		-32137, -32176, -32213, -32249, -32284, -32318, -32350, -32382,
		-32412, -32441, -32468, -32495, -32520, -32544, -32567, -32588,
		-32609, -32628, -32646, -32662, -32678, -32692, -32705, -32717,
		-32727, -32736, -32744, -32751, -32757, -32761, -32764, -32766,
		-32767, -32766, -32764, -32761, -32757, -32751, -32744, -32736,
		-32727, -32717, -32705, -32692, -32678, -32662, -32646, -32628,
		-32609, -32588, -32567, -32544, -32520, -32495, -32468, -32441,
		-32412, -32382, -32350, -32318, -32284, -32249, -32213, -32176,
		-32137, -32097, -32056, -32014, -31970, -31926, -31880, -31833,
		-31785, -31735, -31684, -31633, -31580, -31525, -31470, -31413,
		-31356, -31297, -31236, -31175, -31113, -31049, -30984, -30918,
		-30851, -30783, -30713, -30643, -30571, -30498, -30424, -30349,
		-30272, -30195, -30116, -30036, -29955, -29873, -29790, -29706,
		-29621, -29534, -29446, -29358, -29268, -29177, -29085, -28992,
		-28897, -28802, -28706, -28608, -28510, -28410, -28309, -28208,
		-28105, -28001, -27896, -27790, -27683, -27575, -27466, -27355,
		-27244, -27132, -27019, -26905, -26789, -26673, -26556, -26437,
		-26318, -26198, -26077, -25954, -25831, -25707, -25582, -25456,
		-25329, -25201, -25072, -24942, -24811, -24679, -24546, -24413,
		-24278, -24143, -24006, -23869, -23731, -23592, -23452, -23311,
		-23169, -23027, -22883, -22739, -22594, -22448, -22301, -22153,
		-22004, -21855, -21705, -21554, -21402, -21249, -21096, -20942,
		-20787, -20631, -20474, -20317, -20159, -20000, -19840, -19680,
		-19519, -19357, -19194, -19031, -18867, -18702, -18537, -18371,
		-18204, -18036, -17868, -17699, -17530, -17360, -17189, -17017,
		-16845, -16672, -16499, -16325, -16150, -15975, -15799, -15623,
		-15446, -15268, -15090, -14911, -14732, -14552, -14372, -14191,
		-14009, -13827, -13645, -13462, -13278, -13094, -12909, -12724,
		-12539, -12353, -12166, -11980, -11792, -11604, -11416, -11227,
		-11038, -10849, -10659, -10469, -10278, -10087,  -9895,  -9703,
		-9511,  -9319,  -9126,  -8932,  -8739,  -8545,  -8351,  -8156,
		-7961,  -7766,  -7571,  -7375,  -7179,  -6982,  -6786,  -6589,
		-6392,  -6195,  -5997,  -5799,  -5601,  -5403,  -5205,  -5006,
		-4807,  -4608,  -4409,  -4210,  -4011,  -3811,  -3611,  -3411,
		-3211,  -3011,  -2811,  -2610,  -2410,  -2209,  -2009,  -1808,
		-1607,  -1406,  -1206,  -1005,   -804,   -603,   -402,   -201,
};

the fft.h looks like this,

Code: 
#ifndef FFT_H
#define FFT_H
#define FIX_MPY(DEST,A,B)       DEST = ((long)(A) * (long)(B))>>15

#define N_WAVE          1024    /* dimension of Sinewave[] */
#define LOG2_N_WAVE     10      /* log2(N_WAVE) */

#ifndef fixed
#define fixed short
#endif

extern fixed Sinewave[N_WAVE];

//function prototypes
fixed fix_mpy(fixed a, fixed b);
int fix_fft(fixed *fr, fixed *fi, int m, int inverse);
#endif

The expected output is:
Input Data
0: 1000, 0
1: 707, 0
2: 0, 0
3: -707, 0
4: -1000, 0
5: -707, 0
6: 0, 0
7: 707, 0

FFT
0: -1, 0
1: 499, -1
2: 0, 0
3: 1, -1
4: 1, 0
5: 1, 1
6: 0, 0
7: 499, 1

IFFT
0: 995, -2
1: 705, 0
2: -1, 0
3: -704, 1
4: -997, 2
5: -707, 0
6: -1, 0
7: 702, -1


But I am getting the following,
Input Data
0: 1000, 0
1: 707, 0
2: 0, 0
3: -707, 0
4: -1000, 0
5: -707, 0
6: 0, 0
7: 707, 0

FFT
0: 0, 0
1: 499, -1
2: 0, 0
3: 0, 0
4: 0, 0
5: 0, 0
6: 0, 0
7: 498, 0

IFFT
0: 997, -1
1: 705, -2
2: 1, 0
3: -703, 0
4: -997, 0
5: -705, 1
6: -1, 0
7: 702, 0


Thanks in advance!!
 

Hi Amit,
I already looked at this solution before posting. I need some debugging in for my code. I do not know where is the problem. The twiddle factor computation, the real, imaginary parts all seems ok. Can you help? Thanks in advance!
 

Hi,

Why do you expect:

*****
FFT
0: -1, 0
1: 499, -1
2: 0, 0
3: 1, -1
4: 1, 0
5: 1, 1
6: 0, 0
7: 499, 1
*****

I`d expect:
*****
FFT
0: 0, 0
1: 500, 0
2: 0, 0
3: 0, 0
4: 0, 0
5: 0, 0
6: 0, 0
7: 500, 0
****

Beause the input is pure sine without DC and without other frequency components.
There may be calculation errors.

So I´d say: what you get is more near what I expect.


Klaus
 

Hi Klaus,
now i realize that i can have some room for value fluctuation. thanks for pointing it out.

--
Tarjina
 

Status
Not open for further replies.

Similar threads

Menu
Log in
Register
Cookies are required to use this site. You must accept them to continue using the site. Learn more…