/* FFT routine             */
/* (c) 2025 cepstrum.co.jp */

/**                                            ***/
/*** you can use better FFT routine (FFTW etc) ***/
/***                                           ***/

float cos_table[FFTLEN/2];
float sin_table[FFTLEN/2];
int   rev_table1[FFTLEN/2];
int   rev_table2[FFTLEN/2];

//----------------------------------------------

int ilog2(int x) {
  int i, j;

  i=x;
  j=0;
  do {
    j=j+1;
    i=i/2;
  } while (i!=1);
  return j;
}

//----------------------------------------------

void make_cos_sin_table(float cos_table[], float sin_table[], int len) {
  float cosine;
  int   i;

  cos_table[0]=1.0;
  cos_table[len/4]=0.0;
  sin_table[0]=0.0;
  sin_table[len/4]=1.0;
  for (i=1; i<len/4; i=i+1) {
    cosine=cos(2.0*PI/len*i);
    cos_table[i]=cosine;
    cos_table[len/2-i]=-cosine;
    sin_table[len/4-i]=cosine;
    sin_table[len/4+i]=cosine;
  }
}

//----------------------------------------------

void make_bit_rev_table(int x[], int y[], int *table_len, int len) {
  int   i;
  float remainder, sum, weight, subtrahend;

  *table_len=0;
  for (i=1; i<=len-2; i=i+1) {
    remainder=i;
    sum=0;
    weight=1;
    subtrahend=len/2;
    while (0<remainder) {
      if (subtrahend<=remainder) {
        sum=sum+weight;
        remainder=remainder-subtrahend;
      }
      weight=2*weight;
      subtrahend=subtrahend/2;
    }
    if (i<sum) {
      x[*table_len]=i;
      y[*table_len]=sum;
      *table_len=*table_len+1;
    }
  }
}

//----------------------------------------------

void bit_reverse(float x[], int table1[], int table2[], int len) {
  float tmp;
  int   i;

  for (i=0; i<len; i=i+1) {
    tmp=x[table1[i]];
    x[table1[i]]=x[table2[i]];
    x[table2[i]]=tmp;
  }
}

//----------------------------------------------

void fft(float ar[], float ai[], int fft_flag, int init) {
  int    stage, i, j;
  int    ix1, ix2, ix_step, ix_difference;
  int    angle, angle_step;
  float  co, si;
  float  wi, wr;
  static int   rev_table_len;

  if (init!=0) {
    make_cos_sin_table(cos_table, sin_table, FFTLEN);
    make_bit_rev_table(rev_table1, rev_table2, &rev_table_len, FFTLEN);
  } else {
    bit_reverse(ar, rev_table1, rev_table2, rev_table_len);
    bit_reverse(ai, rev_table1, rev_table2, rev_table_len);
    angle_step=FFTLEN/2;
    ix_difference=1;
    ix_step=2;
    for (stage=1; stage<=ilog2(FFTLEN); stage=stage+1) {
      angle=0;
      for (i=0; i<ix_difference; i=i+1) {
        ix1=i;
        co=cos_table[angle];
        si=sin_table[angle]*fft_flag;
        for (j=1; j<=angle_step; j=j+1) {
          ix2=ix1+ix_difference;
          wr=ar[ix2]*co-ai[ix2]*si;
          wi=ar[ix2]*si+ai[ix2]*co;
          ar[ix2]=ar[ix1]-wr;
          ai[ix2]=ai[ix1]-wi;
          ar[ix1]=ar[ix1]+wr;
          ai[ix1]=ai[ix1]+wi;
          ix1=ix1+ix_step;
        }
        angle=angle+angle_step;
      }
      angle_step=angle_step/2;
      ix_difference=ix_difference*2;
      ix_step=ix_step*2;
    }
    if (0<fft_flag) {
      for (i=0; i<FFTLEN; i=i+1) {
        ar[i]=ar[i]/FFTLEN;
        ai[i]=ai[i]/FFTLEN;
      }
    }
  }
}

//----------------------------------------------