home *** CD-ROM | disk | FTP | other *** search
/ Gold Fish 2 / goldfish_vol2_cd1.bin / files / misc / sci / gfft / source / okwrite.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-29  |  10.3 KB  |  413 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:        okwrite.c
  13.  * Purpose:     Write the results of a spectrum analysis
  14.  * Author:      Charles Peterson (CPP)
  15.  * History:     23-August-1993 CPP; Created.
  16.  *              6-Aug-1994 CPP (1.10); Add spectrum file identifier #FFT
  17.  *              28-Aug-94 CPP (1.12); fix LowFrequency for FFP
  18.  * Comment:     This sure has gotten complicated.  Sorry.
  19.  */
  20.  
  21. /* #define PARSEVAL_DEBUG */
  22.  
  23. #include <stdio.h>
  24. #include <math.h>
  25.  
  26. #ifdef _AMIGA
  27. #ifdef _M68881
  28. #include <m68881.h>
  29. #endif
  30. #endif
  31.  
  32. #include "gfft.h"
  33. #include "settings.h" /* OkRate, WritePtr, Mean, Amplitude */
  34. #include "errcodes.h"
  35.  
  36. double LowestFrequencyWritten = 0;
  37. double HighestFrequencyWritten = 0;
  38.  
  39. void ok_write (BIN_TYPE *bins, long number_bins, long number_segments)
  40. {
  41.     long i;
  42.     double frequency;
  43.     FINAL_TYPE value;
  44. /*    double two_doubles[2]; */
  45.     double nyquist_frequency = OkRate / (2.0L * Interleave);
  46.     double delta_frequency = 0;
  47.     double value_divisor_temp = 1.0 * number_segments;
  48.     FINAL_TYPE multiply = Multiply;
  49.     FINAL_TYPE value_divisor;
  50.     double *mesh;
  51.     FINAL_TYPE smoothing_value_sum = 0.0L;
  52.     double smoothing_frequency_sum = 0.0L;
  53.     FINAL_TYPE save_value;
  54.     double save_frequency;
  55.     double test_frequency;
  56.     long mesh_index = 0;
  57.     long smoothing_count = 0;
  58.     BOOLEAN first_db_of_zero = TRUE;
  59.     double sum_of_bins;
  60.  
  61.  
  62.     fprintf (WritePtr, "#FFT\n");
  63.  
  64.     if (number_segments == 0) return;
  65.  
  66. /*
  67.  * reset index used internally in calibration
  68.  * because we're starting from lowest frequency (again, maybe)
  69.  */
  70.     calibration_list__reset (&CalibrationList);
  71.  
  72. /*
  73.  * If 0 bins (i.e., the DC one only)
  74.  * avoid divide by zero problems
  75.  */
  76.     if (number_bins)
  77.     {
  78.     delta_frequency = nyquist_frequency / number_bins;
  79.     mesh  = ok_mesh (nyquist_frequency, delta_frequency);
  80.  
  81.     if (Mean)
  82.     {
  83.     /*
  84.          * Mean normalization for power or amplitude.
  85.      * See Press, et. al, page 439.
  86.      * To get mean, we'll have to divide by N^2, then by number_segments
  87.          * Note: N is number_bins * 2 if number_bins > 0, thus
  88.      *  (number_bins * 2)^2 == 4 * number_bins^2
  89.          * 
  90.          * READ THIS: if Parseval option selected (now default)
  91.      *  we recompute this later using Parseval's theorem.  The
  92.      *  results are remarkably similar...showing that I WAS doing
  93.      *  the correct thing here.  The new Parseval method is
  94.      *  more accurate, but the old method is 2.5% faster overall
  95.      *  because the squaring and summation in ok_read is not needed.
  96.          */
  97.         value_divisor_temp = (((4.0 * number_bins) * number_bins) * 
  98.                 number_segments);
  99.     }
  100.     else /* !Mean */
  101.     {
  102.  
  103. /*          value_divisor_temp = 2.0 * number_bins; */
  104.  
  105. /*
  106.  * We correct for overlapping bins, by multiplying by extra factor
  107.  * (number_segments * <segment size> / Sample_Frame_Count)
  108.  */
  109.  
  110. /*        value_divisor_temp = 2.0 * number_bins * number_segments *
  111.         2.0 * number_bins / Sample_Frame_Count */
  112.         
  113. /*
  114.  * which simplifies to the following:
  115.  */
  116.  
  117.         value_divisor_temp = 4.0 * number_bins * number_bins *
  118.           number_segments / Sample_Frame_Count;
  119.     }
  120.     }
  121.  
  122.     if (WindowType)
  123. /*
  124.  * Special extra correction if non-square window is used
  125.  * NOTE: as currently done, this is an estimate
  126.  * based on window's mean squared average gain...
  127.  * Parseval method works better.
  128.  */
  129.     {
  130.     value_divisor_temp *= ok_window__gain2();
  131.     }
  132.  
  133. /*
  134.  * This is now where the normalization normally gets done
  135.  * Using Parseval's theorem, i.e.
  136.  *   sum of power in the time domain equals sum of power in freq domain
  137.  */
  138.     if (Parseval)
  139.     {
  140.     sum_of_bins = bins[0];
  141.     if (number_bins > 0)
  142.     {
  143.         sum_of_bins += bins[number_bins];
  144.     }
  145.     for (i = 1; i <= number_bins-1; i++)
  146.     {
  147.         sum_of_bins += bins[i] * 2.0;
  148.     }
  149.  
  150.     if (Mean && Sample_Sum_Of_Squares != 0.0)
  151.     {
  152.  
  153. #ifdef PARSEVAL_DEBUG
  154.     fprintf (stderr,"Sample MSS: %g, Bin MSS: %g, Old Divisor: %g\n",
  155.          Sample_Sum_Of_Squares / Sample_Frame_Count,
  156.          sum_of_bins / value_divisor_temp,
  157.          value_divisor_temp);
  158. #endif
  159.  
  160.         value_divisor_temp = sum_of_bins * Sample_Frame_Count /
  161.           Sample_Sum_Of_Squares;
  162.  
  163. #ifdef PARSEVAL_DEBUG
  164.     fprintf (stderr,"Corrected Sample Bin MSS: %g\n",
  165.          sum_of_bins / value_divisor_temp);
  166. #endif
  167.  
  168.     }
  169.     else if (!Mean && Sample_Sum_Of_Squares != 0.0)
  170.     {
  171. #ifdef PARSEVAL_DEBUG
  172.         fprintf (stderr,"Original divisor: %g\n",value_divisor_temp);
  173. #endif
  174.  
  175.         value_divisor_temp = sum_of_bins / Sample_Sum_Of_Squares;
  176.  
  177. #ifdef PARSEVAL_DEBUG
  178.        fprintf (stderr,
  179.         "Post correction: Sample SS: %g, Bin SS: %g, Divisor: %g\n",
  180.         Sample_Sum_Of_Squares,
  181.         sum_of_bins / value_divisor_temp,
  182.         value_divisor_temp);
  183. #endif
  184.  
  185.     }
  186.     /* else...all's zeros anyway */
  187.     }
  188.  
  189.     if (PSDensity)
  190.     {
  191.     /*
  192.      * (Re-)Adjust for PSDensity:
  193.      * The effect of this is to allow you to splice together curves
  194.      * computed with different numbers of bins, which is useful in some
  195.      * applications.
  196.      *
  197.      * Ok, we could multiply back in an extra N,
  198.      * but, I've chosen to divide by the bin size
  199.      * giving more meaningful units of (.../Hz) or something like that,
  200.      * I think...
  201.      */
  202.     value_divisor_temp *= delta_frequency;
  203.     }
  204.  
  205.  
  206.     value_divisor = value_divisor_temp;  /* reduce to final precision */
  207.  
  208.     LowestFrequencyWritten = -1.0;
  209. /*
  210.  * Thanks to bug in FFP handing of DBL_MIN, LOWEST_FREQUENCY is now -1
  211.  * and I have to do this extra stuff
  212.  */
  213.     if (LowFrequency == LOWEST_FREQUENCY)
  214.     {
  215.     i = 1;
  216.     }
  217.     else
  218.     {
  219.     i = 0;
  220.     }
  221.     for (; i <= number_bins; i++)
  222.     {
  223.     frequency = delta_frequency * i;  /* Note 888 below if changing */
  224.  
  225.     if (frequency < LowFrequency)
  226.     {
  227.         continue;
  228.     }
  229.     if (frequency > HighFrequency)
  230.     {
  231.         continue;
  232.     }
  233.  
  234.     value = bins[i];
  235.  
  236.     if (i != 0 && i != number_bins)
  237.     {
  238.         value *= 2.0;    /* Correct for one-sidedness */
  239.     }
  240.  
  241.     value /= value_divisor;
  242.  
  243.     if (Pink)
  244.     {
  245.         if (frequency == 0)
  246.         {
  247.         continue;
  248.         }
  249.         /*
  250.          * The following normalizes the spectrum values to about what
  251.          * would be reported by an octave band spectrum analyzer.
  252.          *  ** Warning ** This is a guess, not known to be correct!
  253.          */
  254.         value *= ((frequency + delta_frequency) / delta_frequency);
  255.     }
  256.  
  257.     if (Amplitude)
  258.     {
  259.         value = sqrt (value);
  260.     }
  261.  
  262.     if (multiply != 1.0)
  263.     {
  264.         value *= multiply;
  265.     }
  266.  
  267.     if (mesh && frequency > 0.0L)
  268.     {
  269.         if (frequency > mesh[mesh_index] && smoothing_count == 0)
  270.         {
  271.         while (frequency > mesh[++mesh_index]);  /* Catch up */
  272.         }
  273.         if (frequency < mesh[mesh_index])
  274.         {
  275.         smoothing_value_sum += (SquaredSmoothing) ?
  276.           value * value : value;
  277.         smoothing_frequency_sum += frequency;
  278.         smoothing_count++;
  279.         continue;  /* Nothing to write this time */
  280.         }
  281.         else if (frequency == mesh[mesh_index])
  282.         {
  283.         /* Set up value and frequency for this write */
  284.         smoothing_value_sum += (SquaredSmoothing) ?
  285.           value * value : value;
  286.         smoothing_frequency_sum += frequency;
  287.         value = smoothing_value_sum / ++smoothing_count;
  288.         if (SquaredSmoothing) value = sqrt (value);
  289.         frequency = smoothing_frequency_sum / smoothing_count;
  290.  
  291.         /* Set up for next pass */
  292.         smoothing_value_sum = 0.0;
  293.         smoothing_frequency_sum = 0.0;
  294.         smoothing_count = 0;
  295.         mesh_index++;
  296.         }
  297.         else if (frequency > mesh[mesh_index])
  298.         {
  299.         /* Set up value and frequency for this write */
  300.         save_value = value;
  301.         save_frequency = frequency;
  302.         value = smoothing_value_sum / smoothing_count;
  303.         if (SquaredSmoothing) value = sqrt (value);
  304.         frequency = smoothing_frequency_sum / smoothing_count;
  305.  
  306.         /* Set up for next pass */
  307.         smoothing_value_sum = (SquaredSmoothing) ?
  308.           save_value * save_value : save_value;
  309.         smoothing_frequency_sum = save_frequency;
  310.         smoothing_count = 1;
  311.         test_frequency = (i+1) * delta_frequency; /*Note 888 above*/
  312.         if (test_frequency > mesh[++mesh_index])
  313.         {   /* Catch up for next pass */
  314.             while (test_frequency > mesh[++mesh_index]);
  315.         }
  316.         }
  317.     }        
  318.     if (Db)
  319.     {
  320.         if (value > 0.0)
  321.         {
  322.         value = ((Amplitude) ? 20 : 10) * log10 (value);
  323.         }
  324.         else
  325.         {
  326.         if (first_db_of_zero)
  327.         {
  328.             error_message (DB_OF_ZERO);
  329.             first_db_of_zero = FALSE;
  330.         }
  331.         continue;  /* Just skip (unless user decides not to) */
  332.         }
  333.     }
  334.  
  335.     if (CalibrationList)
  336.     {
  337.         CATCH_ERROR
  338.         {
  339.         value = calibration_list__apply (&CalibrationList, 
  340.                          value, frequency);
  341.         }
  342.         ON_ERROR
  343.         {
  344.         if (first_db_of_zero)
  345.         {
  346.             error_message (DB_OF_ZERO);
  347.             first_db_of_zero = FALSE;
  348.         }
  349.         continue;  /* Just skip (unless user decides not to) */
  350.         }
  351.         END_CATCH_ERROR;
  352.     } /* end if CalibrationList */
  353.  
  354.     if (Quantization != MIN_QUANTIZATION)
  355.     {
  356.         long factor = value / Quantization;
  357.         double value1 = (factor-1) * Quantization;
  358.         double value2 = factor * Quantization;
  359.         double value3 = (factor+1) * Quantization;
  360.         double delta1 = value - value1;
  361.         double delta2 = fabs (value - value2);
  362.         double delta3 = value3 - value;
  363.         if (delta2 <= delta1 && delta2 <= delta3)
  364.         {
  365.         value = value2;
  366.         }
  367.         else if (delta1 <= delta3)
  368.         {
  369.         value = value1;
  370.         }
  371.         else
  372.         {
  373.         value = value3;
  374.         }
  375.     }
  376.  
  377.     if (LowestFrequencyWritten < 0) LowestFrequencyWritten = frequency;
  378.     if (!OutputFormat.binary)
  379.     {
  380.         if (Time3D)
  381.         {
  382. #ifdef _FFP
  383.         fprintf (WritePtr, "%-15.8g %-15.8g %-13.7g\n",
  384.              frequency, OkSegmentTime, value);
  385. #else
  386.         fprintf (WritePtr, "%-19.12g %-19.12g %-13.7g\n", 
  387.              frequency, OkSegmentTime, value);
  388. #endif
  389.         }
  390.         else
  391.         {
  392. #ifdef _FFP
  393.         fprintf (WritePtr, "%-15.8g %-13.7g\n", frequency, value);
  394. #else
  395.         fprintf (WritePtr, "%-19.12g %-13.7g\n", frequency, value);
  396. #endif
  397.         }
  398.     }
  399. /*    else  ** This didn't work...need to figure binary formats out **
  400.  *    {
  401.  *        two_doubles[0] = frequency;
  402.  *        two_doubles[1] = value;
  403.  *        fwrite (two_doubles, sizeof(double), (size_t) 2, WritePtr);
  404.  *    }
  405.  */        
  406.     }
  407.     if (mesh) gfree (mesh);
  408.     if (Time3D) fprintf (WritePtr, "\n");
  409.     HighestFrequencyWritten = frequency;
  410. }
  411.  
  412.       
  413.