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

tarjina · Mar 30, 2016

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!!

amit.kumar11 · Mar 30, 2016

Hi,

There is a post ..Plz check

https://www.edaboard.com/threads/82466/

Amit

tarjina · Mar 30, 2016

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!

KlausST · Mar 30, 2016

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

tarjina · Mar 30, 2016

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

--
Tarjina

Welcome to EDAboard.com

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

tarjina

Junior Member level 3

amit.kumar11

Member level 5

tarjina

Junior Member level 3

KlausST

Advanced Member level 7

tarjina

Junior Member level 3

Similar threads

Part and Inventory Search

Welcome to EDABoard.com

Sponsor

Connect with us

Online statistics

Forum statistics