# Looking for a simple core FFT algorithm code

Status
Not open for further replies.

#### fala

##### Full Member level 5 Hello, I'm developing an FFT algorithm for impedance calculations. I searched the net and found very long and sophisticated codes and libraries for that like https://www.fftw.org I don't want to use pre-written libraries. I just want to know the very basic idea so I can find out what differences with DFT it has so hopefully I can write an algorithm in C with a few dozen lines.
Can someone give me a hint what keyword or what site I should look for the core idea of FFT?

thanks alot.

#### Sal

##### Full Member level 4 ### fala

Points: 2

#### cydi

##### Full Member level 3 simple fft

Check out DSP.com, and get the e-book Scientists and engineers guide to DSP. Yuo will get pseudo-code in it. Its pretty simple.

### fala

Points: 2

#### amit_8561

##### Member level 2 simple FFT

you can search for some help on dspguide.com

#### echo47 simple FFT

Here's a simple (low performance) 16-bit complex FFT. I forget where I found it, and I may have rewritten some of it. You can change NFFT to some other power of two. It assumes 32-bit int and 16-bit short. The output array overwrites the input array. The array's real and imaginary values are interleaved i0, j0, i1, j1, i2, j2, ...

Code:
#include <math.h>
#define PI 3.141592653589793238462643
#define NFFT 256

void fft(short *x)
{
static short once=1, sint[3 * NFFT / 4];
int n1, n2, ie, i, j, k, xt, yt;
short *xi, *xl, *sinus;

if (once)
{
for (i=0; i < 3*NFFT/4; i++)    /* sine table */
sint[i] = 0x7FFF * sin(2 * PI * i / NFFT);
once = 0;
}

for (n2=NFFT, ie=1; ie<NFFT; ie=ie<<1)
{
n1 = n2;
n2 = n2 >> 1;
sinus = sint;
for (j=0; j < n2; j++)
{
for (i=j<<1; i < 2*NFFT; i+=n1<<1)
{
xi = x + i;
xl = xi + n1;
xt = (xi - xl) >> 1;
yt = (xi - xl) >> 1;
*(int*)xi = (((xi + xl) << 15) & 0xFFFF0000) | (unsigned short)((xi + xl) >> 1);
*(int*)xl = (((xt * sinus[NFFT/4] + yt * sinus) << 1) & 0xFFFF0000) | (unsigned short)((yt * sinus[NFFT/4] - xt * sinus) >> 15);
}
sinus = sinus + ie;
}
}

for (j=i=0; i < 2*NFFT; i+=2)   /* bit reverse */
{
if (i < j)
{
ie = *(int*)(x+i);
*(int*)(x+i) = *(int*)(x+j);
*(int*)(x+j) = ie;
}
k = NFFT;
while (j >= k)
{
j -= k;
k = k >> 1;
}
j += k;
}
}

### fala

Points: 2

#### cedance Re: simple FFT

here is another 1D fft, i had to write for my image processing project. as for excellent explanations of algorithms (many), both in DIT and DIF, check out Inside fft blackbox book from CRC press.

Code:
void hiwi::fft1(Vec_IO_DP &inpR, Vec_IO_DP &inpI, int Inv)

{

int NumberOfProblems = 1;
int ProblemSize = inpR.size();
int N = inpR.size();
int Distance = 1;
bool NotSwitchInput = 1;
int i;
int JFirst = 0;
int J, K, JTwiddle;
double WR, WI;
Vec_DP wR(N), wI(N), aR(N), aI(N), bR(N), bI(N);

for(i=0; i<N; i++) {
wR[i] = 0;
wI[i] = 0;
aR[i] = inpR[i];
aI[i] = inpI[i];
bR[i] = 0;
bI[i] = 0;
}
for(i=0; i<N; i++) {
wR[i] = cos(2* pi * i / N);
if(Inv == 1)
wI[i] = -sin(2* pi * i / N);
else
wI[i] = +sin(2* pi * i / N);
}
while(ProblemSize > 1) {
if(NotSwitchInput) {
for(JFirst = 0; JFirst < NumberOfProblems; JFirst++) {
J = JFirst;
JTwiddle = 0;
K = JFirst;
while( J < N - 1 ) {
WR = wR[JTwiddle];
WI = wI[JTwiddle];
bR[J] = aR[K] + aR[K + N/2];
bI[J] = aI[K] + aI[K + N/2];
bR[J + Distance] = (aR[K] - aR[K + N/2]) * WR - (aI[K] - aI[K + N/2]) * WI;
bI[J + Distance] = (aR[K] - aR[K + N/2]) * WI + (aI[K] - aI[K + N/2]) * WR;
JTwiddle += NumberOfProblems;
J += 2 * NumberOfProblems;
K += NumberOfProblems;
}
}
NotSwitchInput = 0;
}
else {
for(JFirst = 0; JFirst < NumberOfProblems;JFirst++) {
J = JFirst;
JTwiddle = 0;
K = JFirst;
while(J < N - 1) {
WR = wR[JTwiddle];
WI = wI[JTwiddle];
aR[J] = bR[K] + bR[K + N/2];
aI[J] = bI[K] + bI[K + N/2];
aR[J + Distance] = (bR[K] - bR[K+N/2]) * WR - (bI[K] - bI[K + N/2]) * WI;
aI[J + Distance] = (bR[K] - bR[K+N/2]) * WI + (bI[K] - bI[K + N/2]) * WR;
JTwiddle += NumberOfProblems;
J += 2 * NumberOfProblems;
K += NumberOfProblems;
}
}
NotSwitchInput = 1;
}
NumberOfProblems *= 2;
ProblemSize /= 2;
Distance *= 2;
}
double factor = log((double) N) / log((double) 2);
if ( fmod(factor,(double) 2) == 0) {
for(i = 0; i<N; i++) {
inpR[i] = aR[i];
inpI[i] = aI[i];
}
}
else {
for(i = 0; i<N; i++) {
inpR[i] = bR[i];
inpI[i] = bI[i];
}
}

}

the first and the second parameters are the input vector real and imaginary part and the third parameter tells whether to calculate fft or ifft. a 1 calculates fft and 0 or any other number calculates ifft.

2D fft should be straight forward from this as u can take the 1D fft of the rows and from the resulting matrix, take 1D fft of columns. hope this might be of use to some one. however, the function input specified here uses matrix templates which i wrote. so, specify the i/p and vector/matrix declaration yourself, as to how u usually do it.

regards,
cedance.

Points: 2