diff options
Diffstat (limited to 'apps/plugins/pdbox/PDa/src/d_imayer_fft.c')
| -rw-r--r-- | apps/plugins/pdbox/PDa/src/d_imayer_fft.c | 1032 |
1 files changed, 1032 insertions, 0 deletions
diff --git a/apps/plugins/pdbox/PDa/src/d_imayer_fft.c b/apps/plugins/pdbox/PDa/src/d_imayer_fft.c new file mode 100644 index 0000000..d8e9e9f --- /dev/null +++ b/apps/plugins/pdbox/PDa/src/d_imayer_fft.c @@ -0,0 +1,1032 @@ +/* +** Algorithm: complex multiplication +** +** Performance: Bad accuracy, great speed. +** +** The simplest, way of generating trig values. Note, this method is +** not stable, and errors accumulate rapidly. +** +** Computation: 2 *, 1 + per value. +*/ + + +#include "m_fixed.h" + + +#define TRIG_FAST + +#ifdef TRIG_FAST +char mtrig_algorithm[] = "Simple"; +#define TRIG_VARS \ + t_fixed t_c,t_s; +#define TRIG_INIT(k,c,s) \ + { \ + t_c = fcostab[k]; \ + t_s = fsintab[k]; \ + c = itofix(1); \ + s = 0; \ + } +#define TRIG_NEXT(k,c,s) \ + { \ + t_fixed t = c; \ + c = mult(t,t_c) - mult(s,t_s); \ + s = mult(t,t_s) + mult(s,t_c); \ + } +#define TRIG_23(k,c1,s1,c2,s2,c3,s3) \ + { \ + c2 = mult(c1,c1) - mult(s1,s1); \ + s2 = (mult(c1,s1)<<2); \ + c3 = mult(c2,c1) - mult(s2,s1); \ + s3 = mult(c2,s1) + mult(s2,c1); \ + } +#endif +#define TRIG_RESET(k,c,s) + +/* +** Algorithm: O. Buneman's trig generator. +** +** Performance: Good accuracy, mediocre speed. +** +** Manipulates a log(n) table to stably create n evenly spaced trig +** values. The newly generated values lay halfway between two +** known values, and are calculated by appropriatelly scaling the +** average of the known trig values appropriatelly according to the +** angle between them. This process is inherently stable; and is +** about as accurate as most trig library functions. The errors +** caused by this code segment are primarily due to the less stable +** method to calculate values for 2t and s 3t. Note: I believe this +** algorithm is patented(!), see readme file for more details. +** +** Computation: 1 +, 1 *, much integer math, per trig value +** +*/ + +#ifdef TRIG_GOOD +#define TRIG_VARS \ + int t_lam=0; \ + double coswrk[TRIG_TABLE_SIZE],sinwrk[TRIG_TABLE_SIZE]; +#define TRIG_INIT(k,c,s) \ + { \ + int i; \ + for (i=0 ; i<=k ; i++) \ + {coswrk[i]=fcostab[i];sinwrk[i]=fsintab[i];} \ + t_lam = 0; \ + c = 1; \ + s = 0; \ + } +#define TRIG_NEXT(k,c,s) \ + { \ + int i,j; \ + (t_lam)++; \ + for (i=0 ; !((1<<i)&t_lam) ; i++); \ + i = k-i; \ + s = sinwrk[i]; \ + c = coswrk[i]; \ + if (i>1) \ + { \ + for (j=k-i+2 ; (1<<j)&t_lam ; j++); \ + j = k - j; \ + sinwrk[i] = halsec[i] * (sinwrk[i-1] + sinwrk[j]); \ + coswrk[i] = halsec[i] * (coswrk[i-1] + coswrk[j]); \ + } \ + } +#endif + + +#define TRIG_TAB_SIZE 22 + +static long long halsec[TRIG_TAB_SIZE]= {1,2,3}; + +#define FFTmult(x,y) mult(x,y) + + + + +#include "d_imayer_tables.h" + +#define SQRT2 ftofix(1.414213562373095048801688724209698) + + +void imayer_fht(t_fixed *fz, int n) +{ + int k,k1,k2,k3,k4,kx; + t_fixed *fi,*fn,*gi; + TRIG_VARS; + + for (k1=1,k2=0;k1<n;k1++) + { + t_fixed aa; + for (k=n>>1; (!((k2^=k)&k)); k>>=1); + if (k1>k2) + { + aa=fz[k1];fz[k1]=fz[k2];fz[k2]=aa; + } + } + for ( k=0 ; (1<<k)<n ; k++ ); + k &= 1; + if (k==0) + { + for (fi=fz,fn=fz+n;fi<fn;fi+=4) + { + t_fixed f0,f1,f2,f3; + f1 = fi[0 ]-fi[1 ]; + f0 = fi[0 ]+fi[1 ]; + f3 = fi[2 ]-fi[3 ]; + f2 = fi[2 ]+fi[3 ]; + fi[2 ] = (f0-f2); + fi[0 ] = (f0+f2); + fi[3 ] = (f1-f3); + fi[1 ] = (f1+f3); + } + } + else + { + for (fi=fz,fn=fz+n,gi=fi+1;fi<fn;fi+=8,gi+=8) + { + t_fixed bs1,bc1,bs2,bc2,bs3,bc3,bs4,bc4, + bg0,bf0,bf1,bg1,bf2,bg2,bf3,bg3; + bc1 = fi[0 ] - gi[0 ]; + bs1 = fi[0 ] + gi[0 ]; + bc2 = fi[2 ] - gi[2 ]; + bs2 = fi[2 ] + gi[2 ]; + bc3 = fi[4 ] - gi[4 ]; + bs3 = fi[4 ] + gi[4 ]; + bc4 = fi[6 ] - gi[6 ]; + bs4 = fi[6 ] + gi[6 ]; + bf1 = (bs1 - bs2); + bf0 = (bs1 + bs2); + bg1 = (bc1 - bc2); + bg0 = (bc1 + bc2); + bf3 = (bs3 - bs4); + bf2 = (bs3 + bs4); + bg3 = FFTmult(SQRT2,bc4); + bg2 = FFTmult(SQRT2,bc3); + fi[4 ] = bf0 - bf2; + fi[0 ] = bf0 + bf2; + fi[6 ] = bf1 - bf3; + fi[2 ] = bf1 + bf3; + gi[4 ] = bg0 - bg2; + gi[0 ] = bg0 + bg2; + gi[6 ] = bg1 - bg3; + gi[2 ] = bg1 + bg3; + } + } + if (n<16) return; + + do + { + t_fixed s1,c1; + int ii; + k += 2; + k1 = 1 << k; + k2 = k1 << 1; + k4 = k2 << 1; + k3 = k2 + k1; + kx = k1 >> 1; + fi = fz; + gi = fi + kx; + fn = fz + n; + do + { + t_fixed g0,f0,f1,g1,f2,g2,f3,g3; + f1 = fi[0 ] - fi[k1]; + f0 = fi[0 ] + fi[k1]; + f3 = fi[k2] - fi[k3]; + f2 = fi[k2] + fi[k3]; + fi[k2] = f0 - f2; + fi[0 ] = f0 + f2; + fi[k3] = f1 - f3; + fi[k1] = f1 + f3; + g1 = gi[0 ] - gi[k1]; + g0 = gi[0 ] + gi[k1]; + g3 = FFTmult(SQRT2, gi[k3]); + g2 = FFTmult(SQRT2, gi[k2]); + gi[k2] = g0 - g2; + gi[0 ] = g0 + g2; + gi[k3] = g1 - g3; + gi[k1] = g1 + g3; + gi += k4; + fi += k4; + } while (fi<fn); + TRIG_INIT(k,c1,s1); + for (ii=1;ii<kx;ii++) + { + t_fixed c2,s2; + TRIG_NEXT(k,c1,s1); + c2 = FFTmult(c1,c1) - FFTmult(s1,s1); + s2 = 2*FFTmult(c1,s1); + fn = fz + n; + fi = fz +ii; + gi = fz +k1-ii; + do + { + t_fixed a,b,g0,f0,f1,g1,f2,g2,f3,g3; + b = FFTmult(s2,fi[k1]) - FFTmult(c2,gi[k1]); + a = FFTmult(c2,fi[k1]) + FFTmult(s2,gi[k1]); + f1 = fi[0 ] - a; + f0 = fi[0 ] + a; + g1 = gi[0 ] - b; + g0 = gi[0 ] + b; + b = FFTmult(s2,fi[k3]) - FFTmult(c2,gi[k3]); + a = FFTmult(c2,fi[k3]) + FFTmult(s2,gi[k3]); + f3 = fi[k2] - a; + f2 = fi[k2] + a; + g3 = gi[k2] - b; + g2 = gi[k2] + b; + b = FFTmult(s1,f2) - FFTmult(c1,g3); + a = FFTmult(c1,f2) + FFTmult(s1,g3); + fi[k2] = f0 - a; + fi[0 ] = f0 + a; + gi[k3] = g1 - b; + gi[k1] = g1 + b; + b = FFTmult(c1,g2) - FFTmult(s1,f3); + a = FFTmult(s1,g2) + FFTmult(c1,f3); + gi[k2] = g0 - a; + gi[0 ] = g0 + a; + fi[k3] = f1 - b; + fi[k1] = f1 + b; + gi += k4; + fi += k4; + } while (fi<fn); + } + TRIG_RESET(k,c1,s1); + } while (k4<n); +} + + +void imayer_fft(int n, t_fixed *real, t_fixed *imag) +{ + t_fixed a,b,c,d; + t_fixed q,r,s,t; + int i,j,k; + for (i=1,j=n-1,k=n/2;i<k;i++,j--) { + a = real[i]; b = real[j]; q=a+b; r=a-b; + c = imag[i]; d = imag[j]; s=c+d; t=c-d; + real[i] = (q+t)>>1; real[j] = (q-t)>>1; + imag[i] = (s-r)>>1; imag[j] = (s+r)>>1; + } + imayer_fht(real,n); + imayer_fht(imag,n); +} + +void imayer_ifft(int n, t_fixed *real, t_fixed *imag) +{ + t_fixed a,b,c,d; + t_fixed q,r,s,t; + int i,j,k; + imayer_fht(real,n); + imayer_fht(imag,n); + for (i=1,j=n-1,k=n/2;i<k;i++,j--) { + a = real[i]; b = real[j]; q=a+b; r=a-b; + c = imag[i]; d = imag[j]; s=c+d; t=c-d; + imag[i] = (s+r)>>1; imag[j] = (s-r)>>1; + real[i] = (q-t)>>1; real[j] = (q+t)>>1; + } +} + +void imayer_realfft(int n, t_fixed *real) +{ + t_fixed a,b,c,d; + int i,j,k; + imayer_fht(real,n); + for (i=1,j=n-1,k=n/2;i<k;i++,j--) { + a = real[i]; + b = real[j]; + real[j] = (a-b)>>1; + real[i] = (a+b)>>1; + } +} + +void imayer_realifft(int n, t_fixed *real) +{ + t_fixed a,b,c,d; + int i,j,k; + for (i=1,j=n-1,k=n/2;i<k;i++,j--) { + a = real[i]; + b = real[j]; + real[j] = (a-b); + real[i] = (a+b); + } + imayer_fht(real,n); +} + + + +#ifdef TEST + +static double dfcostab[TRIG_TAB_SIZE]= + { + .00000000000000000000000000000000000000000000000000, + .70710678118654752440084436210484903928483593768847, + .92387953251128675612818318939678828682241662586364, + .98078528040323044912618223613423903697393373089333, + .99518472667219688624483695310947992157547486872985, + .99879545620517239271477160475910069444320361470461, + .99969881869620422011576564966617219685006108125772, + .99992470183914454092164649119638322435060646880221, + .99998117528260114265699043772856771617391725094433, + .99999529380957617151158012570011989955298763362218, + .99999882345170190992902571017152601904826792288976, + .99999970586288221916022821773876567711626389934930, + .99999992646571785114473148070738785694820115568892, + .99999998161642929380834691540290971450507605124278, + .99999999540410731289097193313960614895889430318945, + .99999999885102682756267330779455410840053741619428, + .99999999971275670684941397221864177608908945791828, + .99999999992818917670977509588385049026048033310951 + }; +static double dfsintab[TRIG_TAB_SIZE]= + { + 1.0000000000000000000000000000000000000000000000000, + .70710678118654752440084436210484903928483593768846, + .38268343236508977172845998403039886676134456248561, + .19509032201612826784828486847702224092769161775195, + .09801714032956060199419556388864184586113667316749, + .04906767432741801425495497694268265831474536302574, + .02454122852291228803173452945928292506546611923944, + .01227153828571992607940826195100321214037231959176, + .00613588464915447535964023459037258091705788631738, + .00306795676296597627014536549091984251894461021344, + .00153398018628476561230369715026407907995486457522, + .00076699031874270452693856835794857664314091945205, + .00038349518757139558907246168118138126339502603495, + .00019174759731070330743990956198900093346887403385, + .00009587379909597734587051721097647635118706561284, + .00004793689960306688454900399049465887274686668768, + .00002396844980841821872918657716502182009476147489, + .00001198422490506970642152156159698898480473197753 + }; + +static double dhalsec[TRIG_TAB_SIZE]= + { + 0, + 0, + .54119610014619698439972320536638942006107206337801, + .50979557910415916894193980398784391368261849190893, + .50241928618815570551167011928012092247859337193963, + .50060299823519630134550410676638239611758632599591, + .50015063602065098821477101271097658495974913010340, + .50003765191554772296778139077905492847503165398345, + .50000941253588775676512870469186533538523133757983, + .50000235310628608051401267171204408939326297376426, + .50000058827484117879868526730916804925780637276181, + .50000014706860214875463798283871198206179118093251, + .50000003676714377807315864400643020315103490883972, + .50000000919178552207366560348853455333939112569380, + .50000000229794635411562887767906868558991922348920, + .50000000057448658687873302235147272458812263401372, + .50000000014362164661654736863252589967935073278768, + .50000000003590541164769084922906986545517021050714 + }; + + +#include <stdio.h> + + +int writetables() +{ + int i; + + printf("/* Tables for fixed point lookup with %d bit precision*/\n\n",fix1); + + printf("static int fsintab[TRIG_TAB_SIZE]= {\n"); + + for (i=0;i<TRIG_TAB_SIZE-1;i++) + printf("%ld, \n",ftofix(dfsintab[i])); + printf("%ld }; \n",ftofix(dfsintab[i])); + + + printf("\n\nstatic int fcostab[TRIG_TAB_SIZE]= {\n"); + + for (i=0;i<TRIG_TAB_SIZE-1;i++) + printf("%ld, \n",ftofix(dfcostab[i])); + printf("%ld }; \n",ftofix(dfcostab[i])); + +} + + +#define OUTPUT \ + fprintf(stderr,"Integer - Float\n");\ + for (i=0;i<NP;i++)\ + fprintf(stderr,"%f %f --> %f %f\n",fixtof(r[i]),fixtof(im[i]),\ + fr[i],fim[i]);\ + fprintf(stderr,"\n\n"); + + + +int main() +{ + int i; + t_fixed r[256]; + t_fixed im[256]; + float fr[256]; + float fim[256]; + + +#if 1 + writetables(); + exit(0); +#endif + + +#if 0 + t_fixed c1,s1,c2,s2,c3,s3; + int k; + int i; + + TRIG_VARS; + + for (k=2;k<10;k+=2) { + TRIG_INIT(k,c1,s1); + for (i=0;i<8;i++) { + TRIG_NEXT(k,c1,s1); + TRIG_23(k,c1,s1,c2,s2,c3,s3); + printf("TRIG value k=%d,%d val1 = %f %f\n",k,i,fixtof(c1),fixtof(s1)); + } + } +#endif + + + +#if 1 + + #define NP 16 + + for (i=0;i<NP;i++) { + fr[i] = 0.; + r[i] = 0.; + fim[i] = 0.; + im[i] = 0; + } + +#if 0 + for (i=0;i<NP;i++) { + if (i&1) { + fr[i] = 0.1*i; + r[i] = ftofix(0.1*i); + } + else { + fr[i] = 0.; + r[i] = 0.; + } + } +#endif +#if 0 + for (i=0;i<NP;i++) { + fr[i] = 0.1; + r[i] = ftofix(0.1); + } +#endif + + r[1] = ftofix(0.1); + fr[1] = 0.1; + + + + /* TEST RUN */ + + OUTPUT + + imayer_fft(NP,r,im); + mayer_fft(NP,fr,fim); +// imayer_fht(r,NP); +// mayer_fht(fr,NP); + +#if 1 + for (i=0;i<NP;i++) { + r[i] = mult(ftofix(0.01),r[i]); + fr[i] = 0.01*fr[i]; + } +#endif + + OUTPUT + + + imayer_fft(NP,r,im); + mayer_fft(NP,fr,fim); + + OUTPUT + + +#endif + + +} + +#endif +/* +** Algorithm: complex multiplication +** +** Performance: Bad accuracy, great speed. +** +** The simplest, way of generating trig values. Note, this method is +** not stable, and errors accumulate rapidly. +** +** Computation: 2 *, 1 + per value. +*/ + + +#include "m_fixed.h" + + +#define TRIG_FAST + +#ifdef TRIG_FAST +char mtrig_algorithm[] = "Simple"; +#define TRIG_VARS \ + t_fixed t_c,t_s; +#define TRIG_INIT(k,c,s) \ + { \ + t_c = fcostab[k]; \ + t_s = fsintab[k]; \ + c = itofix(1); \ + s = 0; \ + } +#define TRIG_NEXT(k,c,s) \ + { \ + t_fixed t = c; \ + c = mult(t,t_c) - mult(s,t_s); \ + s = mult(t,t_s) + mult(s,t_c); \ + } +#define TRIG_23(k,c1,s1,c2,s2,c3,s3) \ + { \ + c2 = mult(c1,c1) - mult(s1,s1); \ + s2 = (mult(c1,s1)<<2); \ + c3 = mult(c2,c1) - mult(s2,s1); \ + s3 = mult(c2,s1) + mult(s2,c1); \ + } +#endif +#define TRIG_RESET(k,c,s) + +/* +** Algorithm: O. Buneman's trig generator. +** +** Performance: Good accuracy, mediocre speed. +** +** Manipulates a log(n) table to stably create n evenly spaced trig +** values. The newly generated values lay halfway between two +** known values, and are calculated by appropriatelly scaling the +** average of the known trig values appropriatelly according to the +** angle between them. This process is inherently stable; and is +** about as accurate as most trig library functions. The errors +** caused by this code segment are primarily due to the less stable +** method to calculate values for 2t and s 3t. Note: I believe this +** algorithm is patented(!), see readme file for more details. +** +** Computation: 1 +, 1 *, much integer math, per trig value +** +*/ + +#ifdef TRIG_GOOD +#define TRIG_VARS \ + int t_lam=0; \ + double coswrk[TRIG_TABLE_SIZE],sinwrk[TRIG_TABLE_SIZE]; +#define TRIG_INIT(k,c,s) \ + { \ + int i; \ + for (i=0 ; i<=k ; i++) \ + {coswrk[i]=fcostab[i];sinwrk[i]=fsintab[i];} \ + t_lam = 0; \ + c = 1; \ + s = 0; \ + } +#define TRIG_NEXT(k,c,s) \ + { \ + int i,j; \ + (t_lam)++; \ + for (i=0 ; !((1<<i)&t_lam) ; i++); \ + i = k-i; \ + s = sinwrk[i]; \ + c = coswrk[i]; \ + if (i>1) \ + { \ + for (j=k-i+2 ; (1<<j)&t_lam ; j++); \ + j = k - j; \ + sinwrk[i] = halsec[i] * (sinwrk[i-1] + sinwrk[j]); \ + coswrk[i] = halsec[i] * (coswrk[i-1] + coswrk[j]); \ + } \ + } +#endif + + +#define TRIG_TAB_SIZE 22 + +static long long halsec[TRIG_TAB_SIZE]= {1,2,3}; + +#define FFTmult(x,y) mult(x,y) + + + + +#include "d_imayer_tables.h" + +#define SQRT2 ftofix(1.414213562373095048801688724209698) + + +void imayer_fht(t_fixed *fz, int n) +{ + int k,k1,k2,k3,k4,kx; + t_fixed *fi,*fn,*gi; + TRIG_VARS; + + for (k1=1,k2=0;k1<n;k1++) + { + t_fixed aa; + for (k=n>>1; (!((k2^=k)&k)); k>>=1); + if (k1>k2) + { + aa=fz[k1];fz[k1]=fz[k2];fz[k2]=aa; + } + } + for ( k=0 ; (1<<k)<n ; k++ ); + k &= 1; + if (k==0) + { + for (fi=fz,fn=fz+n;fi<fn;fi+=4) + { + t_fixed f0,f1,f2,f3; + f1 = fi[0 ]-fi[1 ]; + f0 = fi[0 ]+fi[1 ]; + f3 = fi[2 ]-fi[3 ]; + f2 = fi[2 ]+fi[3 ]; + fi[2 ] = (f0-f2); + fi[0 ] = (f0+f2); + fi[3 ] = (f1-f3); + fi[1 ] = (f1+f3); + } + } + else + { + for (fi=fz,fn=fz+n,gi=fi+1;fi<fn;fi+=8,gi+=8) + { + t_fixed bs1,bc1,bs2,bc2,bs3,bc3,bs4,bc4, + bg0,bf0,bf1,bg1,bf2,bg2,bf3,bg3; + bc1 = fi[0 ] - gi[0 ]; + bs1 = fi[0 ] + gi[0 ]; + bc2 = fi[2 ] - gi[2 ]; + bs2 = fi[2 ] + gi[2 ]; + bc3 = fi[4 ] - gi[4 ]; + bs3 = fi[4 ] + gi[4 ]; + bc4 = fi[6 ] - gi[6 ]; + bs4 = fi[6 ] + gi[6 ]; + bf1 = (bs1 - bs2); + bf0 = (bs1 + bs2); + bg1 = (bc1 - bc2); + bg0 = (bc1 + bc2); + bf3 = (bs3 - bs4); + bf2 = (bs3 + bs4); + bg3 = FFTmult(SQRT2,bc4); + bg2 = FFTmult(SQRT2,bc3); + fi[4 ] = bf0 - bf2; + fi[0 ] = bf0 + bf2; + fi[6 ] = bf1 - bf3; + fi[2 ] = bf1 + bf3; + gi[4 ] = bg0 - bg2; + gi[0 ] = bg0 + bg2; + gi[6 ] = bg1 - bg3; + gi[2 ] = bg1 + bg3; + } + } + if (n<16) return; + + do + { + t_fixed s1,c1; + int ii; + k += 2; + k1 = 1 << k; + k2 = k1 << 1; + k4 = k2 << 1; + k3 = k2 + k1; + kx = k1 >> 1; + fi = fz; + gi = fi + kx; + fn = fz + n; + do + { + t_fixed g0,f0,f1,g1,f2,g2,f3,g3; + f1 = fi[0 ] - fi[k1]; + f0 = fi[0 ] + fi[k1]; + f3 = fi[k2] - fi[k3]; + f2 = fi[k2] + fi[k3]; + fi[k2] = f0 - f2; + fi[0 ] = f0 + f2; + fi[k3] = f1 - f3; + fi[k1] = f1 + f3; + g1 = gi[0 ] - gi[k1]; + g0 = gi[0 ] + gi[k1]; + g3 = FFTmult(SQRT2, gi[k3]); + g2 = FFTmult(SQRT2, gi[k2]); + gi[k2] = g0 - g2; + gi[0 ] = g0 + g2; + gi[k3] = g1 - g3; + gi[k1] = g1 + g3; + gi += k4; + fi += k4; + } while (fi<fn); + TRIG_INIT(k,c1,s1); + for (ii=1;ii<kx;ii++) + { + t_fixed c2,s2; + TRIG_NEXT(k,c1,s1); + c2 = FFTmult(c1,c1) - FFTmult(s1,s1); + s2 = 2*FFTmult(c1,s1); + fn = fz + n; + fi = fz +ii; + gi = fz +k1-ii; + do + { + t_fixed a,b,g0,f0,f1,g1,f2,g2,f3,g3; + b = FFTmult(s2,fi[k1]) - FFTmult(c2,gi[k1]); + a = FFTmult(c2,fi[k1]) + FFTmult(s2,gi[k1]); + f1 = fi[0 ] - a; + f0 = fi[0 ] + a; + g1 = gi[0 ] - b; + g0 = gi[0 ] + b; + b = FFTmult(s2,fi[k3]) - FFTmult(c2,gi[k3]); + a = FFTmult(c2,fi[k3]) + FFTmult(s2,gi[k3]); + f3 = fi[k2] - a; + f2 = fi[k2] + a; + g3 = gi[k2] - b; + g2 = gi[k2] + b; + b = FFTmult(s1,f2) - FFTmult(c1,g3); + a = FFTmult(c1,f2) + FFTmult(s1,g3); + fi[k2] = f0 - a; + fi[0 ] = f0 + a; + gi[k3] = g1 - b; + gi[k1] = g1 + b; + b = FFTmult(c1,g2) - FFTmult(s1,f3); + a = FFTmult(s1,g2) + FFTmult(c1,f3); + gi[k2] = g0 - a; + gi[0 ] = g0 + a; + fi[k3] = f1 - b; + fi[k1] = f1 + b; + gi += k4; + fi += k4; + } while (fi<fn); + } + TRIG_RESET(k,c1,s1); + } while (k4<n); +} + + +void imayer_fft(int n, t_fixed *real, t_fixed *imag) +{ + t_fixed a,b,c,d; + t_fixed q,r,s,t; + int i,j,k; + for (i=1,j=n-1,k=n/2;i<k;i++,j--) { + a = real[i]; b = real[j]; q=a+b; r=a-b; + c = imag[i]; d = imag[j]; s=c+d; t=c-d; + real[i] = (q+t)>>1; real[j] = (q-t)>>1; + imag[i] = (s-r)>>1; imag[j] = (s+r)>>1; + } + imayer_fht(real,n); + imayer_fht(imag,n); +} + +void imayer_ifft(int n, t_fixed *real, t_fixed *imag) +{ + t_fixed a,b,c,d; + t_fixed q,r,s,t; + int i,j,k; + imayer_fht(real,n); + imayer_fht(imag,n); + for (i=1,j=n-1,k=n/2;i<k;i++,j--) { + a = real[i]; b = real[j]; q=a+b; r=a-b; + c = imag[i]; d = imag[j]; s=c+d; t=c-d; + imag[i] = (s+r)>>1; imag[j] = (s-r)>>1; + real[i] = (q-t)>>1; real[j] = (q+t)>>1; + } +} + +void imayer_realfft(int n, t_fixed *real) +{ + t_fixed a,b,c,d; + int i,j,k; + imayer_fht(real,n); + for (i=1,j=n-1,k=n/2;i<k;i++,j--) { + a = real[i]; + b = real[j]; + real[j] = (a-b)>>1; + real[i] = (a+b)>>1; + } +} + +void imayer_realifft(int n, t_fixed *real) +{ + t_fixed a,b,c,d; + int i,j,k; + for (i=1,j=n-1,k=n/2;i<k;i++,j--) { + a = real[i]; + b = real[j]; + real[j] = (a-b); + real[i] = (a+b); + } + imayer_fht(real,n); +} + + + +#ifdef TEST + +static double dfcostab[TRIG_TAB_SIZE]= + { + .00000000000000000000000000000000000000000000000000, + .70710678118654752440084436210484903928483593768847, + .92387953251128675612818318939678828682241662586364, + .98078528040323044912618223613423903697393373089333, + .99518472667219688624483695310947992157547486872985, + .99879545620517239271477160475910069444320361470461, + .99969881869620422011576564966617219685006108125772, + .99992470183914454092164649119638322435060646880221, + .99998117528260114265699043772856771617391725094433, + .99999529380957617151158012570011989955298763362218, + .99999882345170190992902571017152601904826792288976, + .99999970586288221916022821773876567711626389934930, + .99999992646571785114473148070738785694820115568892, + .99999998161642929380834691540290971450507605124278, + .99999999540410731289097193313960614895889430318945, + .99999999885102682756267330779455410840053741619428, + .99999999971275670684941397221864177608908945791828, + .99999999992818917670977509588385049026048033310951 + }; +static double dfsintab[TRIG_TAB_SIZE]= + { + 1.0000000000000000000000000000000000000000000000000, + .70710678118654752440084436210484903928483593768846, + .38268343236508977172845998403039886676134456248561, + .19509032201612826784828486847702224092769161775195, + .09801714032956060199419556388864184586113667316749, + .04906767432741801425495497694268265831474536302574, + .02454122852291228803173452945928292506546611923944, + .01227153828571992607940826195100321214037231959176, + .00613588464915447535964023459037258091705788631738, + .00306795676296597627014536549091984251894461021344, + .00153398018628476561230369715026407907995486457522, + .00076699031874270452693856835794857664314091945205, + .00038349518757139558907246168118138126339502603495, + .00019174759731070330743990956198900093346887403385, + .00009587379909597734587051721097647635118706561284, + .00004793689960306688454900399049465887274686668768, + .00002396844980841821872918657716502182009476147489, + .00001198422490506970642152156159698898480473197753 + }; + +static double dhalsec[TRIG_TAB_SIZE]= + { + 0, + 0, + .54119610014619698439972320536638942006107206337801, + .50979557910415916894193980398784391368261849190893, + .50241928618815570551167011928012092247859337193963, + .50060299823519630134550410676638239611758632599591, + .50015063602065098821477101271097658495974913010340, + .50003765191554772296778139077905492847503165398345, + .50000941253588775676512870469186533538523133757983, + .50000235310628608051401267171204408939326297376426, + .50000058827484117879868526730916804925780637276181, + .50000014706860214875463798283871198206179118093251, + .50000003676714377807315864400643020315103490883972, + .50000000919178552207366560348853455333939112569380, + .50000000229794635411562887767906868558991922348920, + .50000000057448658687873302235147272458812263401372, + .50000000014362164661654736863252589967935073278768, + .50000000003590541164769084922906986545517021050714 + }; + + +#include <stdio.h> + + +int writetables() +{ + int i; + + printf("/* Tables for fixed point lookup with %d bit precision*/\n\n",fix1); + + printf("static int fsintab[TRIG_TAB_SIZE]= {\n"); + + for (i=0;i<TRIG_TAB_SIZE-1;i++) + printf("%ld, \n",ftofix(dfsintab[i])); + printf("%ld }; \n",ftofix(dfsintab[i])); + + + printf("\n\nstatic int fcostab[TRIG_TAB_SIZE]= {\n"); + + for (i=0;i<TRIG_TAB_SIZE-1;i++) + printf("%ld, \n",ftofix(dfcostab[i])); + printf("%ld }; \n",ftofix(dfcostab[i])); + +} + + +#define OUTPUT \ + fprintf(stderr,"Integer - Float\n");\ + for (i=0;i<NP;i++)\ + fprintf(stderr,"%f %f --> %f %f\n",fixtof(r[i]),fixtof(im[i]),\ + fr[i],fim[i]);\ + fprintf(stderr,"\n\n"); + + + +int main() +{ + int i; + t_fixed r[256]; + t_fixed im[256]; + float fr[256]; + float fim[256]; + + +#if 1 + writetables(); + exit(0); +#endif + + +#if 0 + t_fixed c1,s1,c2,s2,c3,s3; + int k; + int i; + + TRIG_VARS; + + for (k=2;k<10;k+=2) { + TRIG_INIT(k,c1,s1); + for (i=0;i<8;i++) { + TRIG_NEXT(k,c1,s1); + TRIG_23(k,c1,s1,c2,s2,c3,s3); + printf("TRIG value k=%d,%d val1 = %f %f\n",k,i,fixtof(c1),fixtof(s1)); + } + } +#endif + + + +#if 1 + + #define NP 16 + + for (i=0;i<NP;i++) { + fr[i] = 0.; + r[i] = 0.; + fim[i] = 0.; + im[i] = 0; + } + +#if 0 + for (i=0;i<NP;i++) { + if (i&1) { + fr[i] = 0.1*i; + r[i] = ftofix(0.1*i); + } + else { + fr[i] = 0.; + r[i] = 0.; + } + } +#endif +#if 0 + for (i=0;i<NP;i++) { + fr[i] = 0.1; + r[i] = ftofix(0.1); + } +#endif + + r[1] = ftofix(0.1); + fr[1] = 0.1; + + + + /* TEST RUN */ + + OUTPUT + + imayer_fft(NP,r,im); + mayer_fft(NP,fr,fim); +// imayer_fht(r,NP); +// mayer_fht(fr,NP); + +#if 1 + for (i=0;i<NP;i++) { + r[i] = mult(ftofix(0.01),r[i]); + fr[i] = 0.01*fr[i]; + } +#endif + + OUTPUT + + + imayer_fft(NP,r,im); + mayer_fft(NP,fr,fim); + + OUTPUT + + +#endif + + +} + +#endif |