1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72
| void fft_float ( unsigned NumSamples,
int InverseTransform,
float *RealIn,
float *ImagIn,
float *RealOut,
float *ImagOut )
{
unsigned NumBits; /* Number of bits needed to store indices */
unsigned i, j, k, n;
unsigned BlockSize, BlockEnd;
float angle_numerator = 2.0F * 3.141592F;
float delta_angle;
float alpha, beta; /* used in recurrence relation */
float delta_ar;
float tr, ti; /* temp real, temp imaginary */
float ar, ai; /* angle vector real, angle vector imaginary */
if (InverseTransform)
angle_numerator = -angle_numerator;
NumBits = NumberOfBitsNeeded ( NumSamples );
for ( i=0; i < NumSamples; i++ )
{
j = ReverseBits ( i, NumBits );
RealOut[j] = RealIn[i];
ImagOut[j] = ImagIn ? ImagIn[i] : 0.0F;
}
BlockEnd = 1;
for ( BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1 )
{
delta_angle = angle_numerator / (float)BlockSize;
alpha = (float)sin(0.5 * delta_angle);
alpha = 2.0F * alpha * alpha;
beta = (float)sin (delta_angle);
for ( i=0; i < NumSamples; i += BlockSize )
{
ar = 1.0; /* cos(0) */
ai = 0.0; /* sin(0) */
for ( j=i, n=0; n < BlockEnd; j++, n++ )
{
k = j + BlockEnd;
tr = ar*RealOut[k] - ai*ImagOut[k];
ti = ar*ImagOut[k] + ai*RealOut[k];
RealOut[k] = RealOut[j] - tr;
ImagOut[k] = ImagOut[j] - ti;
RealOut[j] += tr;
ImagOut[j] += ti;
delta_ar = alpha*ar + beta*ai;
ai -= (alpha*ai - beta*ar);
ar -= delta_ar;
}
}
BlockEnd = BlockSize;
}
if ( InverseTransform )
{
float denom = (float)NumSamples;
for ( i=0; i < NumSamples; i++ )
{
RealOut[i] /= denom;
ImagOut[i] /= denom;
}
}
} |
Partager