home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / C / FFT142.ZIP / IFFT.C < prev    next >
Encoding:
C/C++ Source or Header  |  1991-06-06  |  9.5 KB  |  385 lines
  1. /*    name...
  2.         ifft
  3.  
  4.     bugs...
  5.         Should do transforms on several functions (-b option).
  6.  
  7.     history...
  8.         -- 1.21 --
  9.         1 Jun 87    -a option suppresses abscissas on output.
  10.         -- 1.20 --
  11.         11 May 87    Using "float" rather than "double" variables
  12.         8 May 87    Reading into doubles and performing range check
  13.                     before converting to floats.
  14.         -- 1.1 --
  15.         6 Apr 87    Default abscissa step is 1 (rather than 0!)
  16.                     Some printing removed.
  17.         -- 1.0 --
  18.         20 Mar 87    default output to STDOUT
  19.         8 Feb 87    split off from fft program.
  20.  
  21.     performance...
  22.         128 point transform in 14.43 sec with 7.5 MHz V-20 and no 8087.
  23.         4096 point transform in 1:53.15 with 6 MHz 80286 and 4 MHz 80287.
  24. */
  25. #include <stdio.h>
  26. #include <math.h>
  27.  
  28. #define VERSION "1.21"
  29.  
  30.  
  31. #define BIGFLOAT 6.e38    /* no floats larger than this  */
  32. #define ENTRIES 4098
  33. #define MAXLABELS 50
  34. #define BUFSIZE 200
  35. #define FLOAT float    /* change to "double" to restore higher precision */
  36.  
  37. FILE *ifile=stdin, *ofile=stdout;
  38.  
  39. char buffer[128], iname[35], oname[35];
  40. char buf[BUFSIZE];
  41. char *label_array[MAXLABELS];
  42. char **p_text=label_array;
  43.  
  44.  
  45. int m,            /* number of input values to be retained */
  46. n,                /* number of data points to be Fourier transformed */
  47. *nnn,            /* points to n */
  48. breaking,        /* nonzero if finding separate transforms for several functions */
  49. labels,            /* number of labels stored */
  50. automatic_abscissas, /* nonzero if abscissas to be calculated */
  51. abscissa_arguments,
  52. keeping=0,        /* if nonzero, the number of input values to be kept */
  53. f_arguments=0,
  54. x_arguments=0,
  55. index_array[MAXLABELS],    /* indexes into x and y */
  56. *p_data=index_array;
  57.  
  58. FLOAT *x,        /* frequency values */
  59. *y;                /* magnitudes (on input) or voltages (on output) */
  60. double abscissa,
  61. origin,            /* abscissa value to be treated as time zero */
  62. abscissa_step=1.,
  63. fmin, fmax,
  64. xmin, xmax;        /* first & last data points to be used */
  65.  
  66. main(argc, argv) int argc; char **argv;
  67. {    int i;
  68.  
  69.     {double junk;
  70.     sscanf("3.4","%lf",&junk);        /* workaround for DESMET sscanf bug */
  71.     }
  72.  
  73.     argc = args(argc, argv);        /* fetch switches */
  74.  
  75.     read_data(argc, argv);
  76.     pad();
  77.  
  78. /*    fprintf(stderr, "\n\n   computing inverse fft\n\n");  */
  79.     ifft(&n, x, y);
  80.  
  81.     i=1;
  82.     while(1)
  83.         {if(!automatic_abscissas)  fprintf(ofile, "%15.6g", x[i] + origin);
  84.         fprintf(ofile, "%15.6g", y[i]);
  85.         if(++i>n) break;
  86.         fprintf(ofile, "\n");
  87.         }
  88.     fprintf(ofile, " \"\"\n");
  89.     if(ofile!=stdout)
  90.         {fclose(ofile);
  91.         }
  92. }
  93.  
  94. pad()
  95. {    int m;
  96.     double freq, df;
  97.     for (m=2; m<ENTRIES; m *= 2)
  98.         if (n<=m+1) break;
  99.     if(n==m+1) return;
  100. /*    printf("transforming %d points after padding\n", m+1);  */
  101.     freq=x[n];
  102.     df=freq-x[n-1];
  103.     while(n<=m)
  104.         {n++;
  105.         freq += df;
  106.         x[n]=freq;
  107.         y[2*n-1]=y[2*n]=0.;
  108.         }
  109. }
  110.  
  111. listt(y, i, j) FLOAT y[]; int i, j;    /* display y[i] through y[j]  */
  112. {    while (i<=j) {printf("%4d %15.9f \n", i, y[i]); ++i;}
  113. }
  114.  
  115.  
  116. read_data(argc, argv) int argc; char **argv;
  117. {    int i, j, length, nums;
  118.     double xx, yyr, yyi, d, *pd, sum;
  119.     char *s, *t, *stop;
  120.     char *malloc();
  121.     char *strchr();
  122.  
  123.     x=(FLOAT *)malloc(ENTRIES*sizeof(FLOAT));
  124.     y=(FLOAT *)malloc(2*ENTRIES*sizeof(FLOAT));
  125.     if(x==0 || y==0) {fprintf(stderr, "can\'t allocate buffer"); exit(1);}
  126.     argc--; argv++;
  127.     if(strchr(argv[0], '?')) help();
  128. /*            open input file         */
  129.     if(argc && *argv[0]!='-')
  130.         {ifile=fopen(argv[0], "r");
  131.         if(ifile==0) {fprintf(stderr, "file %s not found\n", argv[0]); exit(1);}
  132.         strcpy(iname, argv[0]);
  133.         argc--; argv++;
  134.         }
  135. /*            open output file         */
  136.     if(argc && *argv[0]!='-')
  137.         {strcpy(oname, argv[0]);
  138.         argc--; argv++;
  139.         unlink(oname);
  140.         if((ofile=fopen(oname, "w"))==0)
  141.             {fprintf(stderr, "can\'t open output file %s", oname);
  142.             exit(1);
  143.             }        
  144.         }
  145.  
  146.     fprint_cmd(ofile, "; ifft %s\n");
  147.  
  148.     p_data[0]=-1;
  149.     i=1;        /* note x[0] and y[0] aren't defined */
  150.     while(i<ENTRIES)
  151.         {if(fgets(buf, BUFSIZE, ifile)==0) break;
  152.         t=buf;
  153.         while(*t && isspace(*t)) t++;
  154.         if(*t == 0) continue;        /* skip blank lines */
  155.         buf[strlen(buf)-1]=0;        /* zero out the line feed */
  156.         if(buf[0]==';')                /* skip comment */
  157.             {fprintf(ofile, "%s\n", buf); continue;
  158.             }
  159.         if(t = strchr(buf, ';')) *t=0;    /* zap same-line comment */
  160.         if(automatic_abscissas)
  161.             {xx=abscissa;
  162.             abscissa+=abscissa_step;
  163.             sscanf(buf, "%lf %lf", &yyr, &yyi);
  164.             }
  165.         else
  166.             {sscanf(buf, "%lf %lf %lf", &xx, &yyr, &yyi);
  167.             }
  168.         range_check(xx); range_check(yyr); range_check(yyi);
  169.         x[i]=xx; y[i+i-1]=yyr; y[i+i]=yyi; /* convert doubles to floats here */
  170.         s=buf;                      /* start looking for label */
  171.         nums=3;
  172.         if(automatic_abscissas) nums--;
  173.         while(nums--)                    /* skip the numbers */
  174.             {while(*s==' ')s++;
  175.             while(*s && (*s!=' '))s++;
  176.             }
  177.         while(*s==' ')s++;
  178.         if((length=strlen(s))&&(labels<MAXLABELS))
  179.             {if(*s=='\"')           /* label in quotes */
  180.                 {t=s+1;
  181.                 while(*t && (*t!='\"')) t++;
  182.                 t++;
  183.                 }
  184.             else                    /* label without quotes */
  185.                 {t=s;
  186.                 while(*t && (*t!=' '))t++;
  187.                 }
  188.             *t=0;
  189.             length=t-s;
  190.             p_data[labels]=i;
  191.             p_text[labels]=(char *)malloc(length+1);
  192.             if(p_text[labels]) strcpy(p_text[labels++], s);
  193.             }
  194.         i++;
  195.         }
  196.     n=i-1;
  197.     if(breaking && (!labels || p_data[labels-1]!=n-1))
  198.         {p_data[labels]=i-1;
  199.         if(p_text[labels]=(char *)malloc(1)) *p_text[labels++]=0;
  200.         }
  201. }
  202.  
  203. /* check whether number is too big for a float */
  204. range_check(x) double x;
  205. {    if(fabs(x)>BIGFLOAT)
  206.         {printf("input number too large: %f\n", x);
  207.         exit(1);
  208.         }
  209. }
  210.  
  211. int streq(a,b) char *a,*b;
  212. {    while(*a)
  213.         {if(*a!=*b)return 0;
  214.         a++; b++;
  215.         }
  216.     return (*b==0);
  217. }
  218.  
  219. /* get_parameter - process one command line option
  220.         (return # parameters used) */
  221. get_parameter(argc, argv) int argc; char **argv;
  222. {    int i;
  223.     if(streq(*argv, "-a"))
  224.         {i=get_double(argc, argv, 1, &abscissa_step, &abscissa, &abscissa);
  225.         abscissa_arguments=i-1;
  226.         automatic_abscissas=1;
  227.         return i;
  228.         }
  229.     else if(streq(*argv, "-z"))
  230.         {origin=0.;
  231.         i=get_double(argc, argv, 1, &origin, &fmax, &fmax);
  232.         return i;
  233.         }
  234.  
  235.     else gripe(argv);
  236. }
  237.  
  238. char *message[]={
  239. "ifft   version ", VERSION,
  240. " - calculate inverse fast Fourier transform \n",
  241. "                     of frequency domain function\n",
  242. "usage: ifft  [infile  [outfile]]  [options]\n",
  243. "options are:\n",
  244. "  -a  [step]     automatic abscissas \n",
  245. "  -z  val        add <val> to each abscissa value\n",
  246. /*
  247. "  -b             break input after each label, \n",
  248. "                 find separate transforms\n",
  249. "  -f  min [max]  minimum and maximum frequencies\n",
  250. "  -t  min [max]  minimum and maximum times\n",
  251. */
  252. 0
  253. };
  254.  
  255.  
  256. /*    name...
  257.         fft
  258.  
  259.     purpose...
  260.         perform forward or inverse FFT
  261.  
  262.     history...
  263.         May 84  translated from FORTRAN
  264.         12 Jun 84  diagnostic printouts shortened
  265. */
  266.  
  267.     int ip, n2, i, j, k, nu, nn;
  268.     int k1n2, l, nu1;
  269.     double arg, p, q, s, c;
  270.     double ninv, time, dt, freq, df;
  271.     double xplus, xminus, yplus, yminus;
  272.     double uplus, uminus, vplus, vminus;
  273.  
  274. /*    format frequency domain data so its
  275.     inverse fast Fourier transform can be calculated
  276.  
  277.     on input:
  278.     n is one more than a power of 2:  n == 3, 9, 17, 33, 65, ...
  279.     x[1...n]  has  f(0)  f(1)  f(2)  ...  f(n-1)
  280.     y[1...2n] has  Qr(0) Qi(0)   Qr(1) Qi(1)  ...  Qr(n-1) Qi(n-1)
  281.  
  282.     on output:
  283.     n has been increased:  n = 2*n-1
  284.     x[1...n] has  t(0)  t(1)  t(2)  ...  t(n-1)
  285.     y[1...n] has  q(0)  q(1)  q(2)  ...  q(n-1)
  286. */
  287. ifft(nnn, x, y) int *nnn; FLOAT x[], y[];
  288. {
  289.     n= *nnn;
  290.     if(n<= 0) {fprintf(stderr, "ifft: bad value of n: %d", n); return 1;}
  291.     df= x[2]-x[1]; dt= .5/x[n];
  292.     if((fabs(df*dt*(n-1)-.5)>.01) || (fabs(x[n]-x[n-1]-df)>.01*df))
  293.         {fprintf(stderr, "ifft: frequencies not equally spaced");
  294.         exit(1);
  295.         }
  296.     nn= 2; nu= 0;
  297.     while(++nu<= 20)
  298.         {if(nn==n-1)break;
  299.         nn= nn+nn;
  300.         }
  301.     if(nu>20)
  302.         {fprintf(stderr, "ifft: n isn\'t 1 more than a power of 2");
  303.         return n;
  304.         }
  305.     ip= 2; i= 0;
  306.     while(++i<= n) {x[i]= y[ip]; y[i]= y[ip-1]; ip=ip+2;}
  307.     n2= n/2+1; ip= n; arg= 0.;
  308.     n= n-1;   /* original n, minus one */
  309.     ninv= 3.14159265358979/n;
  310.     i= 0;
  311.     while(++i<= n2)
  312.         {/* printf("i= %d ", i); */
  313.         s= sin(arg); c= cos(arg);
  314.         yplus= y[i]+y[ip]; yminus= y[i]-y[ip];
  315.         xplus= x[i]+x[ip]; xminus= x[i]-x[ip];
  316.         p= s*yminus+c*xplus; q= s*xplus-c*yminus;
  317.         x[i]= q-xminus; x[ip]= q+xminus;
  318.         y[i]= yplus-p;  y[ip]= yplus+p;
  319.         --ip;
  320.         arg= arg+ninv;   /* faster than: arg=ninv*i; */
  321.         }
  322.  
  323.     fft2(y, x, n, nu);
  324.  
  325. /*    fprintf(stderr, "ifft: transform calculated\n");
  326.     i=0; while(++i<=10) {printf("%6d %15.6f %15.6f\n", i, y[i], x[i]);}
  327. */
  328.     ip=n;
  329.     while(ip)
  330.         {y[ip+ip]= -df*x[ip]; y[ip+ip-1]=df*y[ip];
  331.         --ip;
  332.         }
  333.     time=0.; i=0;
  334.     *nnn=n=n+n;  /* twice (original n, minus 1) */
  335.     while(++i<=n) {x[i]=time; time=time+dt;}
  336. }
  337. fft2(xreal, ximag, n, nu)
  338. FLOAT xreal[], ximag[];
  339. int n, nu;
  340. {    double treal, timag;
  341.     FLOAT *xr1, *xr2, *xi1, *xi2;
  342.  
  343.     n2=n/2; nu1=nu-1; k=0; l=0;
  344.     while(++l<=nu)
  345.         {while(k<n)
  346.             {p=ibitr(k>>nu1, nu);
  347.             arg=3.14159265358979*2*p/n;
  348.             c=cos(arg); s=sin(arg); i=0;
  349.             while(++i<=n2)
  350.                 {++k; k1n2=k+n2;
  351.                 /* initialize four pointers */
  352.                 xr1=xreal+k; xr2=xreal+k1n2;
  353.                 xi1=ximag+k; xi2=ximag+k1n2;
  354.                 treal= *xr2*c+ *xi2*s;
  355.                 timag= *xi2*c- *xr2*s;
  356.                 *xr2= *xr1-treal;
  357.                 *xi2= *xi1-timag;
  358.                 *xr1= *xr1+treal;
  359.                 *xi1= *xi1+timag;
  360.                 }
  361.             k=k+n2;
  362.             }
  363.         k=0; --nu1; n2=n2/2;
  364.         }
  365.     k=0;
  366.     while(++k<=n)
  367.         {i=ibitr(k-1, nu)+1;
  368.         if(i>k)
  369.             {treal=xreal[k];   timag=ximag[k];
  370.             xreal[k]=xreal[i]; ximag[k]=ximag[i];
  371.             xreal[i]=treal;    ximag[i]=timag;
  372.             }
  373.         }
  374. }
  375.  
  376. /* reverse the last nu bits of j */
  377. ibitr(j, nu) int j, nu;
  378. {    int ib;
  379.     ib=0;
  380.     while(nu--){ib=(ib<<1)+(j&1); j=j>>1;}
  381.     return ib;
  382. }
  383.  
  384.  
  385.