home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_300 / 315_01 / fft.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-05-14  |  13.5 KB  |  566 lines
  1.  
  2. /* int_scale(), dbl_scale(), and long_scale() moved from herc.c */
  3. /* 10/89, because herc.c is being discontinued, and msherc.com is */
  4. /* being used for hercules graphcis support.  T Clune */
  5.  
  6. /* fft.c is a file of functions to support fourier transforms. */
  7. /* Written 3/89 by T Clune.  Copyright (c) 1989, Eye Research Institute, */
  8. /* 20 Staniford St, Boston, MA 02114.  All Rights Reserved. */
  9.  
  10. #include <stdio.h>
  11. #include <stdlib.h>
  12. #include <math.h>
  13. #include <malloc.h>
  14. #include <conio.h>
  15. #include "fft.h"
  16.  
  17.  
  18.  
  19. /* data_multiply() multiplies two data files of equal length together */
  20.  
  21. void data_multiply(file1, file2)
  22. char file1[], file2[];
  23. {
  24.     FILE *f1, *f2, *f3;
  25.     char string[80];
  26.     int i, numpts1, numpts2;
  27.     double d1, d2, *d_out, minval, maxval;
  28.  
  29.  
  30.     f1=fopen(file1,"r");
  31.     if(f1==NULL)
  32.     {
  33.     printf("Error opening data file.  Exiting routine.\n");
  34.     printf("Press any key to continue...\n");
  35.     mgetch();
  36.     return;
  37.     }
  38.  
  39.     f2=fopen(file2,"r");
  40.     if(f2==NULL)
  41.     {
  42.     printf("Error opening data file.  Exiting routine.\n");
  43.     printf("Press any key to continue...\n");
  44.     mgetch();
  45.     return;
  46.     }
  47.  
  48.     fscanf(f1,"%d %*lf %*lf", &numpts1);
  49.     fscanf(f2,"%d %*lf %*lf", &numpts2);
  50.     if(numpts1 != numpts2)
  51.     {
  52.     printf("Different size input files.  Exiting routine.\n");
  53.     printf("Press any key to continue...\n");
  54.     mgetch();
  55.     return;
  56.     }
  57.  
  58.     while(kbhit()) ;
  59.     printf("Enter filespec for output file\n");
  60.     mgets(string);
  61.     f3=fopen(string,"w");
  62.  
  63.     if(f3==NULL)
  64.     {
  65.     printf("Error opening output file.  Exiting routine.\n");
  66.     printf("Press any key to continue...\n");
  67.     mgetch();
  68.     return;
  69.     }
  70.     d_out=(double *)malloc(sizeof(double)*numpts1);
  71.     if(d_out==NULL)
  72.     {
  73.     printf("malloc() error in data_multiply(). Program aborting.\n");
  74.     exit(-1);
  75.     }
  76.  
  77.  
  78.     for(i=0;i<numpts1;i++)
  79.     {
  80.     fscanf(f1,"%lf", &d1);
  81.     fscanf(f2,"%lf",&d2);
  82.     d_out[i]=d1*d2;
  83.     }
  84.  
  85.     minmax(d_out,numpts1,&minval,&maxval);
  86.     fprintf(f3,"%d %f %f\n", numpts1, minval, maxval);
  87.  
  88.     for(i=0;i<numpts1;i++)
  89.     fprintf(f3,"%f\n",d_out[i]);
  90.  
  91.     fclose(f1);
  92.     fclose(f2);
  93.     fclose(f3);
  94.     free(d_out);
  95.  
  96. }
  97.  
  98.  
  99.  
  100.  
  101. /* dbl_scale() scales an extended precision floating point number */
  102. /* to an integer value.  It is useful for scaling floating point  */
  103. /* numbers for graphing on the screen.                            */
  104. /* moved from herc.c 10/89 by T Clune */
  105.  
  106.  
  107. int dbl_scale(val, min, max, out_min, out_max)
  108.  
  109. double val, min, max;  /* val= number to be scaled; min= smallest possible */
  110.     /* value for val; max=largest possible value */
  111. int out_min, out_max;  /* out_min=integer to return if val=min; */
  112.         /* out_max=integer to return if val=max.  NOTE THAT out_min */
  113.         /* MAY BE LARGER THAN out_max.  If you want to change screen */
  114.         /* values from 4th quadrant to 1st quadrant, the y-value */
  115.         /* out_min would be larger than the y-value out_max, while */
  116.         /* x-value out_max would be larger than the x-value out_min */
  117.  
  118. {
  119.  
  120.     int out_val, largest, smallest; /* val returned, larger of out_min */
  121.         /* and out_max, smaller of the same */
  122.  
  123.     double scale_factor;    /* ratio of in and out ranges */
  124.         /* scale the input value */
  125.  
  126.     scale_factor = (out_max - out_min) / (max - min);
  127.     out_val = (int)((val - min) * scale_factor + out_min);
  128.  
  129.  
  130.     if(out_min>out_max)     /* check for max/min inversion */
  131.     {
  132.     smallest = out_max;
  133.     largest = out_min;
  134.  
  135.     }
  136.     else
  137.     {
  138.     smallest = out_min;
  139.     largest = out_max;
  140.     }
  141.  
  142.     if(out_val > largest)      /* check for out-of-bounds values */
  143.     out_val = largest;     /* since val is not checked to insure */
  144.     if(out_val < smallest)      /* that it is between min and max, this */
  145.     out_val = smallest;     /* is necessary.  By checking at the last */
  146.                 /* step, any arithmetic round-off accumulation */
  147.                 /* is also trapped */
  148.  
  149.     return out_val;          /* return the scaled value */
  150.  
  151. }
  152.  
  153.  
  154.  
  155.  
  156.  
  157. /* fft() accepts the time-domain data x,y (n points, where n is a power of 2) */
  158. /* and converts to frequency domain rectangular data if flag==FORWARD, */
  159. /* or accepts frequency-domain data in x,y and converts to time domain if */
  160. /* flag==INVERSE.  FORWARD and INVERSE are defined in fft.h */
  161.  
  162. void fft(x, y, n, flag)
  163. double x[], y[];
  164. int n, flag;
  165. {
  166.    int maxpower, arg, pnt0, pnt1;
  167.    int i, j=0, a, b, k;
  168.    double prodreal, prodimag, harm, temp;
  169.    double *cos_coef;
  170.    double *sin_coef;
  171.  
  172.    cos_coef=(double*)calloc(n,sizeof(double));
  173.    sin_coef=(double*)calloc(n,sizeof(double));
  174.  
  175.    if((cos_coef==NULL) || (sin_coef==NULL))
  176.    {
  177.     printf("malloc() error in fft function\n");
  178.     if(cos_coef != NULL)
  179.         free(cos_coef);
  180.     if(sin_coef != NULL)
  181.         free(sin_coef);
  182.     return;
  183.     }
  184.  
  185.    if(flag==INVERSE)
  186.    {
  187.       for (i=0;i<n;i++)
  188.       {
  189.      x[i] /=n;
  190.      y[i] /=n;
  191.       }
  192.    }
  193.  
  194.    for(i=0;i<(n-1);i++)
  195.    {
  196.       if(i<j)
  197.       {
  198.      temp=x[i];
  199.      x[i]=x[j];
  200.      x[j]=temp;
  201.      temp=y[i];
  202.      y[i]=y[j];
  203.      y[j]=temp;
  204.       }
  205.       k=n/2;
  206.       while(k<=j)
  207.       {
  208.      j -=k;
  209.      k /=2;
  210.       }
  211.       j +=k;
  212.    }
  213.    maxpower=0;
  214.    i=n;
  215.    while(i>1)
  216.    {
  217.       maxpower++;
  218.       i /=2;
  219.    }
  220.    harm=2*PI/n;
  221.    for(i=0;i<n;i++)
  222.    {
  223.       sin_coef[i]=flag*sin(harm*i);
  224.       cos_coef[i]=cos(harm*i);
  225.    }
  226.    a=2;
  227.    b=1;
  228.  
  229.    for(j=0;j<maxpower;j++)
  230.    {
  231.       pnt0=n/a;
  232.       pnt1=0;
  233.       for(k=0;k<b;k++)
  234.       {
  235.      i=k;
  236.      while(i<n)
  237.      {
  238.         arg=i+b;
  239.         if(k==0)
  240.         {
  241.            prodreal=x[arg];
  242.            prodimag=y[arg];
  243.         }
  244.         else
  245.         {
  246.            prodreal=x[arg]*cos_coef[pnt1]-y[arg]*sin_coef[pnt1];
  247.            prodimag=x[arg]*sin_coef[pnt1]+y[arg]*cos_coef[pnt1];
  248.         }
  249.         x[arg]=x[i]-prodreal;
  250.         y[arg]=y[i]-prodimag;
  251.         x[i] +=prodreal;
  252.         y[i] +=prodimag;
  253.         i +=a;
  254.      }
  255.      pnt1 +=pnt0;
  256.       }
  257.       a *=2;
  258.       b *=2;
  259.    }
  260.    free(cos_coef);
  261.    free(sin_coef);
  262. }
  263.  
  264.  
  265.  
  266.  
  267.  
  268. /* get_filter() builds a high_pass or low_pass filter.  The arguments are: */
  269. /* FILTER is the array for the filter coefficients (frequency-domain */
  270. /* AMPLITUDE coefficients [0 to 1]), N is the number of elements in FILTER */
  271. /* SIGN is -1 for lowpass, +1 for highpass; NYQUIST is the nyquist freq */
  272. /* for the FILTER data set, CUTOFF is the filter cutoff frequency */
  273. /* ROLLOFF is the rolloff in dBells per filter_unit; */
  274. /* FILTER_UNIT is 2.0 for rolloff as dB/octave or 10.0 for dB/decade; */
  275. /* get_filter() will create positive/negative format amplitude filter */
  276. /* coefficients ONLY. The function returns the coefficients in FILTER */
  277.  
  278. double * get_filter(filter, n, sign, nyquist, cutoff, rolloff, filter_units)
  279. double filter[];
  280. int n, sign;
  281. double nyquist, cutoff, rolloff, filter_units;
  282. {
  283.     double freq_interval, active_freq;
  284.     double half_power, exp_var, octaves;
  285.     double filter_factor;
  286.     int i,m;
  287.  
  288.     half_power=sqrt(0.5);
  289.     filter_factor = log10(filter_units);
  290.     m=n/2;
  291.     freq_interval=nyquist/m;
  292.  
  293.     if(sign== (-1))
  294.     filter[0]=1.0;
  295.     else
  296.     filter[0]=0.0;
  297.  
  298.     for(i=1;i<=m;i++)
  299.     {
  300.     active_freq=freq_interval*i;
  301.  
  302.         /* number of octaves: */
  303.     octaves=log10(active_freq/cutoff)/filter_factor;
  304.         /* 20=constant for dB formula as amplitude */
  305.     exp_var=rolloff*octaves*sign/20.0;
  306.     exp_var=pow(10.0,exp_var);
  307.     filter[i]=exp_var*half_power;
  308.     if(filter[i]>1.0)
  309.         filter[i]=1.0;
  310.     if(filter[i]<MINVAL)
  311.     filter[i]=0.0;
  312.     if(sign==(-1) && filter[i]==0.0)
  313.         for( ;i<=m;i++)
  314.         filter[i]=0.0;
  315.     if(sign==1 && filter[i]==1.0)
  316.         for( ;i<=m;i++)
  317.         filter[i]=1.0;
  318.     }
  319.     for(i=1;i<m;i++)
  320.     filter[m+i]=filter[m-i];
  321.     unscramble_transform(filter,n);
  322.  
  323.     return filter;
  324.  
  325. }
  326.  
  327.  
  328.  
  329.  
  330.  
  331. /* int_scale() is a routine to scale integer values to fit graph */
  332. /* moved from herc.c 10/89 by T Clune */
  333.  
  334. int int_scale(in_val, in_min, in_max, out_min, out_max)
  335.  
  336.         /* NOTICE THAT ALL PARAMETERS ARE INT */
  337.         /* in_val = number to be scaled */
  338.         /* in_min = smallest possible in_val */
  339.         /* in_max = largest possible in_val */
  340.         /* out_min = smallest acceptable scaling value */
  341.         /* out_max = largest acceptable scaling value */
  342.         /* to scale x-values to the screen, out_min will */
  343.         /* usually be larger than out_max.  This has the effect */
  344.         /* of making the screen appear to be first-quadrant */
  345.         /* instead of absolute value of fourth-quadrant. */
  346.  
  347.  
  348. int in_val, in_min, in_max, out_min, out_max;
  349.  
  350. {
  351.     int out_val;                /* return value variable */
  352.     int largest, smallest;      /* out_min, out_max ordered by size */
  353.     double scale_factor;        /* ratio of out:in ranges */
  354.  
  355.  
  356.  
  357.     scale_factor = (double)(out_max - out_min) / (in_max - in_min);
  358.      /* without the DOUBLE cast, scale_factor values would be INT */
  359.      /* because all rvalues in the expression are INT */
  360.  
  361.     out_val = (int)((in_val - in_min)*scale_factor + out_min);
  362.      /* out_val will be INT without the cast, but why take chances? */
  363.  
  364.     if(out_min > out_max)   /* define largest and smallest return values */
  365.     {
  366.     smallest = out_max;
  367.     largest = out_min;
  368.     }
  369.     else
  370.     {
  371.     smallest = out_min;
  372.     largest = out_max;
  373.     }
  374.  
  375.         /* make sure integerizing has not put out_val out of bounds */
  376.     if(out_val > largest)
  377.         out_val = largest;
  378.  
  379.     if(out_val < smallest)
  380.         out_val = smallest;
  381.  
  382.  
  383.     return out_val;
  384.  
  385. }
  386.  
  387.  
  388.  
  389.  
  390.  
  391.  
  392.  
  393. /* long_scale() is a routine to scale histogram long integer values */
  394. /* to fit on a screen graph.  N.B.: Returns int, not long.          */
  395. /* moved from herc.c 10/89 by T Clune */
  396.  
  397. int long_scale(in_val, in_min, in_max, out_min, out_max)
  398.  
  399.         /* in_val = number to be scaled */
  400.         /* in_min = smallest possible in_val */
  401.         /* in_max = largest possible in_val */
  402.         /* out_min = smallest acceptable scaling value */
  403.         /* out_max = largest acceptable scaling value */
  404.         /* to scale y-values to the screen, out_min will */
  405.         /* usually be larger than out_max.  This has the effect */
  406.         /* of making the screen appear to be first-quadrant */
  407.         /* instead of absolute value of fourth-quadrant. */
  408.  
  409.  
  410. long int in_val, in_min, in_max;
  411. int out_min, out_max;
  412.  
  413. {
  414.     int out_val;                /* return value variable */
  415.     int largest, smallest;      /* out_min, out_max ordered by size */
  416.     double scale_factor;        /* ratio of out:in ranges */
  417.  
  418.  
  419.     if(in_val>=in_max)  /* kludgey patch for histarray[0] & [255] anomalies */
  420.     return(out_max);
  421.  
  422.     scale_factor = (double)(out_max - out_min) / (in_max - in_min);
  423.      /* without the DOUBLE cast, scale_factor values would be INT */
  424.      /* because all r-values in the expression are INT */
  425.  
  426.     out_val = (int)((in_val - in_min)*scale_factor + out_min);
  427.      /* out_val will be INT without the cast, but why take chances? */
  428.  
  429.     if(out_min > out_max)   /* define largest and smallest return values */
  430.     {
  431.     smallest = out_max;
  432.     largest = out_min;
  433.     }
  434.     else
  435.     {
  436.     smallest = out_min;
  437.     largest = out_max;
  438.     }
  439.  
  440.         /* make sure integerizing has not put out_val out of bounds */
  441.     if(out_val > largest)
  442.         out_val = largest;
  443.  
  444.     if(out_val < smallest)
  445.         out_val = smallest;
  446.  
  447.  
  448.     return out_val;
  449.  
  450. }
  451.  
  452.  
  453.  
  454.  
  455.  
  456. /* minmax() finds the min and max of the data set array.  n is the number of */
  457. /* points in the data set */
  458.  
  459. void minmax(array,n,minval,maxval)
  460. double array[];
  461. int n;
  462. double *minval, *maxval;
  463. {
  464.     int i;
  465.  
  466.     *minval= *maxval=array[0];
  467.     for(i=1;i<n;i++)
  468.     {
  469.     if(array[i]< *minval)
  470.         *minval=array[i];
  471.     if(array[i]> *maxval)
  472.         *maxval=array[i];
  473.     }
  474.  
  475. }
  476.  
  477.  
  478.  
  479.  
  480.  
  481.  
  482.  
  483. /* polar_to_xy() converts the amplitude,phase data into x,y data. */
  484. /* n is the number of data points. */
  485.  
  486. void polar_to_xy(amp, phase, x, y, n)
  487. double amp[], phase[], x[], y[];
  488. int n;
  489. {
  490.     int i;
  491.     double t_amp, t_phase;
  492.     for(i=0;i<n;i++)
  493.     {
  494.     t_amp=amp[i];   /* temp variables allow amp, phase to be x,y */
  495.     t_phase=phase[i];
  496.     x[i]=t_amp*cos(t_phase);
  497.     y[i]=t_amp*sin(t_phase);
  498.     }
  499.  
  500. }
  501.  
  502.  
  503.  
  504.  
  505. /* unscramble_transform() places the negative frequency components at the */
  506. /* beginning of the data array, so the x-axis of the graph is continuous */
  507. /* calling this function twice will restore the data array to its original */
  508. /* order.   If both amplitude & phase or x & y are used, call this function */
  509. /* twice, once with each array as an argument.  */
  510.  
  511. void unscramble_transform(array, n)
  512. double array[];
  513. int n;
  514. {
  515.     int i, j;
  516.     double temp;
  517.  
  518.     /* place negative freqs at beginning of fft array. */
  519.     j=n/2;
  520.     for(i=0;i<j;i++)
  521.     {
  522.     temp=array[i];
  523.     array[i]=array[i+j];
  524.     array[i+j]=temp;
  525.     }
  526.  
  527. }
  528.  
  529.  
  530.  
  531.  
  532. /* xy_to_polar() converts from rectangular coords to polar coords */
  533. /* n is the number of points in each array.  MINVAL is the threshold */
  534. /* value for an amplitude to be considered real.  It is used to decide */
  535. /* whether the PHASE value is real or just noise.  If AMP[i] is <= MINVAL */
  536. /* PHASE[i] is set to zero. xy_to_polar() allows x,y and amp,phase */
  537. /* to be the same arrays. phase is an angle from 0 to 2*PI. */
  538.  
  539. void xy_to_polar(x, y, amp, phase, n, minval)
  540. double x[], y[], amp[], phase[];
  541. int n;
  542. double minval;
  543. {
  544.     int i;
  545.     double tx, ty;
  546.  
  547.     for(i=0;i<n;i++)
  548.     {
  549.     tx=x[i];    /* the temp variables allow AMP,PHASE to be X,Y */
  550.     ty=y[i];
  551.     amp[i]=sqrt(tx*tx+ty*ty);
  552.     if(amp[i]<=minval)
  553.         phase[i]=0.0;
  554.     else
  555.         phase[i]=atan2(ty,tx);
  556.  
  557.         /* shift negative angles to 180 to 360 for compatibility with */
  558.         /* polar_to_xy() conversion format */
  559.  
  560.     if(phase[i]<0.0)
  561.         phase[i]=2*PI+phase[i];
  562.  
  563.     }
  564.  
  565. }
  566.