home *** CD-ROM | disk | FTP | other *** search
/ Gold Fish 2 / goldfish_vol2_cd1.bin / files / misc / sci / gfft / source / okspctrm.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-05-25  |  11.1 KB  |  407 lines
  1. /***************************************************************************
  2.  *          Copyright (C) 1994  Charles P. Peterson                  *
  3.  *         4007 Enchanted Sun, San Antonio, Texas 78244-1254             *
  4.  *              Email: Charles_P_Peterson@fcircus.sat.tx.us                *
  5.  *                                                                         *
  6.  *          This is free software with NO WARRANTY.                  *
  7.  *          See gfft.c, or run program itself, for details.              *
  8.  *              Support is available for a fee.                      *
  9.  ***************************************************************************
  10.  *
  11.  * Program:     gfft--General FFT analysis
  12.  * File:        okspctrm.c
  13.  * Purpose:     OK! Do an fft spectrum analysis
  14.  * Author:      Charles Peterson (CPP)
  15.  * History:     24-July-1993 CPP; Created.
  16.  * Comment:
  17.  *
  18.  *    As used here, a 'spectrum analysis' is one which starts with real
  19.  *    samples, and ends up with real magnitudes of some sort (possibly with
  20.  *    some sort of 'windowing' applied to the input which may be analyzed
  21.  *    in 'segments' which might be padded and/or overlapped, and possibly
  22.  *    with normalization applied to the output).  The output magnitudes may
  23.  *    either be a 'power' spectrum or an 'amplitude' or 'voltage' spectrum
  24.  *    (the latter being the positive square root of the former.)  *
  25.  *
  26.  *    (Ordinary fft, in contrast, may begin with complex input values and
  27.  *    usually returns complex output values.  Usually the input to an
  28.  *    ordinary fft is not segmented, and doing so may present additional
  29.  *    problems.  Most frequently, ordinary fft is actually part of some
  30.  *    larger problem, such as convolution or spectrum analysis.)
  31.  *
  32.  *    This function has gotten a little out of hand, primarily because of
  33.  *    all the possible combinations of options that are applied here.
  34.  *    Possibly 'Interleave' might be left out next time unless someone
  35.  *    notes a use for it.  (It didn't turn out to be as useful as I
  36.  *    expected.  When curves made with different amounts of interleave are
  37.  *    spliced together, they don't meet, even in 'PSDensity' mode.)  Likewise
  38.  *    for 'Pad,' an option which I've never seen do anything but make
  39.  *    matters much worse.  But somehow, some people seem to think that it's
  40.  *    useful, or at least harmless.  It isn't either, in my experience.
  41.  *    Leaving out those two options would reduce this size of this function
  42.  *    substantially.  
  43.  *
  44.  *    But, don't expect any such streamlining to make a signficant (or
  45.  *    even measureable) difference in performance.  My lprof tests indicate
  46.  *    that this entire beast takes less than 1% of the total time under
  47.  *    typical conditions.  Considering that there is real work going on
  48.  *    here--even if you don't use the more esoteric options--that is quite
  49.  *    reasonable.  (As expected, cfft takes around 70% of the time under
  50.  *    typical conditions, which isn't unreasonable either.)
  51.  *
  52.  *    Anyway, if you want to see my best programming and comment writing,
  53.  *    go to cfft.c, the function which actually does the fft, though it
  54.  *    has gotten harrier since I added an accelerated mode.
  55.  *
  56.  */
  57.  
  58. #define USE_MEMMOVE 
  59. /* Also uses MEMCPY, where applicable */
  60. /* Measured with lprof, only about 0.3% improvement overall. */
  61. /* But, even with 1 bin, didn't make matters worse either. */
  62.  
  63. #include <time.h>  /* time and difftime */
  64. #include <string.h> /* memmove and memcpy */
  65.  
  66. #include "gfft.h"
  67. #include "settings.h"
  68.  
  69. double LoopTimeElapsed = 0.0;
  70.  
  71. static BIN_TYPE *bins = NULL;
  72. static ULONG number_bins = 0;
  73. static int segment_count = 0;
  74.  
  75. void do_re_output (void)
  76. {
  77.     if (bins && number_bins && segment_count)
  78.     {
  79.     ok_write (bins, (long) number_bins, (long) segment_count);
  80.     }
  81.     else
  82.     {
  83.     error_message (CANT_RE_OUTPUT);
  84.     RAISE_ERROR (NOTHING_SPECIAL);    /* longjmp outa here */
  85.     }
  86. }
  87.  
  88.  
  89. ULONG ok_spectrum (BOOLEAN do_it_for_real)
  90. {
  91.     static float *indata = NULL;
  92.     static float *overlapdata = NULL;
  93.     static float *interleavedata = NULL;  /* e.g. Aseg Bseg Aseg Bseg ... */
  94.     ULONG actually_read = 0;
  95.     ULONG number_samples = 1;
  96.     ULONG number_interleave_samples = 1;
  97.     ULONG half_samples = 0;
  98.     int more_data = TRUE;
  99.     int error_in_loop = FALSE;
  100.     int overlap = Overlap;     /* May be changed if overlap impossible */
  101.     int last_partial_segment_used = FALSE;
  102.     ULONG loop_count = 0;
  103.     ULONG i;
  104.     ULONG j;
  105.     ULONG ai;
  106.     ULONG aj;
  107.     time_t time_beginning;
  108.     time_t time_ending;
  109.  
  110.     number_bins = 0;
  111.  
  112.     Sample_Sum_Of_Squares = 0.0;
  113.     Sample_Frame_Count = 0;
  114.  
  115.     if (NumberBins == MAX_BINS)
  116.     {
  117.     ULONG total_actually_read;
  118.     total_actually_read = ok_read (NULL, NULL);
  119.     if (!total_actually_read)
  120.     {
  121.         error_message (NO_DATA_PRESENT);
  122.         RAISE_ERROR (NOTHING_SPECIAL);  /* longjmp outa here! */
  123.     }
  124.     else
  125.     {
  126.         ULONG next_power;
  127.         ULONG total_in_leaf = total_actually_read / Interleave;
  128.         if (!total_in_leaf)
  129.         {
  130.         error_message (INSUFFICIENT_DATA);
  131.         RAISE_ERROR (NOTHING_SPECIAL);  /* longjmp outa here! */
  132.         }
  133.         next_power = get_pos_power_2 (total_in_leaf);
  134.         if (Pad || next_power == total_in_leaf)
  135.         {
  136.         number_bins = next_power / 2; /* pad if required*/
  137.         }
  138.         else
  139.         {
  140.         number_bins = next_power / 4; /* truncate if required */
  141.         }
  142.         number_samples = number_bins * 2 * Interleave;
  143.         bins_d_message (total_actually_read, number_bins);
  144.     }
  145.     }
  146.     else
  147.     {
  148.     number_bins = (ULONG) NumberBins;
  149.     if (number_bins > 0)
  150.     {
  151.         number_samples = number_bins * 2 * Interleave;
  152.     }
  153.     else
  154.     {
  155.         number_samples = 1;
  156.     }
  157.     }
  158. /*
  159.  * Now we have number_samples and number_bins
  160.  */
  161.     if (!do_it_for_real) return number_bins;
  162.  
  163.     segment_count = 0;
  164.  
  165.     time_beginning = time (NULL);
  166.  
  167.     half_samples = number_samples / 2;
  168.     if (number_samples < 2)
  169.     {
  170. /*
  171.  * Can't overlap segments of length 1
  172.  */
  173.     overlap = FALSE;
  174.     }
  175.     if (bins)
  176.     {
  177.     gfree (bins);
  178.     }
  179.     if (indata)
  180.     {
  181.     gfree (indata);
  182.     }
  183.     if (overlapdata)
  184.     {
  185.     gfree (overlapdata);
  186.     overlapdata = NULL;
  187.     }
  188.     if (interleavedata)
  189.     {
  190.     gfree (interleavedata);
  191.     interleavedata = NULL;
  192.     }
  193.       
  194. /*
  195.  * Note that there is one extra bin for 0 Hz; bins must be pre-zeroed
  196.  */
  197.     bins = gcalloc ((number_bins + 1), (long) sizeof(BIN_TYPE), 
  198.             NOTHING_SPECIAL);
  199.     indata = gmalloc (number_samples * sizeof(float), 
  200.               NOTHING_SPECIAL);
  201.     if (overlap)
  202.     {
  203.     overlapdata = gmalloc (number_samples * sizeof(float), 
  204.                NOTHING_SPECIAL);
  205.     }
  206.     if (Interleave > 1)
  207.     {
  208.     number_interleave_samples = number_samples / Interleave;
  209.     if (!number_interleave_samples)
  210.     {
  211.         error_message (INSUFFICIENT_DATA);
  212.         RAISE_ERROR (NOTHING_SPECIAL);  /* longjmp outa here! */
  213.     }
  214.     interleavedata = gmalloc (number_interleave_samples * sizeof (float),
  215.                  NOTHING_SPECIAL);
  216.     }
  217.  
  218. /*
  219.  * Now that these preliminaries are out of the way, we can
  220.  * begin looping over input data
  221.  */
  222.     time_beginning = time (NULL);
  223.  
  224.     while (more_data)
  225.     {
  226.     actually_read = ok_read (indata, number_samples);
  227.     if (!actually_read)
  228.     {
  229.         more_data = FALSE;
  230.         if (!loop_count)
  231.         {
  232.         error_message (NO_DATA_PRESENT);
  233.         error_in_loop = TRUE;
  234.         }
  235.         break;
  236.     }
  237.     if (actually_read < number_samples)  /* Deal with partial segment */
  238.     {
  239.         more_data = FALSE;
  240.         if (overlap && loop_count)
  241.         {
  242.         for (i = 0, j = actually_read; j < number_samples; 
  243.              i++, j++)
  244.         {
  245.             overlapdata[i] = overlapdata[j];
  246.         }
  247.         for (j = 0; j < actually_read; i++, j++)
  248.         {
  249.             overlapdata[i] = indata[j];
  250.         }
  251.         if (Interleave > 1)    /* Overlap and Interleave */
  252.         {
  253.             for (ai = 0; ai < Interleave; ai++)
  254.             {
  255.             float *a_data = interleavedata;
  256.             for (aj = ai; aj < number_samples; aj += Interleave)
  257.             {
  258.                 *a_data++ = overlapdata[aj];
  259.             }
  260.             ok_window (interleavedata, number_interleave_samples);
  261.             ok_rfft (interleavedata, number_interleave_samples);
  262.             ok_sigma (interleavedata, bins, number_bins);
  263.             segment_count++;
  264.             }
  265.         }
  266.         else  /* Overlap only...no Interleave */
  267.         {
  268.             ok_window (overlapdata, number_samples);
  269.             ok_rfft (overlapdata, number_samples);
  270.             ok_sigma (overlapdata, bins, number_bins);
  271.             segment_count++;
  272.         }
  273.         overlap = FALSE;  /* overlap buffer now used up */
  274.         last_partial_segment_used = TRUE;
  275.         }
  276.         if (Pad)
  277.         {
  278.         error_message (PADDING_TAILEND);
  279.         for (i = actually_read; i < number_samples; i++)
  280.         {
  281.             indata[i] = 0.0;
  282.         }
  283.         }
  284.         else
  285.         {
  286.         if (!segment_count)
  287.         {
  288.             error_message (INSUFFICIENT_DATA);
  289.             error_in_loop = TRUE;
  290.         }
  291. /*        else if (!last_partial_segment_used)
  292.         {
  293.             error_message (IGNORING_TAILEND);
  294.         } */
  295.         break;
  296.         }
  297.     }
  298. /*
  299.  * Here, we've got a full segment to work with (padded or not)
  300.  */
  301.     if (overlap)
  302.     {
  303.         if (loop_count)  /* Can't process overlap first time around */
  304.         {
  305. #ifdef USE_MEMMOVE
  306.         memmove (overlapdata, &overlapdata[half_samples], 
  307.              half_samples * sizeof (float));
  308.         memcpy (&overlapdata[half_samples], indata,
  309.             half_samples * sizeof (float));
  310. #else
  311.         for (i = 0, j = half_samples; i < half_samples; i++, j++)
  312.         {
  313.             overlapdata[i] = overlapdata[j];
  314.         }
  315.         for (j = 0; j < half_samples; i++, j++)
  316.         {
  317.             overlapdata[i] = indata[j];
  318.         }
  319. #endif
  320.         if (Interleave > 1)   /* Overlap and Interleave */
  321.         {
  322.             for (ai = 0; ai < Interleave; ai++)
  323.             {
  324.             float *a_data = interleavedata;
  325.             for (aj = ai; aj < number_samples; aj += Interleave)
  326.             {
  327.                 *a_data++ = overlapdata[aj];
  328.             }
  329.             ok_window (interleavedata, number_interleave_samples);
  330.             ok_rfft (interleavedata, number_interleave_samples);
  331.             ok_sigma (interleavedata, bins, number_bins);
  332.             segment_count++;
  333.             }
  334.         }
  335.         else  /* Overlap only...no Interleave */
  336.         {
  337.  
  338.             ok_window (overlapdata, number_samples);
  339.             ok_rfft (overlapdata, number_samples);
  340.             ok_sigma (overlapdata, bins, number_bins);
  341.             segment_count++;
  342.         }
  343.         }
  344. #ifdef USE_MEMMOVE
  345.         memcpy (overlapdata, indata, number_samples * sizeof (float));
  346. #else
  347.         for (i = 0; i < number_samples; i++)
  348.         {
  349.         overlapdata[i] = indata[i];
  350.         }
  351. #endif
  352.     }
  353.     if (Interleave > 1)  /* Regular (unoverlapped) alternation */
  354.     {
  355.         for (ai = 0; ai < Interleave; ai++)
  356.         {
  357.         float *a_data = interleavedata;
  358.         for (aj = ai; aj < number_samples; aj += Interleave)
  359.         {
  360.             *a_data++ = indata[aj];
  361.         }
  362.         ok_window (interleavedata, number_interleave_samples);
  363.         ok_rfft (interleavedata, number_interleave_samples);
  364.         ok_sigma (interleavedata, bins, number_bins);
  365.         segment_count++;
  366.         }
  367.     }
  368.     else
  369.     {
  370.         ok_window (indata, number_samples);
  371.         ok_rfft (indata, number_samples);
  372.         ok_sigma (indata, bins, number_bins);
  373.         segment_count++;
  374.     }
  375.     loop_count++;
  376.     }
  377.     time_ending = time (NULL);
  378.     LoopTimeElapsed = difftime (time_ending, time_beginning);
  379.  
  380.     if (!error_in_loop)
  381.     {
  382.     loop_time_message (LoopTimeElapsed);
  383.     ok_write (bins, number_bins, segment_count);
  384.     }
  385.     gfree (indata);
  386.     indata = NULL;
  387. /*    gfree (bins); save for ReOut */
  388. /*    bins = NULL; */
  389.     if (overlapdata)
  390.     {
  391.     gfree (overlapdata);
  392.     overlapdata = NULL;
  393.     }
  394.     if (interleavedata)
  395.     {
  396.     gfree (interleavedata);
  397.     interleavedata = NULL;
  398.     }
  399.     if (error_in_loop)
  400.     {
  401.     RAISE_ERROR (NOTHING_SPECIAL);  /* longjmp outa here! */
  402.     }
  403.     return number_bins;
  404. }
  405.  
  406.     
  407.