Můžeš zkusit kód co jsem napsal kdysi, jenom si změň floaty třeba na 16b inty (akorát nemám už tušení co by mělo být znaménkové):
/*
n. iterace
x0 ------+----+-(/2)------+-----+--(/2)- X0
\ / \ /
/ \ \ /
x2 -(W0)+(-1)+--(/2)------+--X--+--(/2)- X1
\/ \/
/\ /\
x1 ------+----+-(/2)-(W0)-+(-1)-+--(/2)- X2
\ / / \
/ \ / \
x3 -(W0)+(-1)+--(/2)-(W1)-+(-1)-+--(/2)- X3
1. a 2. iterace
x0 -----+----+-------+-----+--(/2)- X0
\ / \ /
/ \ \ /
x2 ----+(-1)+--------+--X--+--(/2)- X1
\/ \/
/\ /\
x1 -----+----+-------+(-1)-+--(/2)- X2
\ / / \
/ \ / \
x3 ----+(-1)+--(-j)--+(-1)-+--(/2)- X3
*/
void fft(float* x, BYTE ni, float* r, float* i)
{
DWORD b = 0;
for(DWORD a = 0; a < (1 << ni); a++) // prehazet vstupni data do R
{
r[a] = x[b];
for(BYTE c = ni - 1; c != 0xFF; c--) // b je pocitadlo s obracenym poradim bitu
{
b = b ^ (1 << c);
if((1 << c) & b) break;
}
}
for(DWORD j = 0; j < (1 << ni); j = j + 4) // 1. a 2. iterace (W = 1 nebo 0, takze neni treba pouzivat tabulku)
{
float a = r[j];
float b = r[j + 1];
r[j] = (a + b + r[j + 2] + r[j + 3]) / 2;
i[j] = 0;
r[j+1] = (a - b) / 2;
i[j+1] = (-r[j + 2] + r[j + 3]) /2;
r[j+2] = (a + b - r[j + 2] - r[j + 3]) / 2;
i[j+2] = 0;
r[j+3] = r[j + 1];
i[j+3] = -i[j + 1];
}
for(BYTE i = 3; i <= ni; i++) // n. iterace
{
for(DWORD j = 0; j < (1 << ni); j = j + (1 << i)) // pro kazdy "svazek" DFT2
{
DWORD _i = 1 << (i - 1);
for(DWORD k = 0; k < (1 << (i - 1)); k++) // pro kazde 2 DFT2 ve svazku
{
float _r = r[j + k];
float _i = i[j + k];
float _r2 = r[j + _i];
float _i2 = i[j + _i];
float _wr = cos[k * (1 << (ni - i)]; // predvypocitane sin/cos tabulky o velikosti 1 << ni
float _wi = msin[k * (1 << (ni - i)];
r[j + k] = (_r + (_wr * _r2) - (_wi * _i2)) / 2;
i[j + k] = (_i + (_wi * _r2) + (_wr * _i2)) / 2;
r[j + _i] = (_r - (_wr * _r2) + (_wi * _i2)) / 2;
i[j + _i] = (_i - (_wi * _r2) - (_wr * _i2)) / 2;
_i ++;
}
}
}
}
Abys dostal amplitudu, musíš si to ještě přepočítat na polární souřadnice a taky osa frekvence co ti vyje je lineární, jak už psal Vroutek.