fft code.txt

int reverse(int n, int v)
/* Odwrenie kolejnoｶci v najmｳodszych bit liczby n - funkcja pomocnicza */
{
  int bit, r = 0;
  for (bit=0; bit<v; bit++) {
    r = (r << 1) | (n & 1);
    n >>= 1;
  }
  return r;
}

/******************************************************************************/

float fft(float *x, float *y, int v, float *SinTab, int M)

/*
     Parametry wejｶciowe:
       x       - N prek rzeczywistego sygnaｳu znormalizowanego i oknowanego
       v       - log2(N)
       SinTab  - tablica wartoｶci funkcji sin(ﾟ) dla ﾟ z przedziaｳu [0,2pi)
       M       - liczba prek funkcji sinus w okresie (rozmiar tablicy SinTab)
     Parametry wyjsciowe:
       y       - N prek logarytmicznego widma FFT sygnaｳu wyraｿonego w [daB],
		 dla cz黌totliwoｶci znormalizowanych fT z przedziaｳu [0,1)
     Zwracana wartoｶ・
       moc sygnaｳu x wyraｿona w [daB] (liczona z twierdzenia Parsevala)
*/

{
  int           p, q, r, grp, bfy;
  register int  lc1, lc2, lc3;
  float         Sin, Cos, xp, yp, xq, yq, sum;

  /* Wyznaczenie liczby punkt FFT */

  #define N (1 << v)

  /* Wyzerowanie cz・ci urojonej (y) */

  for (lc1=0; lc1<N; lc1++)
    y[lc1] = 0.0;

  /* Szybkie przeksztaｳcenie Fouriera wg Cooley'a i Tukey'a
     (obliczenia z podstawieniem, prki widma w odwronym porzｱdku bit),
     por. [Oppenheim, Schafer], wz (6.16) i rys.6-12 */

  grp = 1;
  bfy = N/2;
  for (lc1=0; lc1<v; lc1++) { /* kolejne przebiegi oblicze・z podstawianiem */
    p = 0;
    q = bfy;
    r = 0;
    for (lc2=0; lc2<grp; lc2++) { /* grupy motylk w ramach danego przebiegu */
      Sin = SinTab[r];
      Cos = SinTab[r+M/4];
      r = reverse(reverse(r,v-1)+1,v-1);
      for (lc3=0; lc3<bfy; lc3++) { /* motylki w ramach danej grupy */
	xq = x[q];
	yq = y[q];
	xp = x[p] + xq*Cos + yq*Sin; x[q] = 2*x[p] - xp; x[p] = xp;
	yp = y[p] + yq*Cos - xq*Sin; y[q] = 2*y[p] - yp; y[p] = yp;
	p++;
	q++;
      }
      p += bfy;
      q += bfy;
    }
    grp *= 2;
    bfy /= 2;
  }

  /* Wyznaczenie kwadratu widma amplitudowego */

  for (lc1=0; lc1<N; lc1++)
    x[lc1] = x[lc1]*x[lc1] + y[lc1]*y[lc1];

  /* Tasowanie prek widma przywracajｱce ich naturalne uporzｱdkowanie */

  for (lc1=0; lc1<N; lc1++)
    y[reverse(lc1,v)] = x[lc1];

  /* Logarytmowanie widma [daB] */

  for (lc1=0; lc1<N; lc1++)
    y[lc1] = decabel(y[lc1]);

  /* Obliczenie logarytmu mocy [daB] na podstawie twierdzenia Parsevala */

  sum = 0.0;
  for (lc1=0; lc1<N; lc1++)
    sum += x[lc1];

  return decabel(sum);
}