// БПФ с вычислением фаз
#include < math.h >

#define _N_SAMPLES_ (512) // 32,64,128,256,512,1024 points
#define _F_SAMPLES_ (100.0) // Частота тестового сигнала, Hz
#define _A_SAMPLES_ (100.0) // Амплитуда тестового сигнала

#if !defined M_PI
#define M_PI (3.14159265359)
#endif

float x[ _N_SAMPLES_ ], y[ _N_SAMPLES_ ], F, A, Re, Im, Fi;

//---------------------------------------------------------------
// I=1 --> БПФ
// I=-1 ---> Oбратное БПФ
void BPF( float *x, float *y, int N, signed char I )
{
volatile float c, s, t1, t2, t3, t4, u1, u2, u3;
volatile int i, j, p, l, L, M, M1, K;

L = N;
M = (N / 2);
M1 = (N - 1);

while( L >= 2 )
{
l = (L / 2);
u1 = 1.0000;
u2 = 0.000;
t1 = M_PI / (float)l;
c = cos( t1 );
s = (-1.000) * I * sin( t1 );

for( j = 0; j < l; j++ )
{
for( i = j; i < N; i += L )
{
p = i + l;
t1 = *(x+i) + *(x+p);
t2 = *(y+i) + *(y+p);
t3 = *(x+i) - *(x+p);
t4 = *(y+i) - *(y+p);
*(x+p) = (t3*u1) - (t4*u2);
*(y+p) = (t4*u1) + (t3*u2);
*(x+i) = t1;
*(y+i) = t2;
}/*for*/
u3 = (u1*c) - (u2*s);
u2 = (u2*c) + (u1*s);
u1 = u3;
}/*for*/
L /= 2;
}/*while*/

j = 0;

for( i = 0; i < M1; i++ )
{
if( i > j )
{
t1 = *(x+j);
t2 = *(y+j);
*(x+j) = *(x+i);
*(y+j) = *(y+i);
*(x+i) = t1;
*(y+i) = t2;
}/*if*/

K = M;

while(j >= K)
{
j -= K;
K /= 2;
}/*while*/
j += K;
}/*for*/

}/*BPF*/
//----------------------------------------------------
/* Тест. Моделирование синусоидального тестового сигнала 1 */
/* Один период синуса за N измерений */
/* P - Массив выборок размерностью N */
/* A - Амплитуда (Частота _F_SAMPLES_ не имеет значения)*/
void Sinsignal_For_Test_One_Period( float *P, float A, int N )
{
 int i;
 float r, re, re1, im, im1;

for( i = 0; i < N; i++ )
{
*(P+i) = A * cos( 2.00 * M_PI *(float)i / (float)N );
}/*for*/
}
//----------------------------------------------------
/* Моделирование синусоидального тестового сигнала 2 */
/* F периодов синуса за N измерений */
/* P - Массив выборок размерностью N */
/* F - Частота, Гц */
/* A - Амплитуда */
void Sinsignal_For_Test( float *P, float F, float A, int N )
{
int i;
float r, re, re1, im, im1;

re = cos((2.00 * M_PI * F ) / (float)N);
im = sin((2.00 * M_PI * F ) / (float)N);

re1 = A;
im1 = 0.00;

for( i = 0; i < N; i++ )
{
*(P+i) = re1;
r = re1;
re1 = (r * re) - (im1 * im);
im1 = (im1 * re) + (r * im);

}/*for*/
}
//----------------------------------------------------
__C_task main( void )
{
unsigned int j, N;

N = _N_SAMPLES_; // Points
F = _F_SAMPLES_; // Hz
A = _A_SAMPLES_; // Amplituda

////////////////////////////////
// Задаем входной сигнал
////////////////////////////////

// Пример 1. F периодов синуса за N выборок
//Sinsignal_For_Test( x, F, A, N ); // Генерация sin

// Пример 2. Один период синуса за N выборок
Sinsignal_For_Test_One_Period( x, A, N ); // Генерация sin

// Пример 3. Импульсный сигнал
// !!!!!! Сделать _N_SAMPLES_ = 32.
// _А_SAMPLES_ _F_SAMPLES_ не имеют значения
// for( j=0; j<8; j++ )
// {
// x[j] = 1.000;
// }/*for*/
// Остальные х[9...31] нули. Все y[0...31] нули

////////////////////////////////
// Вызов БПФ
////////////////////////////////
BPF( x, y, N, 1);
////////////////////////////////
for( j = 0; j < N/2; j++ )
{
/////////////// Расчет амплитуд
Re = *(x + j);
Im = *(y + j);
///////////////////////////////////////
A = sqrt( (Re * Re) + (Im * Im) );
if( j )
{
A = (A * 2.000);
}/*if*/
A = ( A / (float)N);
////////////////////////////////////////////////////
/////////////// Расчет фаз
Fi = 0.000;
if( (A < (-0.001)) || (A > 0.001 ) )
{
if( (Re > (-0.000000001)) && (Re < 0.000000001 ) )
{
if( Im > 0.000000001 ) Fi = 1.5708;
if( Im < (-0.000000001) ) Fi = (-1.5708);
}/*if*/
else
{
Fi = atan2( Im , Re );
Fi = ((180.000 * Fi) / M_PI);
}/*else*/
}/*if*/
////////////////////////////////////////////////////
// Вывод на экран
///////////////////////////////////////////////////
printf(" A[%d] = %.3f, Fi = %.2f \n", j, A, Fi );
///////////////////////////////////////////////////
}/*for*/
}/*main*/