97 lines
2.3 KiB
C
97 lines
2.3 KiB
C
|
||
#include "config.h"
|
||
|
||
#include "alcomplex.h"
|
||
#include "math_defs.h"
|
||
|
||
|
||
extern inline ALcomplex complex_add(ALcomplex a, ALcomplex b);
|
||
extern inline ALcomplex complex_sub(ALcomplex a, ALcomplex b);
|
||
extern inline ALcomplex complex_mult(ALcomplex a, ALcomplex b);
|
||
|
||
void complex_fft(ALcomplex *FFTBuffer, ALsizei FFTSize, ALdouble Sign)
|
||
{
|
||
ALsizei i, j, k, mask, step, step2;
|
||
ALcomplex temp, u, w;
|
||
ALdouble arg;
|
||
|
||
/* Bit-reversal permutation applied to a sequence of FFTSize items */
|
||
for(i = 1;i < FFTSize-1;i++)
|
||
{
|
||
for(mask = 0x1, j = 0;mask < FFTSize;mask <<= 1)
|
||
{
|
||
if((i&mask) != 0)
|
||
j++;
|
||
j <<= 1;
|
||
}
|
||
j >>= 1;
|
||
|
||
if(i < j)
|
||
{
|
||
temp = FFTBuffer[i];
|
||
FFTBuffer[i] = FFTBuffer[j];
|
||
FFTBuffer[j] = temp;
|
||
}
|
||
}
|
||
|
||
/* Iterative form of DanielsonLanczos lemma */
|
||
for(i = 1, step = 2;i < FFTSize;i<<=1, step<<=1)
|
||
{
|
||
step2 = step >> 1;
|
||
arg = M_PI / step2;
|
||
|
||
w.Real = cos(arg);
|
||
w.Imag = sin(arg) * Sign;
|
||
|
||
u.Real = 1.0;
|
||
u.Imag = 0.0;
|
||
|
||
for(j = 0;j < step2;j++)
|
||
{
|
||
for(k = j;k < FFTSize;k+=step)
|
||
{
|
||
temp = complex_mult(FFTBuffer[k+step2], u);
|
||
FFTBuffer[k+step2] = complex_sub(FFTBuffer[k], temp);
|
||
FFTBuffer[k] = complex_add(FFTBuffer[k], temp);
|
||
}
|
||
|
||
u = complex_mult(u, w);
|
||
}
|
||
}
|
||
}
|
||
|
||
/*Discrete Hilbert Transform (analytic signal form)*/
|
||
void hilbert(ALsizei size, ALcomplex *InOutBuffer )
|
||
{
|
||
ALsizei k;
|
||
const ALdouble inverse_size = 1.0/(ALfloat)size;
|
||
|
||
for ( k = 0; k < size;k++ )
|
||
InOutBuffer[k].Imag = 0.0;
|
||
|
||
complex_fft( InOutBuffer, size, 1.0 );
|
||
|
||
for( k = 0; k < size; k++ )
|
||
{
|
||
if( k == 0 || k == size/2 )
|
||
{
|
||
InOutBuffer[k].Real *= inverse_size;
|
||
InOutBuffer[k].Imag *= inverse_size;
|
||
}
|
||
|
||
else if ( k >=1 && k < size/2 )
|
||
{
|
||
InOutBuffer[k].Real *= 2.0*inverse_size;
|
||
InOutBuffer[k].Imag *= 2.0*inverse_size;
|
||
}
|
||
|
||
else
|
||
{
|
||
InOutBuffer[k].Real = 0.0;
|
||
InOutBuffer[k].Imag = 0.0;
|
||
}
|
||
}
|
||
|
||
complex_fft( InOutBuffer, size,-1.0 );
|
||
}
|