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