ref: 96270757ecc40d555aa27236c970dda9b4ec1975
dir: /lib/fft/fftlib.c/
/* Public domain FFT implementation in C by John Green */ #include <stdint.h> #include <stdlib.h> #include <math.h> #include "base.h" #ifndef M_PI #define M_PI 3.14159265358979323846 #endif #define LOG2(n) (log(n) / log(2)) #define POW2(m) ((uint32_t) 1 << (m)) /* integer power of 2 for m<32 */ /* fft's with M bigger than this bust primary cache */ #define MCACHE (11 - (sizeof(SPFLOAT) / 8)) /* some math constants to 40 decimal places */ #define MYPI 3.141592653589793238462643383279502884197 /* pi */ #define MYROOT2 1.414213562373095048801688724209698078569 /* sqrt(2) */ #define MYCOSPID8 0.9238795325112867561281831893967882868224 /* cos(pi/8) */ #define MYSINPID8 0.3826834323650897717284599840303988667613 /* sin(pi/8) */ /***************************************************** * routines to initialize tables used by fft routines * *****************************************************/ static void fftCosInit(int M, SPFLOAT *Utbl) { /* Compute Utbl, the cosine table for ffts */ /* of size (pow(2,M)/4 +1) */ /* INPUTS */ /* M = log2 of fft size */ /* OUTPUTS */ /* *Utbl = cosine table */ unsigned int fftN = POW2(M); unsigned int i1; Utbl[0] = 1.0; for (i1 = 1; i1 < fftN/4; i1++) Utbl[i1] = cos((2.0 * M_PI * (SPFLOAT)i1) / (SPFLOAT)fftN); Utbl[fftN/4] = 0.0; } void fftBRInit(int M, int16_t *BRLow) { /* Compute BRLow, the bit reversed table for ffts */ /* of size pow(2,M/2 -1) */ /* INPUTS */ /* M = log2 of fft size */ /* OUTPUTS */ /* *BRLow = bit reversed counter table */ int Mroot_1 = M / 2 - 1; int Nroot_1 = POW2(Mroot_1); int i1; int bitsum; int bitmask; int bit; for (i1 = 0; i1 < Nroot_1; i1++) { bitsum = 0; bitmask = 1; for (bit = 1; bit <= Mroot_1; bitmask <<= 1, bit++) if (i1 & bitmask) bitsum = bitsum + (Nroot_1 >> bit); BRLow[i1] = bitsum; } } /***************** * parts of ffts1 * *****************/ static void bitrevR2(SPFLOAT *ioptr, int M, int16_t *BRLow) { /*** bit reverse and first radix 2 stage of forward or inverse fft ***/ SPFLOAT f0r; SPFLOAT f0i; SPFLOAT f1r; SPFLOAT f1i; SPFLOAT f2r; SPFLOAT f2i; SPFLOAT f3r; SPFLOAT f3i; SPFLOAT f4r; SPFLOAT f4i; SPFLOAT f5r; SPFLOAT f5i; SPFLOAT f6r; SPFLOAT f6i; SPFLOAT f7r; SPFLOAT f7i; SPFLOAT t0r; SPFLOAT t0i; SPFLOAT t1r; SPFLOAT t1i; SPFLOAT *p0r; SPFLOAT *p1r; SPFLOAT *IOP; SPFLOAT *iolimit; int Colstart; int iCol; unsigned int posA; unsigned int posAi; unsigned int posB; unsigned int posBi; const unsigned int Nrems2 = POW2((M + 3) / 2); const unsigned int Nroot_1_ColInc = POW2(M) - Nrems2; const unsigned int Nroot_1 = POW2(M / 2 - 1) - 1; const unsigned int ColstartShift = (M + 1) / 2 + 1; posA = POW2(M); /* 1/2 of POW2(M) complexes */ posAi = posA + 1; posB = posA + 2; posBi = posB + 1; iolimit = ioptr + Nrems2; for (; ioptr < iolimit; ioptr += POW2(M / 2 + 1)) { for (Colstart = Nroot_1; Colstart >= 0; Colstart--) { iCol = Nroot_1; p0r = ioptr + Nroot_1_ColInc + BRLow[Colstart] * 4; IOP = ioptr + (Colstart << ColstartShift); p1r = IOP + BRLow[iCol] * 4; f0r = *(p0r); f0i = *(p0r + 1); f1r = *(p0r + posA); f1i = *(p0r + posAi); for (; iCol > Colstart;) { f2r = *(p0r + 2); f2i = *(p0r + (2 + 1)); f3r = *(p0r + posB); f3i = *(p0r + posBi); f4r = *(p1r); f4i = *(p1r + 1); f5r = *(p1r + posA); f5i = *(p1r + posAi); f6r = *(p1r + 2); f6i = *(p1r + (2 + 1)); f7r = *(p1r + posB); f7i = *(p1r + posBi); t0r = f0r + f1r; t0i = f0i + f1i; f1r = f0r - f1r; f1i = f0i - f1i; t1r = f2r + f3r; t1i = f2i + f3i; f3r = f2r - f3r; f3i = f2i - f3i; f0r = f4r + f5r; f0i = f4i + f5i; f5r = f4r - f5r; f5i = f4i - f5i; f2r = f6r + f7r; f2i = f6i + f7i; f7r = f6r - f7r; f7i = f6i - f7i; *(p1r) = t0r; *(p1r + 1) = t0i; *(p1r + 2) = f1r; *(p1r + (2 + 1)) = f1i; *(p1r + posA) = t1r; *(p1r + posAi) = t1i; *(p1r + posB) = f3r; *(p1r + posBi) = f3i; *(p0r) = f0r; *(p0r + 1) = f0i; *(p0r + 2) = f5r; *(p0r + (2 + 1)) = f5i; *(p0r + posA) = f2r; *(p0r + posAi) = f2i; *(p0r + posB) = f7r; *(p0r + posBi) = f7i; p0r -= Nrems2; f0r = *(p0r); f0i = *(p0r + 1); f1r = *(p0r + posA); f1i = *(p0r + posAi); iCol -= 1; p1r = IOP + BRLow[iCol] * 4; } f2r = *(p0r + 2); f2i = *(p0r + (2 + 1)); f3r = *(p0r + posB); f3i = *(p0r + posBi); t0r = f0r + f1r; t0i = f0i + f1i; f1r = f0r - f1r; f1i = f0i - f1i; t1r = f2r + f3r; t1i = f2i + f3i; f3r = f2r - f3r; f3i = f2i - f3i; *(p0r) = t0r; *(p0r + 1) = t0i; *(p0r + 2) = f1r; *(p0r + (2 + 1)) = f1i; *(p0r + posA) = t1r; *(p0r + posAi) = t1i; *(p0r + posB) = f3r; *(p0r + posBi) = f3i; } } } static void fft2pt(SPFLOAT *ioptr) { /*** RADIX 2 fft ***/ SPFLOAT f0r, f0i, f1r, f1i; SPFLOAT t0r, t0i; /* bit reversed load */ f0r = ioptr[0]; f0i = ioptr[1]; f1r = ioptr[2]; f1i = ioptr[3]; /* Butterflys */ /* f0 - - t0 f1 - 1 - f1 */ t0r = f0r + f1r; t0i = f0i + f1i; f1r = f0r - f1r; f1i = f0i - f1i; /* store result */ ioptr[0] = t0r; ioptr[1] = t0i; ioptr[2] = f1r; ioptr[3] = f1i; } static void fft4pt(SPFLOAT *ioptr) { /*** RADIX 4 fft ***/ SPFLOAT f0r, f0i, f1r, f1i, f2r, f2i, f3r, f3i; SPFLOAT t0r, t0i, t1r, t1i; /* bit reversed load */ f0r = ioptr[0]; f0i = ioptr[1]; f1r = ioptr[4]; f1i = ioptr[5]; f2r = ioptr[2]; f2i = ioptr[3]; f3r = ioptr[6]; f3i = ioptr[7]; /* Butterflys */ /* f0 - - t0 - - f0 f1 - 1 - f1 - - f1 f2 - - f2 - 1 - f2 f3 - 1 - t1 - -i - f3 */ t0r = f0r + f1r; t0i = f0i + f1i; f1r = f0r - f1r; f1i = f0i - f1i; t1r = f2r - f3r; t1i = f2i - f3i; f2r = f2r + f3r; f2i = f2i + f3i; f0r = t0r + f2r; f0i = t0i + f2i; f2r = t0r - f2r; f2i = t0i - f2i; f3r = f1r - t1i; f3i = f1i + t1r; f1r = f1r + t1i; f1i = f1i - t1r; /* store result */ ioptr[0] = f0r; ioptr[1] = f0i; ioptr[2] = f1r; ioptr[3] = f1i; ioptr[4] = f2r; ioptr[5] = f2i; ioptr[6] = f3r; ioptr[7] = f3i; } static void fft8pt(SPFLOAT *ioptr) { /*** RADIX 8 fft ***/ SPFLOAT w0r = (SPFLOAT)(1.0 / MYROOT2); /* cos(pi/4) */ SPFLOAT f0r, f0i, f1r, f1i, f2r, f2i, f3r, f3i; SPFLOAT f4r, f4i, f5r, f5i, f6r, f6i, f7r, f7i; SPFLOAT t0r, t0i, t1r, t1i; const SPFLOAT Two = 2.0; /* bit reversed load */ f0r = ioptr[0]; f0i = ioptr[1]; f1r = ioptr[8]; f1i = ioptr[9]; f2r = ioptr[4]; f2i = ioptr[5]; f3r = ioptr[12]; f3i = ioptr[13]; f4r = ioptr[2]; f4i = ioptr[3]; f5r = ioptr[10]; f5i = ioptr[11]; f6r = ioptr[6]; f6i = ioptr[7]; f7r = ioptr[14]; f7i = ioptr[15]; /* Butterflys */ /* f0 - - t0 - - f0 - - f0 f1 - 1 - f1 - - f1 - - f1 f2 - - f2 - 1 - f2 - - f2 f3 - 1 - t1 - -i - f3 - - f3 f4 - - t0 - - f4 - 1 - t0 f5 - 1 - f5 - - f5 - w3 - f4 f6 - - f6 - 1 - f6 - -i - t1 f7 - 1 - t1 - -i - f7 - iw3- f6 */ t0r = f0r + f1r; t0i = f0i + f1i; f1r = f0r - f1r; f1i = f0i - f1i; t1r = f2r - f3r; t1i = f2i - f3i; f2r = f2r + f3r; f2i = f2i + f3i; f0r = t0r + f2r; f0i = t0i + f2i; f2r = t0r - f2r; f2i = t0i - f2i; f3r = f1r - t1i; f3i = f1i + t1r; f1r = f1r + t1i; f1i = f1i - t1r; t0r = f4r + f5r; t0i = f4i + f5i; f5r = f4r - f5r; f5i = f4i - f5i; t1r = f6r - f7r; t1i = f6i - f7i; f6r = f6r + f7r; f6i = f6i + f7i; f4r = t0r + f6r; f4i = t0i + f6i; f6r = t0r - f6r; f6i = t0i - f6i; f7r = f5r - t1i; f7i = f5i + t1r; f5r = f5r + t1i; f5i = f5i - t1r; t0r = f0r - f4r; t0i = f0i - f4i; f0r = f0r + f4r; f0i = f0i + f4i; t1r = f2r - f6i; t1i = f2i + f6r; f2r = f2r + f6i; f2i = f2i - f6r; f4r = f1r - f5r * w0r - f5i * w0r; f4i = f1i + f5r * w0r - f5i * w0r; f1r = f1r * Two - f4r; f1i = f1i * Two - f4i; f6r = f3r + f7r * w0r - f7i * w0r; f6i = f3i + f7r * w0r + f7i * w0r; f3r = f3r * Two - f6r; f3i = f3i * Two - f6i; /* store result */ ioptr[0] = f0r; ioptr[1] = f0i; ioptr[2] = f1r; ioptr[3] = f1i; ioptr[4] = f2r; ioptr[5] = f2i; ioptr[6] = f3r; ioptr[7] = f3i; ioptr[8] = t0r; ioptr[9] = t0i; ioptr[10] = f4r; ioptr[11] = f4i; ioptr[12] = t1r; ioptr[13] = t1i; ioptr[14] = f6r; ioptr[15] = f6i; } static void bfR2(SPFLOAT *ioptr, int M, int NDiffU) { /*** 2nd radix 2 stage ***/ unsigned int pos; unsigned int posi; unsigned int pinc; unsigned int pnext; unsigned int NSameU; unsigned int SameUCnt; SPFLOAT *pstrt; SPFLOAT *p0r, *p1r, *p2r, *p3r; SPFLOAT f0r, f0i, f1r, f1i, f2r, f2i, f3r, f3i; SPFLOAT f4r, f4i, f5r, f5i, f6r, f6i, f7r, f7i; pinc = NDiffU * 2; /* 2 floats per complex */ pnext = pinc * 4; pos = 2; posi = pos + 1; NSameU = POW2(M) / 4 / NDiffU; /* 4 Us at a time */ pstrt = ioptr; p0r = pstrt; p1r = pstrt + pinc; p2r = p1r + pinc; p3r = p2r + pinc; /* Butterflys */ /* f0 - - f4 f1 - 1 - f5 f2 - - f6 f3 - 1 - f7 */ /* Butterflys */ /* f0 - - f4 f1 - 1 - f5 f2 - - f6 f3 - 1 - f7 */ for (SameUCnt = NSameU; SameUCnt > 0; SameUCnt--) { f0r = *p0r; f1r = *p1r; f0i = *(p0r + 1); f1i = *(p1r + 1); f2r = *p2r; f3r = *p3r; f2i = *(p2r + 1); f3i = *(p3r + 1); f4r = f0r + f1r; f4i = f0i + f1i; f5r = f0r - f1r; f5i = f0i - f1i; f6r = f2r + f3r; f6i = f2i + f3i; f7r = f2r - f3r; f7i = f2i - f3i; *p0r = f4r; *(p0r + 1) = f4i; *p1r = f5r; *(p1r + 1) = f5i; *p2r = f6r; *(p2r + 1) = f6i; *p3r = f7r; *(p3r + 1) = f7i; f0r = *(p0r + pos); f1i = *(p1r + posi); f0i = *(p0r + posi); f1r = *(p1r + pos); f2r = *(p2r + pos); f3i = *(p3r + posi); f2i = *(p2r + posi); f3r = *(p3r + pos); f4r = f0r + f1i; f4i = f0i - f1r; f5r = f0r - f1i; f5i = f0i + f1r; f6r = f2r + f3i; f6i = f2i - f3r; f7r = f2r - f3i; f7i = f2i + f3r; *(p0r + pos) = f4r; *(p0r + posi) = f4i; *(p1r + pos) = f5r; *(p1r + posi) = f5i; *(p2r + pos) = f6r; *(p2r + posi) = f6i; *(p3r + pos) = f7r; *(p3r + posi) = f7i; p0r += pnext; p1r += pnext; p2r += pnext; p3r += pnext; } } static void bfR4(SPFLOAT *ioptr, int M, int NDiffU) { /*** 1 radix 4 stage ***/ unsigned int pos; unsigned int posi; unsigned int pinc; unsigned int pnext; unsigned int pnexti; unsigned int NSameU; unsigned int SameUCnt; SPFLOAT *pstrt; SPFLOAT *p0r, *p1r, *p2r, *p3r; SPFLOAT w1r = 1.0 / MYROOT2; /* cos(pi/4) */ SPFLOAT f0r, f0i, f1r, f1i, f2r, f2i, f3r, f3i; SPFLOAT f4r, f4i, f5r, f5i, f6r, f6i, f7r, f7i; SPFLOAT t1r, t1i; const SPFLOAT Two = 2.0; pinc = NDiffU * 2; /* 2 floats per complex */ pnext = pinc * 4; pnexti = pnext + 1; pos = 2; posi = pos + 1; NSameU = POW2(M) / 4 / NDiffU; /* 4 pts per butterfly */ pstrt = ioptr; p0r = pstrt; p1r = pstrt + pinc; p2r = p1r + pinc; p3r = p2r + pinc; /* Butterflys */ /* f0 - - f0 - - f4 f1 - 1 - f5 - - f5 f2 - - f6 - 1 - f6 f3 - 1 - f3 - -i - f7 */ /* Butterflys */ /* f0 - - f4 - - f4 f1 - -i - t1 - - f5 f2 - - f2 - w1 - f6 f3 - -i - f7 - iw1- f7 */ f0r = *p0r; f1r = *p1r; f2r = *p2r; f3r = *p3r; f0i = *(p0r + 1); f1i = *(p1r + 1); f2i = *(p2r + 1); f3i = *(p3r + 1); f5r = f0r - f1r; f5i = f0i - f1i; f0r = f0r + f1r; f0i = f0i + f1i; f6r = f2r + f3r; f6i = f2i + f3i; f3r = f2r - f3r; f3i = f2i - f3i; for (SameUCnt = NSameU - 1; SameUCnt > 0; SameUCnt--) { f7r = f5r - f3i; f7i = f5i + f3r; f5r = f5r + f3i; f5i = f5i - f3r; f4r = f0r + f6r; f4i = f0i + f6i; f6r = f0r - f6r; f6i = f0i - f6i; f2r = *(p2r + pos); f2i = *(p2r + posi); f1r = *(p1r + pos); f1i = *(p1r + posi); f3i = *(p3r + posi); f0r = *(p0r + pos); f3r = *(p3r + pos); f0i = *(p0r + posi); *p3r = f7r; *p0r = f4r; *(p3r + 1) = f7i; *(p0r + 1) = f4i; *p1r = f5r; *p2r = f6r; *(p1r + 1) = f5i; *(p2r + 1) = f6i; f7r = f2r - f3i; f7i = f2i + f3r; f2r = f2r + f3i; f2i = f2i - f3r; f4r = f0r + f1i; f4i = f0i - f1r; t1r = f0r - f1i; t1i = f0i + f1r; f5r = t1r - f7r * w1r + f7i * w1r; f5i = t1i - f7r * w1r - f7i * w1r; f7r = t1r * Two - f5r; f7i = t1i * Two - f5i; f6r = f4r - f2r * w1r - f2i * w1r; f6i = f4i + f2r * w1r - f2i * w1r; f4r = f4r * Two - f6r; f4i = f4i * Two - f6i; f3r = *(p3r + pnext); f0r = *(p0r + pnext); f3i = *(p3r + pnexti); f0i = *(p0r + pnexti); f2r = *(p2r + pnext); f2i = *(p2r + pnexti); f1r = *(p1r + pnext); f1i = *(p1r + pnexti); *(p2r + pos) = f6r; *(p1r + pos) = f5r; *(p2r + posi) = f6i; *(p1r + posi) = f5i; *(p3r + pos) = f7r; *(p0r + pos) = f4r; *(p3r + posi) = f7i; *(p0r + posi) = f4i; f6r = f2r + f3r; f6i = f2i + f3i; f3r = f2r - f3r; f3i = f2i - f3i; f5r = f0r - f1r; f5i = f0i - f1i; f0r = f0r + f1r; f0i = f0i + f1i; p3r += pnext; p0r += pnext; p1r += pnext; p2r += pnext; } f7r = f5r - f3i; f7i = f5i + f3r; f5r = f5r + f3i; f5i = f5i - f3r; f4r = f0r + f6r; f4i = f0i + f6i; f6r = f0r - f6r; f6i = f0i - f6i; f2r = *(p2r + pos); f2i = *(p2r + posi); f1r = *(p1r + pos); f1i = *(p1r + posi); f3i = *(p3r + posi); f0r = *(p0r + pos); f3r = *(p3r + pos); f0i = *(p0r + posi); *p3r = f7r; *p0r = f4r; *(p3r + 1) = f7i; *(p0r + 1) = f4i; *p1r = f5r; *p2r = f6r; *(p1r + 1) = f5i; *(p2r + 1) = f6i; f7r = f2r - f3i; f7i = f2i + f3r; f2r = f2r + f3i; f2i = f2i - f3r; f4r = f0r + f1i; f4i = f0i - f1r; t1r = f0r - f1i; t1i = f0i + f1r; f5r = t1r - f7r * w1r + f7i * w1r; f5i = t1i - f7r * w1r - f7i * w1r; f7r = t1r * Two - f5r; f7i = t1i * Two - f5i; f6r = f4r - f2r * w1r - f2i * w1r; f6i = f4i + f2r * w1r - f2i * w1r; f4r = f4r * Two - f6r; f4i = f4i * Two - f6i; *(p2r + pos) = f6r; *(p1r + pos) = f5r; *(p2r + posi) = f6i; *(p1r + posi) = f5i; *(p3r + pos) = f7r; *(p0r + pos) = f4r; *(p3r + posi) = f7i; *(p0r + posi) = f4i; } static void bfstages(SPFLOAT *ioptr, int M, SPFLOAT *Utbl, int Ustride, int NDiffU, int StageCnt) { /*** RADIX 8 Stages ***/ unsigned int pos; unsigned int posi; unsigned int pinc; unsigned int pnext; unsigned int NSameU; int Uinc; int Uinc2; int Uinc4; unsigned int DiffUCnt; unsigned int SameUCnt; unsigned int U2toU3; SPFLOAT *pstrt; SPFLOAT *p0r, *p1r, *p2r, *p3r; SPFLOAT *u0r, *u0i, *u1r, *u1i, *u2r, *u2i; SPFLOAT w0r, w0i, w1r, w1i, w2r, w2i, w3r, w3i; SPFLOAT f0r, f0i, f1r, f1i, f2r, f2i, f3r, f3i; SPFLOAT f4r, f4i, f5r, f5i, f6r, f6i, f7r, f7i; SPFLOAT t0r, t0i, t1r, t1i; const SPFLOAT Two = 2.0; pinc = NDiffU * 2; /* 2 floats per complex */ pnext = pinc * 8; pos = pinc * 4; posi = pos + 1; NSameU = POW2(M) / 8 / NDiffU; /* 8 pts per butterfly */ Uinc = (int) NSameU * Ustride; Uinc2 = Uinc * 2; Uinc4 = Uinc * 4; U2toU3 = (POW2(M) / 8) * Ustride; for (; StageCnt > 0; StageCnt--) { u0r = &Utbl[0]; u0i = &Utbl[POW2(M - 2) * Ustride]; u1r = u0r; u1i = u0i; u2r = u0r; u2i = u0i; w0r = *u0r; w0i = *u0i; w1r = *u1r; w1i = *u1i; w2r = *u2r; w2i = *u2i; w3r = *(u2r + U2toU3); w3i = *(u2i - U2toU3); pstrt = ioptr; p0r = pstrt; p1r = pstrt + pinc; p2r = p1r + pinc; p3r = p2r + pinc; /* Butterflys */ /* f0 - - t0 - - f0 - - f0 f1 - w0- f1 - - f1 - - f1 f2 - - f2 - w1- f2 - - f4 f3 - w0- t1 - iw1- f3 - - f5 f4 - - t0 - - f4 - w2- t0 f5 - w0- f5 - - f5 - w3- t1 f6 - - f6 - w1- f6 - iw2- f6 f7 - w0- t1 - iw1- f7 - iw3- f7 */ for (DiffUCnt = NDiffU; DiffUCnt > 0; DiffUCnt--) { f0r = *p0r; f0i = *(p0r + 1); f1r = *p1r; f1i = *(p1r + 1); for (SameUCnt = NSameU - 1; SameUCnt > 0; SameUCnt--) { f2r = *p2r; f2i = *(p2r + 1); f3r = *p3r; f3i = *(p3r + 1); t0r = f0r + f1r * w0r + f1i * w0i; t0i = f0i - f1r * w0i + f1i * w0r; f1r = f0r * Two - t0r; f1i = f0i * Two - t0i; f4r = *(p0r + pos); f4i = *(p0r + posi); f5r = *(p1r + pos); f5i = *(p1r + posi); f6r = *(p2r + pos); f6i = *(p2r + posi); f7r = *(p3r + pos); f7i = *(p3r + posi); t1r = f2r - f3r * w0r - f3i * w0i; t1i = f2i + f3r * w0i - f3i * w0r; f2r = f2r * Two - t1r; f2i = f2i * Two - t1i; f0r = t0r + f2r * w1r + f2i * w1i; f0i = t0i - f2r * w1i + f2i * w1r; f2r = t0r * Two - f0r; f2i = t0i * Two - f0i; f3r = f1r + t1r * w1i - t1i * w1r; f3i = f1i + t1r * w1r + t1i * w1i; f1r = f1r * Two - f3r; f1i = f1i * Two - f3i; t0r = f4r + f5r * w0r + f5i * w0i; t0i = f4i - f5r * w0i + f5i * w0r; f5r = f4r * Two - t0r; f5i = f4i * Two - t0i; t1r = f6r - f7r * w0r - f7i * w0i; t1i = f6i + f7r * w0i - f7i * w0r; f6r = f6r * Two - t1r; f6i = f6i * Two - t1i; f4r = t0r + f6r * w1r + f6i * w1i; f4i = t0i - f6r * w1i + f6i * w1r; f6r = t0r * Two - f4r; f6i = t0i * Two - f4i; f7r = f5r + t1r * w1i - t1i * w1r; f7i = f5i + t1r * w1r + t1i * w1i; f5r = f5r * Two - f7r; f5i = f5i * Two - f7i; t0r = f0r - f4r * w2r - f4i * w2i; t0i = f0i + f4r * w2i - f4i * w2r; f0r = f0r * Two - t0r; f0i = f0i * Two - t0i; t1r = f1r - f5r * w3r - f5i * w3i; t1i = f1i + f5r * w3i - f5i * w3r; f1r = f1r * Two - t1r; f1i = f1i * Two - t1i; *(p0r + pos) = t0r; *(p1r + pos) = t1r; *(p0r + posi) = t0i; *(p1r + posi) = t1i; *p0r = f0r; *p1r = f1r; *(p0r + 1) = f0i; *(p1r + 1) = f1i; p0r += pnext; f0r = *p0r; f0i = *(p0r + 1); p1r += pnext; f1r = *p1r; f1i = *(p1r + 1); f4r = f2r - f6r * w2i + f6i * w2r; f4i = f2i - f6r * w2r - f6i * w2i; f6r = f2r * Two - f4r; f6i = f2i * Two - f4i; f5r = f3r - f7r * w3i + f7i * w3r; f5i = f3i - f7r * w3r - f7i * w3i; f7r = f3r * Two - f5r; f7i = f3i * Two - f5i; *p2r = f4r; *p3r = f5r; *(p2r + 1) = f4i; *(p3r + 1) = f5i; *(p2r + pos) = f6r; *(p3r + pos) = f7r; *(p2r + posi) = f6i; *(p3r + posi) = f7i; p2r += pnext; p3r += pnext; } f2r = *p2r; f2i = *(p2r + 1); f3r = *p3r; f3i = *(p3r + 1); t0r = f0r + f1r * w0r + f1i * w0i; t0i = f0i - f1r * w0i + f1i * w0r; f1r = f0r * Two - t0r; f1i = f0i * Two - t0i; f4r = *(p0r + pos); f4i = *(p0r + posi); f5r = *(p1r + pos); f5i = *(p1r + posi); f6r = *(p2r + pos); f6i = *(p2r + posi); f7r = *(p3r + pos); f7i = *(p3r + posi); t1r = f2r - f3r * w0r - f3i * w0i; t1i = f2i + f3r * w0i - f3i * w0r; f2r = f2r * Two - t1r; f2i = f2i * Two - t1i; f0r = t0r + f2r * w1r + f2i * w1i; f0i = t0i - f2r * w1i + f2i * w1r; f2r = t0r * Two - f0r; f2i = t0i * Two - f0i; f3r = f1r + t1r * w1i - t1i * w1r; f3i = f1i + t1r * w1r + t1i * w1i; f1r = f1r * Two - f3r; f1i = f1i * Two - f3i; if ((int) DiffUCnt == NDiffU / 2) Uinc4 = -Uinc4; u0r += Uinc4; u0i -= Uinc4; u1r += Uinc2; u1i -= Uinc2; u2r += Uinc; u2i -= Uinc; pstrt += 2; t0r = f4r + f5r * w0r + f5i * w0i; t0i = f4i - f5r * w0i + f5i * w0r; f5r = f4r * Two - t0r; f5i = f4i * Two - t0i; t1r = f6r - f7r * w0r - f7i * w0i; t1i = f6i + f7r * w0i - f7i * w0r; f6r = f6r * Two - t1r; f6i = f6i * Two - t1i; f4r = t0r + f6r * w1r + f6i * w1i; f4i = t0i - f6r * w1i + f6i * w1r; f6r = t0r * Two - f4r; f6i = t0i * Two - f4i; f7r = f5r + t1r * w1i - t1i * w1r; f7i = f5i + t1r * w1r + t1i * w1i; f5r = f5r * Two - f7r; f5i = f5i * Two - f7i; w0r = *u0r; w0i = *u0i; w1r = *u1r; w1i = *u1i; if ((int) DiffUCnt <= NDiffU / 2) w0r = -w0r; t0r = f0r - f4r * w2r - f4i * w2i; t0i = f0i + f4r * w2i - f4i * w2r; f0r = f0r * Two - t0r; f0i = f0i * Two - t0i; f4r = f2r - f6r * w2i + f6i * w2r; f4i = f2i - f6r * w2r - f6i * w2i; f6r = f2r * Two - f4r; f6i = f2i * Two - f4i; *(p0r + pos) = t0r; *p2r = f4r; *(p0r + posi) = t0i; *(p2r + 1) = f4i; w2r = *u2r; w2i = *u2i; *p0r = f0r; *(p2r + pos) = f6r; *(p0r + 1) = f0i; *(p2r + posi) = f6i; p0r = pstrt; p2r = pstrt + pinc + pinc; t1r = f1r - f5r * w3r - f5i * w3i; t1i = f1i + f5r * w3i - f5i * w3r; f1r = f1r * Two - t1r; f1i = f1i * Two - t1i; f5r = f3r - f7r * w3i + f7i * w3r; f5i = f3i - f7r * w3r - f7i * w3i; f7r = f3r * Two - f5r; f7i = f3i * Two - f5i; *(p1r + pos) = t1r; *p3r = f5r; *(p1r + posi) = t1i; *(p3r + 1) = f5i; w3r = *(u2r + U2toU3); w3i = *(u2i - U2toU3); *p1r = f1r; *(p3r + pos) = f7r; *(p1r + 1) = f1i; *(p3r + posi) = f7i; p1r = pstrt + pinc; p3r = p2r + pinc; } NSameU /= 8; Uinc /= 8; Uinc2 /= 8; Uinc4 = Uinc * 4; NDiffU *= 8; pinc *= 8; pnext *= 8; pos *= 8; posi = pos + 1; } } static void fftrecurs(SPFLOAT *ioptr, int M, SPFLOAT *Utbl, int Ustride, int NDiffU, int StageCnt) { /* recursive bfstages calls to maximize on chip cache efficiency */ int i1; if (M <= (int) MCACHE) /* fits on chip ? */ bfstages(ioptr, M, Utbl, Ustride, NDiffU, StageCnt); /* RADIX 8 Stages */ else { for (i1 = 0; i1 < 8; i1++) { fftrecurs(&ioptr[i1 * POW2(M - 3) * 2], M - 3, Utbl, 8 * Ustride, NDiffU, StageCnt - 1); /* RADIX 8 Stages */ } bfstages(ioptr, M, Utbl, Ustride, POW2(M - 3), 1); /* RADIX 8 Stage */ } } static void ffts1(SPFLOAT *ioptr, int M, SPFLOAT *Utbl, int16_t *BRLow) { /* Compute in-place complex fft on the rows of the input array */ /* INPUTS */ /* *ioptr = input data array */ /* M = log2 of fft size (ex M=10 for 1024 point fft) */ /* *Utbl = cosine table */ /* *BRLow = bit reversed counter table */ /* OUTPUTS */ /* *ioptr = output data array */ int StageCnt; int NDiffU; switch (M) { case 0: break; case 1: fft2pt(ioptr); /* a 2 pt fft */ break; case 2: fft4pt(ioptr); /* a 4 pt fft */ break; case 3: fft8pt(ioptr); /* an 8 pt fft */ break; default: bitrevR2(ioptr, M, BRLow); /* bit reverse and first radix 2 stage */ StageCnt = (M - 1) / 3; /* number of radix 8 stages */ NDiffU = 2; /* one radix 2 stage already complete */ if ((M - 1 - (StageCnt * 3)) == 1) { bfR2(ioptr, M, NDiffU); /* 1 radix 2 stage */ NDiffU *= 2; } if ((M - 1 - (StageCnt * 3)) == 2) { bfR4(ioptr, M, NDiffU); /* 1 radix 4 stage */ NDiffU *= 4; } if (M <= (int) MCACHE) bfstages(ioptr, M, Utbl, 1, NDiffU, StageCnt); /* RADIX 8 Stages */ else fftrecurs(ioptr, M, Utbl, 1, NDiffU, StageCnt); /* RADIX 8 Stages */ } } /****************** * parts of iffts1 * ******************/ static void scbitrevR2(SPFLOAT *ioptr, int M, int16_t *BRLow, SPFLOAT scale) { /*** scaled bit reverse and first radix 2 stage forward or inverse fft ***/ SPFLOAT f0r; SPFLOAT f0i; SPFLOAT f1r; SPFLOAT f1i; SPFLOAT f2r; SPFLOAT f2i; SPFLOAT f3r; SPFLOAT f3i; SPFLOAT f4r; SPFLOAT f4i; SPFLOAT f5r; SPFLOAT f5i; SPFLOAT f6r; SPFLOAT f6i; SPFLOAT f7r; SPFLOAT f7i; SPFLOAT t0r; SPFLOAT t0i; SPFLOAT t1r; SPFLOAT t1i; SPFLOAT *p0r; SPFLOAT *p1r; SPFLOAT *IOP; SPFLOAT *iolimit; int Colstart; int iCol; unsigned int posA; unsigned int posAi; unsigned int posB; unsigned int posBi; const unsigned int Nrems2 = POW2((M + 3) / 2); const unsigned int Nroot_1_ColInc = POW2(M) - Nrems2; const unsigned int Nroot_1 = POW2(M / 2 - 1) - 1; const unsigned int ColstartShift = (M + 1) / 2 + 1; posA = POW2(M); /* 1/2 of POW2(M) complexes */ posAi = posA + 1; posB = posA + 2; posBi = posB + 1; iolimit = ioptr + Nrems2; for (; ioptr < iolimit; ioptr += POW2(M / 2 + 1)) { for (Colstart = Nroot_1; Colstart >= 0; Colstart--) { iCol = Nroot_1; p0r = ioptr + Nroot_1_ColInc + BRLow[Colstart] * 4; IOP = ioptr + (Colstart << ColstartShift); p1r = IOP + BRLow[iCol] * 4; f0r = *(p0r); f0i = *(p0r + 1); f1r = *(p0r + posA); f1i = *(p0r + posAi); for (; iCol > Colstart;) { f2r = *(p0r + 2); f2i = *(p0r + (2 + 1)); f3r = *(p0r + posB); f3i = *(p0r + posBi); f4r = *(p1r); f4i = *(p1r + 1); f5r = *(p1r + posA); f5i = *(p1r + posAi); f6r = *(p1r + 2); f6i = *(p1r + (2 + 1)); f7r = *(p1r + posB); f7i = *(p1r + posBi); t0r = f0r + f1r; t0i = f0i + f1i; f1r = f0r - f1r; f1i = f0i - f1i; t1r = f2r + f3r; t1i = f2i + f3i; f3r = f2r - f3r; f3i = f2i - f3i; f0r = f4r + f5r; f0i = f4i + f5i; f5r = f4r - f5r; f5i = f4i - f5i; f2r = f6r + f7r; f2i = f6i + f7i; f7r = f6r - f7r; f7i = f6i - f7i; *(p1r) = scale * t0r; *(p1r + 1) = scale * t0i; *(p1r + 2) = scale * f1r; *(p1r + (2 + 1)) = scale * f1i; *(p1r + posA) = scale * t1r; *(p1r + posAi) = scale * t1i; *(p1r + posB) = scale * f3r; *(p1r + posBi) = scale * f3i; *(p0r) = scale * f0r; *(p0r + 1) = scale * f0i; *(p0r + 2) = scale * f5r; *(p0r + (2 + 1)) = scale * f5i; *(p0r + posA) = scale * f2r; *(p0r + posAi) = scale * f2i; *(p0r + posB) = scale * f7r; *(p0r + posBi) = scale * f7i; p0r -= Nrems2; f0r = *(p0r); f0i = *(p0r + 1); f1r = *(p0r + posA); f1i = *(p0r + posAi); iCol -= 1; p1r = IOP + BRLow[iCol] * 4; } f2r = *(p0r + 2); f2i = *(p0r + (2 + 1)); f3r = *(p0r + posB); f3i = *(p0r + posBi); t0r = f0r + f1r; t0i = f0i + f1i; f1r = f0r - f1r; f1i = f0i - f1i; t1r = f2r + f3r; t1i = f2i + f3i; f3r = f2r - f3r; f3i = f2i - f3i; *(p0r) = scale * t0r; *(p0r + 1) = scale * t0i; *(p0r + 2) = scale * f1r; *(p0r + (2 + 1)) = scale * f1i; *(p0r + posA) = scale * t1r; *(p0r + posAi) = scale * t1i; *(p0r + posB) = scale * f3r; *(p0r + posBi) = scale * f3i; } } } /* static void ifft2pt(SPFLOAT *ioptr, SPFLOAT scale) */ /* { */ /* /\*** RADIX 2 ifft ***\/ */ /* SPFLOAT f0r, f0i, f1r, f1i; */ /* SPFLOAT t0r, t0i; */ /* /\* bit reversed load *\/ */ /* f0r = ioptr[0]; */ /* f0i = ioptr[1]; */ /* f1r = ioptr[2]; */ /* f1i = ioptr[3]; */ /* /\* Butterflys *\/ */ /* /\* */ /* f0 - - t0 */ /* f1 - 1 - f1 */ /* *\/ */ /* t0r = f0r + f1r; */ /* t0i = f0i + f1i; */ /* f1r = f0r - f1r; */ /* f1i = f0i - f1i; */ /* /\* store result *\/ */ /* ioptr[0] = scale * t0r; */ /* ioptr[1] = scale * t0i; */ /* ioptr[2] = scale * f1r; */ /* ioptr[3] = scale * f1i; */ /* } */ /* static void ifft4pt(SPFLOAT *ioptr, SPFLOAT scale) */ /* { */ /* /\*** RADIX 4 ifft ***\/ */ /* SPFLOAT f0r, f0i, f1r, f1i, f2r, f2i, f3r, f3i; */ /* SPFLOAT t0r, t0i, t1r, t1i; */ /* /\* bit reversed load *\/ */ /* f0r = ioptr[0]; */ /* f0i = ioptr[1]; */ /* f1r = ioptr[4]; */ /* f1i = ioptr[5]; */ /* f2r = ioptr[2]; */ /* f2i = ioptr[3]; */ /* f3r = ioptr[6]; */ /* f3i = ioptr[7]; */ /* /\* Butterflys *\/ */ /* /\* */ /* f0 - - t0 - - f0 */ /* f1 - 1 - f1 - - f1 */ /* f2 - - f2 - 1 - f2 */ /* f3 - 1 - t1 - i - f3 */ /* *\/ */ /* t0r = f0r + f1r; */ /* t0i = f0i + f1i; */ /* f1r = f0r - f1r; */ /* f1i = f0i - f1i; */ /* t1r = f2r - f3r; */ /* t1i = f2i - f3i; */ /* f2r = f2r + f3r; */ /* f2i = f2i + f3i; */ /* f0r = t0r + f2r; */ /* f0i = t0i + f2i; */ /* f2r = t0r - f2r; */ /* f2i = t0i - f2i; */ /* f3r = f1r + t1i; */ /* f3i = f1i - t1r; */ /* f1r = f1r - t1i; */ /* f1i = f1i + t1r; */ /* /\* store result *\/ */ /* ioptr[0] = scale * f0r; */ /* ioptr[1] = scale * f0i; */ /* ioptr[2] = scale * f1r; */ /* ioptr[3] = scale * f1i; */ /* ioptr[4] = scale * f2r; */ /* ioptr[5] = scale * f2i; */ /* ioptr[6] = scale * f3r; */ /* ioptr[7] = scale * f3i; */ /* } */ /* static void ifft8pt(SPFLOAT *ioptr, SPFLOAT scale) */ /* { */ /* /\*** RADIX 8 ifft ***\/ */ /* SPFLOAT w0r = 1.0 / MYROOT2; /\* cos(pi/4) *\/ */ /* SPFLOAT f0r, f0i, f1r, f1i, f2r, f2i, f3r, f3i; */ /* SPFLOAT f4r, f4i, f5r, f5i, f6r, f6i, f7r, f7i; */ /* SPFLOAT t0r, t0i, t1r, t1i; */ /* const SPFLOAT Two = 2.0; */ /* /\* bit reversed load *\/ */ /* f0r = ioptr[0]; */ /* f0i = ioptr[1]; */ /* f1r = ioptr[8]; */ /* f1i = ioptr[9]; */ /* f2r = ioptr[4]; */ /* f2i = ioptr[5]; */ /* f3r = ioptr[12]; */ /* f3i = ioptr[13]; */ /* f4r = ioptr[2]; */ /* f4i = ioptr[3]; */ /* f5r = ioptr[10]; */ /* f5i = ioptr[11]; */ /* f6r = ioptr[6]; */ /* f6i = ioptr[7]; */ /* f7r = ioptr[14]; */ /* f7i = ioptr[15]; */ /* /\* Butterflys *\/ */ /* /\* */ /* f0 - - t0 - - f0 - - f0 */ /* f1 - 1 - f1 - - f1 - - f1 */ /* f2 - - f2 - 1 - f2 - - f2 */ /* f3 - 1 - t1 - i - f3 - - f3 */ /* f4 - - t0 - - f4 - 1 - t0 */ /* f5 - 1 - f5 - - f5 - w3 - f4 */ /* f6 - - f6 - 1 - f6 - i - t1 */ /* f7 - 1 - t1 - i - f7 - iw3- f6 */ /* *\/ */ /* t0r = f0r + f1r; */ /* t0i = f0i + f1i; */ /* f1r = f0r - f1r; */ /* f1i = f0i - f1i; */ /* t1r = f2r - f3r; */ /* t1i = f2i - f3i; */ /* f2r = f2r + f3r; */ /* f2i = f2i + f3i; */ /* f0r = t0r + f2r; */ /* f0i = t0i + f2i; */ /* f2r = t0r - f2r; */ /* f2i = t0i - f2i; */ /* f3r = f1r + t1i; */ /* f3i = f1i - t1r; */ /* f1r = f1r - t1i; */ /* f1i = f1i + t1r; */ /* t0r = f4r + f5r; */ /* t0i = f4i + f5i; */ /* f5r = f4r - f5r; */ /* f5i = f4i - f5i; */ /* t1r = f6r - f7r; */ /* t1i = f6i - f7i; */ /* f6r = f6r + f7r; */ /* f6i = f6i + f7i; */ /* f4r = t0r + f6r; */ /* f4i = t0i + f6i; */ /* f6r = t0r - f6r; */ /* f6i = t0i - f6i; */ /* f7r = f5r + t1i; */ /* f7i = f5i - t1r; */ /* f5r = f5r - t1i; */ /* f5i = f5i + t1r; */ /* t0r = f0r - f4r; */ /* t0i = f0i - f4i; */ /* f0r = f0r + f4r; */ /* f0i = f0i + f4i; */ /* t1r = f2r + f6i; */ /* t1i = f2i - f6r; */ /* f2r = f2r - f6i; */ /* f2i = f2i + f6r; */ /* f4r = f1r - f5r * w0r + f5i * w0r; */ /* f4i = f1i - f5r * w0r - f5i * w0r; */ /* f1r = f1r * Two - f4r; */ /* f1i = f1i * Two - f4i; */ /* f6r = f3r + f7r * w0r + f7i * w0r; */ /* f6i = f3i - f7r * w0r + f7i * w0r; */ /* f3r = f3r * Two - f6r; */ /* f3i = f3i * Two - f6i; */ /* /\* store result *\/ */ /* ioptr[0] = scale * f0r; */ /* ioptr[1] = scale * f0i; */ /* ioptr[2] = scale * f1r; */ /* ioptr[3] = scale * f1i; */ /* ioptr[4] = scale * f2r; */ /* ioptr[5] = scale * f2i; */ /* ioptr[6] = scale * f3r; */ /* ioptr[7] = scale * f3i; */ /* ioptr[8] = scale * t0r; */ /* ioptr[9] = scale * t0i; */ /* ioptr[10] = scale * f4r; */ /* ioptr[11] = scale * f4i; */ /* ioptr[12] = scale * t1r; */ /* ioptr[13] = scale * t1i; */ /* ioptr[14] = scale * f6r; */ /* ioptr[15] = scale * f6i; */ /* } */ static void ibfR2(SPFLOAT *ioptr, int M, int NDiffU) { /*** 2nd radix 2 stage ***/ unsigned int pos; unsigned int posi; unsigned int pinc; unsigned int pnext; unsigned int NSameU; unsigned int SameUCnt; SPFLOAT *pstrt; SPFLOAT *p0r, *p1r, *p2r, *p3r; SPFLOAT f0r, f0i, f1r, f1i, f2r, f2i, f3r, f3i; SPFLOAT f4r, f4i, f5r, f5i, f6r, f6i, f7r, f7i; pinc = NDiffU * 2; /* 2 floats per complex */ pnext = pinc * 4; pos = 2; posi = pos + 1; NSameU = POW2(M) / 4 / NDiffU; /* 4 Us at a time */ pstrt = ioptr; p0r = pstrt; p1r = pstrt + pinc; p2r = p1r + pinc; p3r = p2r + pinc; /* Butterflys */ /* f0 - - f4 f1 - 1 - f5 f2 - - f6 f3 - 1 - f7 */ /* Butterflys */ /* f0 - - f4 f1 - 1 - f5 f2 - - f6 f3 - 1 - f7 */ for (SameUCnt = NSameU; SameUCnt > 0; SameUCnt--) { f0r = *p0r; f1r = *p1r; f0i = *(p0r + 1); f1i = *(p1r + 1); f2r = *p2r; f3r = *p3r; f2i = *(p2r + 1); f3i = *(p3r + 1); f4r = f0r + f1r; f4i = f0i + f1i; f5r = f0r - f1r; f5i = f0i - f1i; f6r = f2r + f3r; f6i = f2i + f3i; f7r = f2r - f3r; f7i = f2i - f3i; *p0r = f4r; *(p0r + 1) = f4i; *p1r = f5r; *(p1r + 1) = f5i; *p2r = f6r; *(p2r + 1) = f6i; *p3r = f7r; *(p3r + 1) = f7i; f0r = *(p0r + pos); f1i = *(p1r + posi); f0i = *(p0r + posi); f1r = *(p1r + pos); f2r = *(p2r + pos); f3i = *(p3r + posi); f2i = *(p2r + posi); f3r = *(p3r + pos); f4r = f0r - f1i; f4i = f0i + f1r; f5r = f0r + f1i; f5i = f0i - f1r; f6r = f2r - f3i; f6i = f2i + f3r; f7r = f2r + f3i; f7i = f2i - f3r; *(p0r + pos) = f4r; *(p0r + posi) = f4i; *(p1r + pos) = f5r; *(p1r + posi) = f5i; *(p2r + pos) = f6r; *(p2r + posi) = f6i; *(p3r + pos) = f7r; *(p3r + posi) = f7i; p0r += pnext; p1r += pnext; p2r += pnext; p3r += pnext; } } static void ibfR4(SPFLOAT *ioptr, int M, int NDiffU) { /*** 1 radix 4 stage ***/ unsigned int pos; unsigned int posi; unsigned int pinc; unsigned int pnext; unsigned int pnexti; unsigned int NSameU; unsigned int SameUCnt; SPFLOAT *pstrt; SPFLOAT *p0r, *p1r, *p2r, *p3r; SPFLOAT w1r = 1.0 / MYROOT2; /* cos(pi/4) */ SPFLOAT f0r, f0i, f1r, f1i, f2r, f2i, f3r, f3i; SPFLOAT f4r, f4i, f5r, f5i, f6r, f6i, f7r, f7i; SPFLOAT t1r, t1i; const SPFLOAT Two = 2.0; pinc = NDiffU * 2; /* 2 floats per complex */ pnext = pinc * 4; pnexti = pnext + 1; pos = 2; posi = pos + 1; NSameU = POW2(M) / 4 / NDiffU; /* 4 pts per butterfly */ pstrt = ioptr; p0r = pstrt; p1r = pstrt + pinc; p2r = p1r + pinc; p3r = p2r + pinc; /* Butterflys */ /* f0 - - f0 - - f4 f1 - 1 - f5 - - f5 f2 - - f6 - 1 - f6 f3 - 1 - f3 - -i - f7 */ /* Butterflys */ /* f0 - - f4 - - f4 f1 - -i - t1 - - f5 f2 - - f2 - w1 - f6 f3 - -i - f7 - iw1- f7 */ f0r = *p0r; f1r = *p1r; f2r = *p2r; f3r = *p3r; f0i = *(p0r + 1); f1i = *(p1r + 1); f2i = *(p2r + 1); f3i = *(p3r + 1); f5r = f0r - f1r; f5i = f0i - f1i; f0r = f0r + f1r; f0i = f0i + f1i; f6r = f2r + f3r; f6i = f2i + f3i; f3r = f2r - f3r; f3i = f2i - f3i; for (SameUCnt = NSameU - 1; SameUCnt > 0; SameUCnt--) { f7r = f5r + f3i; f7i = f5i - f3r; f5r = f5r - f3i; f5i = f5i + f3r; f4r = f0r + f6r; f4i = f0i + f6i; f6r = f0r - f6r; f6i = f0i - f6i; f2r = *(p2r + pos); f2i = *(p2r + posi); f1r = *(p1r + pos); f1i = *(p1r + posi); f3i = *(p3r + posi); f0r = *(p0r + pos); f3r = *(p3r + pos); f0i = *(p0r + posi); *p3r = f7r; *p0r = f4r; *(p3r + 1) = f7i; *(p0r + 1) = f4i; *p1r = f5r; *p2r = f6r; *(p1r + 1) = f5i; *(p2r + 1) = f6i; f7r = f2r + f3i; f7i = f2i - f3r; f2r = f2r - f3i; f2i = f2i + f3r; f4r = f0r - f1i; f4i = f0i + f1r; t1r = f0r + f1i; t1i = f0i - f1r; f5r = t1r - f7r * w1r - f7i * w1r; f5i = t1i + f7r * w1r - f7i * w1r; f7r = t1r * Two - f5r; f7i = t1i * Two - f5i; f6r = f4r - f2r * w1r + f2i * w1r; f6i = f4i - f2r * w1r - f2i * w1r; f4r = f4r * Two - f6r; f4i = f4i * Two - f6i; f3r = *(p3r + pnext); f0r = *(p0r + pnext); f3i = *(p3r + pnexti); f0i = *(p0r + pnexti); f2r = *(p2r + pnext); f2i = *(p2r + pnexti); f1r = *(p1r + pnext); f1i = *(p1r + pnexti); *(p2r + pos) = f6r; *(p1r + pos) = f5r; *(p2r + posi) = f6i; *(p1r + posi) = f5i; *(p3r + pos) = f7r; *(p0r + pos) = f4r; *(p3r + posi) = f7i; *(p0r + posi) = f4i; f6r = f2r + f3r; f6i = f2i + f3i; f3r = f2r - f3r; f3i = f2i - f3i; f5r = f0r - f1r; f5i = f0i - f1i; f0r = f0r + f1r; f0i = f0i + f1i; p3r += pnext; p0r += pnext; p1r += pnext; p2r += pnext; } f7r = f5r + f3i; f7i = f5i - f3r; f5r = f5r - f3i; f5i = f5i + f3r; f4r = f0r + f6r; f4i = f0i + f6i; f6r = f0r - f6r; f6i = f0i - f6i; f2r = *(p2r + pos); f2i = *(p2r + posi); f1r = *(p1r + pos); f1i = *(p1r + posi); f3i = *(p3r + posi); f0r = *(p0r + pos); f3r = *(p3r + pos); f0i = *(p0r + posi); *p3r = f7r; *p0r = f4r; *(p3r + 1) = f7i; *(p0r + 1) = f4i; *p1r = f5r; *p2r = f6r; *(p1r + 1) = f5i; *(p2r + 1) = f6i; f7r = f2r + f3i; f7i = f2i - f3r; f2r = f2r - f3i; f2i = f2i + f3r; f4r = f0r - f1i; f4i = f0i + f1r; t1r = f0r + f1i; t1i = f0i - f1r; f5r = t1r - f7r * w1r - f7i * w1r; f5i = t1i + f7r * w1r - f7i * w1r; f7r = t1r * Two - f5r; f7i = t1i * Two - f5i; f6r = f4r - f2r * w1r + f2i * w1r; f6i = f4i - f2r * w1r - f2i * w1r; f4r = f4r * Two - f6r; f4i = f4i * Two - f6i; *(p2r + pos) = f6r; *(p1r + pos) = f5r; *(p2r + posi) = f6i; *(p1r + posi) = f5i; *(p3r + pos) = f7r; *(p0r + pos) = f4r; *(p3r + posi) = f7i; *(p0r + posi) = f4i; } static void ibfstages(SPFLOAT *ioptr, int M, SPFLOAT *Utbl, int Ustride, int NDiffU, int StageCnt) { /*** RADIX 8 Stages ***/ unsigned int pos; unsigned int posi; unsigned int pinc; unsigned int pnext; unsigned int NSameU; int Uinc; int Uinc2; int Uinc4; unsigned int DiffUCnt; unsigned int SameUCnt; unsigned int U2toU3; SPFLOAT *pstrt; SPFLOAT *p0r, *p1r, *p2r, *p3r; SPFLOAT *u0r, *u0i, *u1r, *u1i, *u2r, *u2i; SPFLOAT w0r, w0i, w1r, w1i, w2r, w2i, w3r, w3i; SPFLOAT f0r, f0i, f1r, f1i, f2r, f2i, f3r, f3i; SPFLOAT f4r, f4i, f5r, f5i, f6r, f6i, f7r, f7i; SPFLOAT t0r, t0i, t1r, t1i; const SPFLOAT Two = 2.0; pinc = NDiffU * 2; /* 2 floats per complex */ pnext = pinc * 8; pos = pinc * 4; posi = pos + 1; NSameU = POW2(M) / 8 / NDiffU; /* 8 pts per butterfly */ Uinc = (int) NSameU * Ustride; Uinc2 = Uinc * 2; Uinc4 = Uinc * 4; U2toU3 = (POW2(M) / 8) * Ustride; for (; StageCnt > 0; StageCnt--) { u0r = &Utbl[0]; u0i = &Utbl[POW2(M - 2) * Ustride]; u1r = u0r; u1i = u0i; u2r = u0r; u2i = u0i; w0r = *u0r; w0i = *u0i; w1r = *u1r; w1i = *u1i; w2r = *u2r; w2i = *u2i; w3r = *(u2r + U2toU3); w3i = *(u2i - U2toU3); pstrt = ioptr; p0r = pstrt; p1r = pstrt + pinc; p2r = p1r + pinc; p3r = p2r + pinc; /* Butterflys */ /* f0 - - t0 - - f0 - - f0 f1 - w0- f1 - - f1 - - f1 f2 - - f2 - w1- f2 - - f4 f3 - w0- t1 - iw1- f3 - - f5 f4 - - t0 - - f4 - w2- t0 f5 - w0- f5 - - f5 - w3- t1 f6 - - f6 - w1- f6 - iw2- f6 f7 - w0- t1 - iw1- f7 - iw3- f7 */ for (DiffUCnt = NDiffU; DiffUCnt > 0; DiffUCnt--) { f0r = *p0r; f0i = *(p0r + 1); f1r = *p1r; f1i = *(p1r + 1); for (SameUCnt = NSameU - 1; SameUCnt > 0; SameUCnt--) { f2r = *p2r; f2i = *(p2r + 1); f3r = *p3r; f3i = *(p3r + 1); t0r = f0r + f1r * w0r - f1i * w0i; t0i = f0i + f1r * w0i + f1i * w0r; f1r = f0r * Two - t0r; f1i = f0i * Two - t0i; f4r = *(p0r + pos); f4i = *(p0r + posi); f5r = *(p1r + pos); f5i = *(p1r + posi); f6r = *(p2r + pos); f6i = *(p2r + posi); f7r = *(p3r + pos); f7i = *(p3r + posi); t1r = f2r - f3r * w0r + f3i * w0i; t1i = f2i - f3r * w0i - f3i * w0r; f2r = f2r * Two - t1r; f2i = f2i * Two - t1i; f0r = t0r + f2r * w1r - f2i * w1i; f0i = t0i + f2r * w1i + f2i * w1r; f2r = t0r * Two - f0r; f2i = t0i * Two - f0i; f3r = f1r + t1r * w1i + t1i * w1r; f3i = f1i - t1r * w1r + t1i * w1i; f1r = f1r * Two - f3r; f1i = f1i * Two - f3i; t0r = f4r + f5r * w0r - f5i * w0i; t0i = f4i + f5r * w0i + f5i * w0r; f5r = f4r * Two - t0r; f5i = f4i * Two - t0i; t1r = f6r - f7r * w0r + f7i * w0i; t1i = f6i - f7r * w0i - f7i * w0r; f6r = f6r * Two - t1r; f6i = f6i * Two - t1i; f4r = t0r + f6r * w1r - f6i * w1i; f4i = t0i + f6r * w1i + f6i * w1r; f6r = t0r * Two - f4r; f6i = t0i * Two - f4i; f7r = f5r + t1r * w1i + t1i * w1r; f7i = f5i - t1r * w1r + t1i * w1i; f5r = f5r * Two - f7r; f5i = f5i * Two - f7i; t0r = f0r - f4r * w2r + f4i * w2i; t0i = f0i - f4r * w2i - f4i * w2r; f0r = f0r * Two - t0r; f0i = f0i * Two - t0i; t1r = f1r - f5r * w3r + f5i * w3i; t1i = f1i - f5r * w3i - f5i * w3r; f1r = f1r * Two - t1r; f1i = f1i * Two - t1i; *(p0r + pos) = t0r; *(p0r + posi) = t0i; *p0r = f0r; *(p0r + 1) = f0i; p0r += pnext; f0r = *p0r; f0i = *(p0r + 1); *(p1r + pos) = t1r; *(p1r + posi) = t1i; *p1r = f1r; *(p1r + 1) = f1i; p1r += pnext; f1r = *p1r; f1i = *(p1r + 1); f4r = f2r - f6r * w2i - f6i * w2r; f4i = f2i + f6r * w2r - f6i * w2i; f6r = f2r * Two - f4r; f6i = f2i * Two - f4i; f5r = f3r - f7r * w3i - f7i * w3r; f5i = f3i + f7r * w3r - f7i * w3i; f7r = f3r * Two - f5r; f7i = f3i * Two - f5i; *p2r = f4r; *(p2r + 1) = f4i; *(p2r + pos) = f6r; *(p2r + posi) = f6i; p2r += pnext; *p3r = f5r; *(p3r + 1) = f5i; *(p3r + pos) = f7r; *(p3r + posi) = f7i; p3r += pnext; } f2r = *p2r; f2i = *(p2r + 1); f3r = *p3r; f3i = *(p3r + 1); t0r = f0r + f1r * w0r - f1i * w0i; t0i = f0i + f1r * w0i + f1i * w0r; f1r = f0r * Two - t0r; f1i = f0i * Two - t0i; f4r = *(p0r + pos); f4i = *(p0r + posi); f5r = *(p1r + pos); f5i = *(p1r + posi); f6r = *(p2r + pos); f6i = *(p2r + posi); f7r = *(p3r + pos); f7i = *(p3r + posi); t1r = f2r - f3r * w0r + f3i * w0i; t1i = f2i - f3r * w0i - f3i * w0r; f2r = f2r * Two - t1r; f2i = f2i * Two - t1i; f0r = t0r + f2r * w1r - f2i * w1i; f0i = t0i + f2r * w1i + f2i * w1r; f2r = t0r * Two - f0r; f2i = t0i * Two - f0i; f3r = f1r + t1r * w1i + t1i * w1r; f3i = f1i - t1r * w1r + t1i * w1i; f1r = f1r * Two - f3r; f1i = f1i * Two - f3i; if ((int) DiffUCnt == NDiffU / 2) Uinc4 = -Uinc4; u0r += Uinc4; u0i -= Uinc4; u1r += Uinc2; u1i -= Uinc2; u2r += Uinc; u2i -= Uinc; pstrt += 2; t0r = f4r + f5r * w0r - f5i * w0i; t0i = f4i + f5r * w0i + f5i * w0r; f5r = f4r * Two - t0r; f5i = f4i * Two - t0i; t1r = f6r - f7r * w0r + f7i * w0i; t1i = f6i - f7r * w0i - f7i * w0r; f6r = f6r * Two - t1r; f6i = f6i * Two - t1i; f4r = t0r + f6r * w1r - f6i * w1i; f4i = t0i + f6r * w1i + f6i * w1r; f6r = t0r * Two - f4r; f6i = t0i * Two - f4i; f7r = f5r + t1r * w1i + t1i * w1r; f7i = f5i - t1r * w1r + t1i * w1i; f5r = f5r * Two - f7r; f5i = f5i * Two - f7i; w0r = *u0r; w0i = *u0i; w1r = *u1r; w1i = *u1i; if ((int) DiffUCnt <= NDiffU / 2) w0r = -w0r; t0r = f0r - f4r * w2r + f4i * w2i; t0i = f0i - f4r * w2i - f4i * w2r; f0r = f0r * Two - t0r; f0i = f0i * Two - t0i; f4r = f2r - f6r * w2i - f6i * w2r; f4i = f2i + f6r * w2r - f6i * w2i; f6r = f2r * Two - f4r; f6i = f2i * Two - f4i; *(p0r + pos) = t0r; *p2r = f4r; *(p0r + posi) = t0i; *(p2r + 1) = f4i; w2r = *u2r; w2i = *u2i; *p0r = f0r; *(p2r + pos) = f6r; *(p0r + 1) = f0i; *(p2r + posi) = f6i; p0r = pstrt; p2r = pstrt + pinc + pinc; t1r = f1r - f5r * w3r + f5i * w3i; t1i = f1i - f5r * w3i - f5i * w3r; f1r = f1r * Two - t1r; f1i = f1i * Two - t1i; f5r = f3r - f7r * w3i - f7i * w3r; f5i = f3i + f7r * w3r - f7i * w3i; f7r = f3r * Two - f5r; f7i = f3i * Two - f5i; *(p1r + pos) = t1r; *p3r = f5r; *(p1r + posi) = t1i; *(p3r + 1) = f5i; w3r = *(u2r + U2toU3); w3i = *(u2i - U2toU3); *p1r = f1r; *(p3r + pos) = f7r; *(p1r + 1) = f1i; *(p3r + posi) = f7i; p1r = pstrt + pinc; p3r = p2r + pinc; } NSameU /= 8; Uinc /= 8; Uinc2 /= 8; Uinc4 = Uinc * 4; NDiffU *= 8; pinc *= 8; pnext *= 8; pos *= 8; posi = pos + 1; } } static void ifftrecurs(SPFLOAT *ioptr, int M, SPFLOAT *Utbl, int Ustride, int NDiffU, int StageCnt) { /* recursive bfstages calls to maximize on chip cache efficiency */ int i1; if (M <= (int) MCACHE) ibfstages(ioptr, M, Utbl, Ustride, NDiffU, StageCnt); /* RADIX 8 Stages */ else { for (i1 = 0; i1 < 8; i1++) { ifftrecurs(&ioptr[i1 * POW2(M - 3) * 2], M - 3, Utbl, 8 * Ustride, NDiffU, StageCnt - 1); /* RADIX 8 Stages */ } ibfstages(ioptr, M, Utbl, Ustride, POW2(M - 3), 1); /* RADIX 8 Stage */ } } /* Paul: not currently used, so commented out */ /* static void iffts1(SPFLOAT *ioptr, int M, SPFLOAT *Utbl, int16_t *BRLow) */ /* { */ /* /\* Compute in-place inverse complex fft on the rows of the input array *\/ */ /* /\* INPUTS *\/ */ /* /\* *ioptr = input data array *\/ */ /* /\* M = log2 of fft size *\/ */ /* /\* *Utbl = cosine table *\/ */ /* /\* *BRLow = bit reversed counter table *\/ */ /* /\* OUTPUTS *\/ */ /* /\* *ioptr = output data array *\/ */ /* int StageCnt; */ /* int NDiffU; */ /* const SPFLOAT scale = 1.0 / POW2(M); */ /* switch (M) { */ /* case 0: */ /* break; */ /* case 1: */ /* ifft2pt(ioptr, scale); /\* a 2 pt fft *\/ */ /* break; */ /* case 2: */ /* ifft4pt(ioptr, scale); /\* a 4 pt fft *\/ */ /* break; */ /* case 3: */ /* ifft8pt(ioptr, scale); /\* an 8 pt fft *\/ */ /* break; */ /* default: */ /* /\* bit reverse and first radix 2 stage *\/ */ /* scbitrevR2(ioptr, M, BRLow, scale); */ /* StageCnt = (M - 1) / 3; /\* number of radix 8 stages *\/ */ /* NDiffU = 2; /\* one radix 2 stage already complete *\/ */ /* if ((M - 1 - (StageCnt * 3)) == 1) { */ /* ibfR2(ioptr, M, NDiffU); /\* 1 radix 2 stage *\/ */ /* NDiffU *= 2; */ /* } */ /* if ((M - 1 - (StageCnt * 3)) == 2) { */ /* ibfR4(ioptr, M, NDiffU); /\* 1 radix 4 stage *\/ */ /* NDiffU *= 4; */ /* } */ /* if (M <= (int) MCACHE) */ /* ibfstages(ioptr, M, Utbl, 1, NDiffU, StageCnt); /\* RADIX 8 Stages *\/ */ /* else */ /* ifftrecurs(ioptr, M, Utbl, 1, NDiffU, StageCnt); /\* RADIX 8 Stages *\/ */ /* } */ /* } */ /****************** * parts of rffts1 * ******************/ static void rfft1pt(SPFLOAT *ioptr) { /*** RADIX 2 rfft ***/ SPFLOAT f0r, f0i; SPFLOAT t0r, t0i; /* bit reversed load */ f0r = ioptr[0]; f0i = ioptr[1]; /* finish rfft */ t0r = f0r + f0i; t0i = f0r - f0i; /* store result */ ioptr[0] = t0r; ioptr[1] = t0i; } static void rfft2pt(SPFLOAT *ioptr) { /*** RADIX 4 rfft ***/ SPFLOAT f0r, f0i, f1r, f1i; SPFLOAT t0r, t0i; /* bit reversed load */ f0r = ioptr[0]; f0i = ioptr[1]; f1r = ioptr[2]; f1i = ioptr[3]; /* Butterflys */ /* f0 - - t0 f1 - 1 - f1 */ t0r = f0r + f1r; t0i = f0i + f1i; f1r = f0r - f1r; f1i = f1i - f0i; /* finish rfft */ f0r = t0r + t0i; f0i = t0r - t0i; /* store result */ ioptr[0] = f0r; ioptr[1] = f0i; ioptr[2] = f1r; ioptr[3] = f1i; } static void rfft4pt(SPFLOAT *ioptr) { /*** RADIX 8 rfft ***/ SPFLOAT f0r, f0i, f1r, f1i, f2r, f2i, f3r, f3i; SPFLOAT t0r, t0i, t1r, t1i; SPFLOAT w0r = 1.0 / MYROOT2; /* cos(pi/4) */ const SPFLOAT Two = 2.0; const SPFLOAT scale = 0.5; /* bit reversed load */ f0r = ioptr[0]; f0i = ioptr[1]; f1r = ioptr[4]; f1i = ioptr[5]; f2r = ioptr[2]; f2i = ioptr[3]; f3r = ioptr[6]; f3i = ioptr[7]; /* Butterflys */ /* f0 - - t0 - - f0 f1 - 1 - f1 - - f1 f2 - - f2 - 1 - f2 f3 - 1 - t1 - -i - f3 */ t0r = f0r + f1r; t0i = f0i + f1i; f1r = f0r - f1r; f1i = f0i - f1i; t1r = f2r - f3r; t1i = f2i - f3i; f2r = f2r + f3r; f2i = f2i + f3i; f0r = t0r + f2r; f0i = t0i + f2i; f2r = t0r - f2r; f2i = f2i - t0i; /* neg for rfft */ f3r = f1r - t1i; f3i = f1i + t1r; f1r = f1r + t1i; f1i = f1i - t1r; /* finish rfft */ t0r = f0r + f0i; /* compute Re(x[0]) */ t0i = f0r - f0i; /* compute Re(x[N/2]) */ t1r = f1r + f3r; t1i = f1i - f3i; f0r = f1i + f3i; f0i = f3r - f1r; f1r = t1r + w0r * f0r + w0r * f0i; f1i = t1i - w0r * f0r + w0r * f0i; f3r = Two * t1r - f1r; f3i = f1i - Two * t1i; /* store result */ ioptr[4] = f2r; ioptr[5] = f2i; ioptr[0] = t0r; ioptr[1] = t0i; ioptr[2] = scale * f1r; ioptr[3] = scale * f1i; ioptr[6] = scale * f3r; ioptr[7] = scale * f3i; } static void rfft8pt(SPFLOAT *ioptr) { /*** RADIX 16 rfft ***/ SPFLOAT w0r = 1.0 / MYROOT2; /* cos(pi/4) */ SPFLOAT w1r = MYCOSPID8; /* cos(pi/8) */ SPFLOAT w1i = MYSINPID8; /* sin(pi/8) */ SPFLOAT f0r, f0i, f1r, f1i, f2r, f2i, f3r, f3i; SPFLOAT f4r, f4i, f5r, f5i, f6r, f6i, f7r, f7i; SPFLOAT t0r, t0i, t1r, t1i; const SPFLOAT Two = 2.0; const SPFLOAT scale = 0.5; /* bit reversed load */ f0r = ioptr[0]; f0i = ioptr[1]; f1r = ioptr[8]; f1i = ioptr[9]; f2r = ioptr[4]; f2i = ioptr[5]; f3r = ioptr[12]; f3i = ioptr[13]; f4r = ioptr[2]; f4i = ioptr[3]; f5r = ioptr[10]; f5i = ioptr[11]; f6r = ioptr[6]; f6i = ioptr[7]; f7r = ioptr[14]; f7i = ioptr[15]; /* Butterflys */ /* f0 - - t0 - - f0 - - f0 f1 - 1 - f1 - - f1 - - f1 f2 - - f2 - 1 - f2 - - f2 f3 - 1 - t1 - -i - f3 - - f3 f4 - - t0 - - f4 - 1 - t0 f5 - 1 - f5 - - f5 - w3 - f4 f6 - - f6 - 1 - f6 - -i - t1 f7 - 1 - t1 - -i - f7 - iw3- f6 */ t0r = f0r + f1r; t0i = f0i + f1i; f1r = f0r - f1r; f1i = f0i - f1i; t1r = f2r - f3r; t1i = f2i - f3i; f2r = f2r + f3r; f2i = f2i + f3i; f0r = t0r + f2r; f0i = t0i + f2i; f2r = t0r - f2r; f2i = t0i - f2i; f3r = f1r - t1i; f3i = f1i + t1r; f1r = f1r + t1i; f1i = f1i - t1r; t0r = f4r + f5r; t0i = f4i + f5i; f5r = f4r - f5r; f5i = f4i - f5i; t1r = f6r - f7r; t1i = f6i - f7i; f6r = f6r + f7r; f6i = f6i + f7i; f4r = t0r + f6r; f4i = t0i + f6i; f6r = t0r - f6r; f6i = t0i - f6i; f7r = f5r - t1i; f7i = f5i + t1r; f5r = f5r + t1i; f5i = f5i - t1r; t0r = f0r - f4r; t0i = f4i - f0i; /* neg for rfft */ f0r = f0r + f4r; f0i = f0i + f4i; t1r = f2r - f6i; t1i = f2i + f6r; f2r = f2r + f6i; f2i = f2i - f6r; f4r = f1r - f5r * w0r - f5i * w0r; f4i = f1i + f5r * w0r - f5i * w0r; f1r = f1r * Two - f4r; f1i = f1i * Two - f4i; f6r = f3r + f7r * w0r - f7i * w0r; f6i = f3i + f7r * w0r + f7i * w0r; f3r = f3r * Two - f6r; f3i = f3i * Two - f6i; /* finish rfft */ f5r = f0r + f0i; /* compute Re(x[0]) */ f5i = f0r - f0i; /* compute Re(x[N/2]) */ f0r = f2r + t1r; f0i = f2i - t1i; f7r = f2i + t1i; f7i = t1r - f2r; f2r = f0r + w0r * f7r + w0r * f7i; f2i = f0i - w0r * f7r + w0r * f7i; t1r = Two * f0r - f2r; t1i = f2i - Two * f0i; f0r = f1r + f6r; f0i = f1i - f6i; f7r = f1i + f6i; f7i = f6r - f1r; f1r = f0r + w1r * f7r + w1i * f7i; f1i = f0i - w1i * f7r + w1r * f7i; f6r = Two * f0r - f1r; f6i = f1i - Two * f0i; f0r = f3r + f4r; f0i = f3i - f4i; f7r = f3i + f4i; f7i = f4r - f3r; f3r = f0r + w1i * f7r + w1r * f7i; f3i = f0i - w1r * f7r + w1i * f7i; f4r = Two * f0r - f3r; f4i = f3i - Two * f0i; /* store result */ ioptr[8] = t0r; ioptr[9] = t0i; ioptr[0] = f5r; ioptr[1] = f5i; ioptr[4] = scale * f2r; ioptr[5] = scale * f2i; ioptr[12] = scale * t1r; ioptr[13] = scale * t1i; ioptr[2] = scale * f1r; ioptr[3] = scale * f1i; ioptr[6] = scale * f3r; ioptr[7] = scale * f3i; ioptr[10] = scale * f4r; ioptr[11] = scale * f4i; ioptr[14] = scale * f6r; ioptr[15] = scale * f6i; } static void frstage(SPFLOAT *ioptr, int M, SPFLOAT *Utbl) { /* Finish RFFT */ unsigned int pos; unsigned int posi; unsigned int diffUcnt; SPFLOAT *p0r, *p1r; SPFLOAT *u0r, *u0i; SPFLOAT w0r, w0i; SPFLOAT f0r, f0i, f1r, f1i, f4r, f4i, f5r, f5i; SPFLOAT t0r, t0i, t1r, t1i; const SPFLOAT Two = 2.0; pos = POW2(M - 1); posi = pos + 1; p0r = ioptr; p1r = ioptr + pos / 2; u0r = Utbl + POW2(M - 3); w0r = *u0r; f0r = *(p0r); f0i = *(p0r + 1); f4r = *(p0r + pos); f4i = *(p0r + posi); f1r = *(p1r); f1i = *(p1r + 1); f5r = *(p1r + pos); f5i = *(p1r + posi); t0r = Two * f0r + Two * f0i; /* compute Re(x[0]) */ t0i = Two * f0r - Two * f0i; /* compute Re(x[N/2]) */ t1r = f4r + f4r; t1i = -f4i - f4i; f0r = f1r + f5r; f0i = f1i - f5i; f4r = f1i + f5i; f4i = f5r - f1r; f1r = f0r + w0r * f4r + w0r * f4i; f1i = f0i - w0r * f4r + w0r * f4i; f5r = Two * f0r - f1r; f5i = f1i - Two * f0i; *(p0r) = t0r; *(p0r + 1) = t0i; *(p0r + pos) = t1r; *(p0r + posi) = t1i; *(p1r) = f1r; *(p1r + 1) = f1i; *(p1r + pos) = f5r; *(p1r + posi) = f5i; u0r = Utbl + 1; u0i = Utbl + (POW2(M - 2) - 1); w0r = *u0r; w0i = *u0i; p0r = (ioptr + 2); p1r = (ioptr + (POW2(M - 2) - 1) * 2); /* Butterflys */ /* f0 - t0 - - f0 f5 - t1 - w0 - f5 f1 - t0 - - f1 f4 - t1 -iw0 - f4 */ for (diffUcnt = POW2(M - 3) - 1; diffUcnt > 0; diffUcnt--) { f0r = *(p0r); f0i = *(p0r + 1); f5r = *(p1r + pos); f5i = *(p1r + posi); f1r = *(p1r); f1i = *(p1r + 1); f4r = *(p0r + pos); f4i = *(p0r + posi); t0r = f0r + f5r; t0i = f0i - f5i; t1r = f0i + f5i; t1i = f5r - f0r; f0r = t0r + w0r * t1r + w0i * t1i; f0i = t0i - w0i * t1r + w0r * t1i; f5r = Two * t0r - f0r; f5i = f0i - Two * t0i; t0r = f1r + f4r; t0i = f1i - f4i; t1r = f1i + f4i; t1i = f4r - f1r; f1r = t0r + w0i * t1r + w0r * t1i; f1i = t0i - w0r * t1r + w0i * t1i; f4r = Two * t0r - f1r; f4i = f1i - Two * t0i; *(p0r) = f0r; *(p0r + 1) = f0i; *(p1r + pos) = f5r; *(p1r + posi) = f5i; w0r = *++u0r; w0i = *--u0i; *(p1r) = f1r; *(p1r + 1) = f1i; *(p0r + pos) = f4r; *(p0r + posi) = f4i; p0r += 2; p1r -= 2; } } static void rffts1(SPFLOAT *ioptr, int M, SPFLOAT *Utbl, int16_t *BRLow) { /* Compute in-place real fft on the rows of the input array */ /* The result is the complex spectra of the positive frequencies */ /* except the location for the first complex number contains the real */ /* values for DC and Nyquest */ /* INPUTS */ /* *ioptr = real input data array */ /* M = log2 of fft size */ /* *Utbl = cosine table */ /* *BRLow = bit reversed counter table */ /* OUTPUTS */ /* *ioptr = output data array in the following order */ /* Re(x[0]), Re(x[N/2]), Re(x[1]), Im(x[1]), Re(x[2]), Im(x[2]), */ /* ... Re(x[N/2-1]), Im(x[N/2-1]). */ SPFLOAT scale; int StageCnt; int NDiffU; M = M - 1; switch (M) { case -1: break; case 0: rfft1pt(ioptr); /* a 2 pt fft */ break; case 1: rfft2pt(ioptr); /* a 4 pt fft */ break; case 2: rfft4pt(ioptr); /* an 8 pt fft */ break; case 3: rfft8pt(ioptr); /* a 16 pt fft */ break; default: scale = 0.5; /* bit reverse and first radix 2 stage */ scbitrevR2(ioptr, M, BRLow, scale); StageCnt = (M - 1) / 3; /* number of radix 8 stages */ NDiffU = 2; /* one radix 2 stage already complete */ if ((M - 1 - (StageCnt * 3)) == 1) { bfR2(ioptr, M, NDiffU); /* 1 radix 2 stage */ NDiffU *= 2; } if ((M - 1 - (StageCnt * 3)) == 2) { bfR4(ioptr, M, NDiffU); /* 1 radix 4 stage */ NDiffU *= 4; } if (M <= (int) MCACHE) bfstages(ioptr, M, Utbl, 2, NDiffU, StageCnt); /* RADIX 8 Stages */ else fftrecurs(ioptr, M, Utbl, 2, NDiffU, StageCnt); /* RADIX 8 Stages */ frstage(ioptr, M + 1, Utbl); } } /******************* * parts of riffts1 * *******************/ static void rifft1pt(SPFLOAT *ioptr, SPFLOAT scale) { /*** RADIX 2 rifft ***/ SPFLOAT f0r, f0i; SPFLOAT t0r, t0i; /* bit reversed load */ f0r = ioptr[0]; f0i = ioptr[1]; /* finish rfft */ t0r = f0r + f0i; t0i = f0r - f0i; /* store result */ ioptr[0] = scale * t0r; ioptr[1] = scale * t0i; } static void rifft2pt(SPFLOAT *ioptr, SPFLOAT scale) { /*** RADIX 4 rifft ***/ SPFLOAT f0r, f0i, f1r, f1i; SPFLOAT t0r, t0i; const SPFLOAT Two = 2.0; /* bit reversed load */ t0r = ioptr[0]; t0i = ioptr[1]; f1r = Two * ioptr[2]; f1i = Two * ioptr[3]; /* start rifft */ f0r = t0r + t0i; f0i = t0r - t0i; /* Butterflys */ /* f0 - - t0 f1 - 1 - f1 */ t0r = f0r + f1r; t0i = f0i - f1i; f1r = f0r - f1r; f1i = f0i + f1i; /* store result */ ioptr[0] = scale * t0r; ioptr[1] = scale * t0i; ioptr[2] = scale * f1r; ioptr[3] = scale * f1i; } static void rifft4pt(SPFLOAT *ioptr, SPFLOAT scale) { /*** RADIX 8 rifft ***/ SPFLOAT f0r, f0i, f1r, f1i, f2r, f2i, f3r, f3i; SPFLOAT t0r, t0i, t1r, t1i; SPFLOAT w0r = 1.0 / MYROOT2; /* cos(pi/4) */ const SPFLOAT Two = 2.0; /* bit reversed load */ t0r = ioptr[0]; t0i = ioptr[1]; f2r = ioptr[2]; f2i = ioptr[3]; f1r = Two * ioptr[4]; f1i = Two * ioptr[5]; f3r = ioptr[6]; f3i = ioptr[7]; /* start rfft */ f0r = t0r + t0i; /* compute Re(x[0]) */ f0i = t0r - t0i; /* compute Re(x[N/2]) */ t1r = f2r + f3r; t1i = f2i - f3i; t0r = f2r - f3r; t0i = f2i + f3i; f2r = t1r - w0r * t0r - w0r * t0i; f2i = t1i + w0r * t0r - w0r * t0i; f3r = Two * t1r - f2r; f3i = f2i - Two * t1i; /* Butterflys */ /* f0 - - t0 - - f0 f1 - 1 - f1 - - f1 f2 - - f2 - 1 - f2 f3 - 1 - t1 - i - f3 */ t0r = f0r + f1r; t0i = f0i - f1i; f1r = f0r - f1r; f1i = f0i + f1i; t1r = f2r - f3r; t1i = f2i - f3i; f2r = f2r + f3r; f2i = f2i + f3i; f0r = t0r + f2r; f0i = t0i + f2i; f2r = t0r - f2r; f2i = t0i - f2i; f3r = f1r + t1i; f3i = f1i - t1r; f1r = f1r - t1i; f1i = f1i + t1r; /* store result */ ioptr[0] = scale * f0r; ioptr[1] = scale * f0i; ioptr[2] = scale * f1r; ioptr[3] = scale * f1i; ioptr[4] = scale * f2r; ioptr[5] = scale * f2i; ioptr[6] = scale * f3r; ioptr[7] = scale * f3i; } static void rifft8pt(SPFLOAT *ioptr, SPFLOAT scale) { /*** RADIX 16 rifft ***/ SPFLOAT w0r = (SPFLOAT) (1.0 / MYROOT2); /* cos(pi/4) */ SPFLOAT w1r = MYCOSPID8; /* cos(pi/8) */ SPFLOAT w1i = MYSINPID8; /* sin(pi/8) */ SPFLOAT f0r, f0i, f1r, f1i, f2r, f2i, f3r, f3i; SPFLOAT f4r, f4i, f5r, f5i, f6r, f6i, f7r, f7i; SPFLOAT t0r, t0i, t1r, t1i; const SPFLOAT Two = 2.0; /* bit reversed load */ t0r = ioptr[0]; t0i = ioptr[1]; f4r = ioptr[2]; f4i = ioptr[3]; f2r = ioptr[4]; f2i = ioptr[5]; f6r = ioptr[6]; f6i = ioptr[7]; f1r = Two * ioptr[8]; f1i = Two * ioptr[9]; f5r = ioptr[10]; f5i = ioptr[11]; f3r = ioptr[12]; f3i = ioptr[13]; f7r = ioptr[14]; f7i = ioptr[15]; /* start rfft */ f0r = t0r + t0i; /* compute Re(x[0]) */ f0i = t0r - t0i; /* compute Re(x[N/2]) */ t0r = f2r + f3r; t0i = f2i - f3i; t1r = f2r - f3r; t1i = f2i + f3i; f2r = t0r - w0r * t1r - w0r * t1i; f2i = t0i + w0r * t1r - w0r * t1i; f3r = Two * t0r - f2r; f3i = f2i - Two * t0i; t0r = f4r + f7r; t0i = f4i - f7i; t1r = f4r - f7r; t1i = f4i + f7i; f4r = t0r - w1i * t1r - w1r * t1i; f4i = t0i + w1r * t1r - w1i * t1i; f7r = Two * t0r - f4r; f7i = f4i - Two * t0i; t0r = f6r + f5r; t0i = f6i - f5i; t1r = f6r - f5r; t1i = f6i + f5i; f6r = t0r - w1r * t1r - w1i * t1i; f6i = t0i + w1i * t1r - w1r * t1i; f5r = Two * t0r - f6r; f5i = f6i - Two * t0i; /* Butterflys */ /* f0 - - t0 - - f0 - - f0 f1* - 1 - f1 - - f1 - - f1 f2 - - f2 - 1 - f2 - - f2 f3 - 1 - t1 - i - f3 - - f3 f4 - - t0 - - f4 - 1 - t0 f5 - 1 - f5 - - f5 - w3 - f4 f6 - - f6 - 1 - f6 - i - t1 f7 - 1 - t1 - i - f7 - iw3- f6 */ t0r = f0r + f1r; t0i = f0i - f1i; f1r = f0r - f1r; f1i = f0i + f1i; t1r = f2r - f3r; t1i = f2i - f3i; f2r = f2r + f3r; f2i = f2i + f3i; f0r = t0r + f2r; f0i = t0i + f2i; f2r = t0r - f2r; f2i = t0i - f2i; f3r = f1r + t1i; f3i = f1i - t1r; f1r = f1r - t1i; f1i = f1i + t1r; t0r = f4r + f5r; t0i = f4i + f5i; f5r = f4r - f5r; f5i = f4i - f5i; t1r = f6r - f7r; t1i = f6i - f7i; f6r = f6r + f7r; f6i = f6i + f7i; f4r = t0r + f6r; f4i = t0i + f6i; f6r = t0r - f6r; f6i = t0i - f6i; f7r = f5r + t1i; f7i = f5i - t1r; f5r = f5r - t1i; f5i = f5i + t1r; t0r = f0r - f4r; t0i = f0i - f4i; f0r = f0r + f4r; f0i = f0i + f4i; t1r = f2r + f6i; t1i = f2i - f6r; f2r = f2r - f6i; f2i = f2i + f6r; f4r = f1r - f5r * w0r + f5i * w0r; f4i = f1i - f5r * w0r - f5i * w0r; f1r = f1r * Two - f4r; f1i = f1i * Two - f4i; f6r = f3r + f7r * w0r + f7i * w0r; f6i = f3i - f7r * w0r + f7i * w0r; f3r = f3r * Two - f6r; f3i = f3i * Two - f6i; /* store result */ ioptr[0] = scale * f0r; ioptr[1] = scale * f0i; ioptr[2] = scale * f1r; ioptr[3] = scale * f1i; ioptr[4] = scale * f2r; ioptr[5] = scale * f2i; ioptr[6] = scale * f3r; ioptr[7] = scale * f3i; ioptr[8] = scale * t0r; ioptr[9] = scale * t0i; ioptr[10] = scale * f4r; ioptr[11] = scale * f4i; ioptr[12] = scale * t1r; ioptr[13] = scale * t1i; ioptr[14] = scale * f6r; ioptr[15] = scale * f6i; } static void ifrstage(SPFLOAT *ioptr, int M, SPFLOAT *Utbl) { /* Start RIFFT */ unsigned int pos; unsigned int posi; unsigned int diffUcnt; SPFLOAT *p0r, *p1r; SPFLOAT *u0r, *u0i; SPFLOAT w0r, w0i; SPFLOAT f0r, f0i, f1r, f1i, f4r, f4i, f5r, f5i; SPFLOAT t0r, t0i, t1r, t1i; const SPFLOAT Two = 2.0; pos = POW2(M - 1); posi = pos + 1; p0r = ioptr; p1r = ioptr + pos / 2; u0r = Utbl + POW2(M - 3); w0r = *u0r; f0r = *(p0r); f0i = *(p0r + 1); f4r = *(p0r + pos); f4i = *(p0r + posi); f1r = *(p1r); f1i = *(p1r + 1); f5r = *(p1r + pos); f5i = *(p1r + posi); t0r = f0r + f0i; t0i = f0r - f0i; t1r = f4r + f4r; t1i = -f4i - f4i; f0r = f1r + f5r; f0i = f1i - f5i; f4r = f1r - f5r; f4i = f1i + f5i; f1r = f0r - w0r * f4r - w0r * f4i; f1i = f0i + w0r * f4r - w0r * f4i; f5r = Two * f0r - f1r; f5i = f1i - Two * f0i; *(p0r) = t0r; *(p0r + 1) = t0i; *(p0r + pos) = t1r; *(p0r + posi) = t1i; *(p1r) = f1r; *(p1r + 1) = f1i; *(p1r + pos) = f5r; *(p1r + posi) = f5i; u0r = Utbl + 1; u0i = Utbl + (POW2(M - 2) - 1); w0r = *u0r; w0i = *u0i; p0r = (ioptr + 2); p1r = (ioptr + (POW2(M - 2) - 1) * 2); /* Butterflys */ /* f0 - t0 - f0 f1 - t1 -w0- f1 f2 - t0 - f2 f3 - t1 -iw0- f3 */ for (diffUcnt = POW2(M - 3) - 1; diffUcnt > 0; diffUcnt--) { f0r = *(p0r); f0i = *(p0r + 1); f5r = *(p1r + pos); f5i = *(p1r + posi); f1r = *(p1r); f1i = *(p1r + 1); f4r = *(p0r + pos); f4i = *(p0r + posi); t0r = f0r + f5r; t0i = f0i - f5i; t1r = f0r - f5r; t1i = f0i + f5i; f0r = t0r - w0i * t1r - w0r * t1i; f0i = t0i + w0r * t1r - w0i * t1i; f5r = Two * t0r - f0r; f5i = f0i - Two * t0i; t0r = f1r + f4r; t0i = f1i - f4i; t1r = f1r - f4r; t1i = f1i + f4i; f1r = t0r - w0r * t1r - w0i * t1i; f1i = t0i + w0i * t1r - w0r * t1i; f4r = Two * t0r - f1r; f4i = f1i - Two * t0i; *(p0r) = f0r; *(p0r + 1) = f0i; *(p1r + pos) = f5r; *(p1r + posi) = f5i; w0r = *++u0r; w0i = *--u0i; *(p1r) = f1r; *(p1r + 1) = f1i; *(p0r + pos) = f4r; *(p0r + posi) = f4i; p0r += 2; p1r -= 2; } } static void riffts1(SPFLOAT *ioptr, int M, SPFLOAT *Utbl, int16_t *BRLow) { /* Compute in-place real ifft on the rows of the input array */ /* data order as from rffts1 */ /* INPUTS */ /* *ioptr = input data array in the following order */ /* M = log2 of fft size */ /* Re(x[0]), Re(x[N/2]), Re(x[1]), Im(x[1]), */ /* Re(x[2]), Im(x[2]), ... Re(x[N/2-1]), Im(x[N/2-1]). */ /* *Utbl = cosine table */ /* *BRLow = bit reversed counter table */ /* OUTPUTS */ /* *ioptr = real output data array */ SPFLOAT scale; int StageCnt; int NDiffU; scale = (SPFLOAT)(1.0 / (double)((int)POW2(M))); M = M - 1; switch (M) { case -1: break; case 0: rifft1pt(ioptr, scale); /* a 2 pt fft */ break; case 1: rifft2pt(ioptr, scale); /* a 4 pt fft */ break; case 2: rifft4pt(ioptr, scale); /* an 8 pt fft */ break; case 3: rifft8pt(ioptr, scale); /* a 16 pt fft */ break; default: ifrstage(ioptr, M + 1, Utbl); /* bit reverse and first radix 2 stage */ scbitrevR2(ioptr, M, BRLow, scale); StageCnt = (M - 1) / 3; /* number of radix 8 stages */ NDiffU = 2; /* one radix 2 stage already complete */ if ((M - 1 - (StageCnt * 3)) == 1) { ibfR2(ioptr, M, NDiffU); /* 1 radix 2 stage */ NDiffU *= 2; } if ((M - 1 - (StageCnt * 3)) == 2) { ibfR4(ioptr, M, NDiffU); /* 1 radix 4 stage */ NDiffU *= 4; } if (M <= (int) MCACHE) ibfstages(ioptr, M, Utbl, 2, NDiffU, StageCnt); /* RADIX 8 Stages */ else ifftrecurs(ioptr, M, Utbl, 2, NDiffU, StageCnt); /* RADIX 8 Stages */ } }