home *** CD-ROM | disk | FTP | other *** search
/ Gold Fish 2 / goldfish_vol2_cd1.bin / files / misc / sci / gfft / source / calibrat.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-05  |  6.1 KB  |  273 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:        calibrat.c
  13.  * Purpose:     calibration(s) 
  14.  * Author:      Charles Peterson (CPP)
  15.  * History:     15-March-1994 CPP; Created.
  16.  *              5-Aug-94 CPP (1.05); add new calibrations to end of list
  17.  * Comment:     Simple 'relative' calibration with linear interpolation.
  18.  *              Any number of calibration(s) may be in effect.
  19.  *              NOT complex deconvolution (maybe in a few years).
  20.  *              Linear interpolation, and flat extension are used
  21.  *              (Flat extension helps with fuzzy endpoints).
  22.  */
  23.  
  24. #include <stdio.h>
  25. #include <stdlib.h>  /* strtod */
  26. #include <string.h>  /* strtok */
  27. #include <math.h>    /* log10, pow */
  28.  
  29. #include "gfft.h"
  30. #include "settings.h"
  31.  
  32. #define CAL_BUFSIZ 133
  33.  
  34. int calibration_list__count (struct cal_st **clistp)
  35. {
  36.     int count = 0;
  37.     struct cal_st *clist = *clistp;
  38.  
  39.     while (clist)
  40.     {
  41.     count++;
  42.     clist = clist->next;
  43.     }
  44.     return count;
  45. }
  46.  
  47.  
  48. void calibration_list__cancel (struct cal_st **clistp)
  49. {
  50.     struct cal_st *clist = *clistp;
  51. /*
  52.  * Free tail of list recursively, then free this node
  53.  */
  54.     if (clist)
  55.     {
  56.     if (clist->next)
  57.     {
  58.         calibration_list__cancel (&clist->next);
  59.     }
  60.     gfree (clist->frequencies);
  61.     gfree (clist->amplitudes);
  62. /*    gfree (clist->filename); potentially dangerous */
  63.     gfree (clist);
  64.     }
  65.     *clistp = NULL;
  66. }
  67.  
  68.  
  69. void calibration_list__add (struct cal_st **clistp,
  70.                 char *cname, 
  71.                 FILE *cfile,
  72.                 BOOLEAN db_scale)
  73. {
  74.     struct cal_st *clist = *clistp;
  75.     char cal_buf[CAL_BUFSIZ];
  76.     char *token;
  77.     char *beyondpointer;
  78.     double *freq_array = NULL;
  79.     float *ampl_array = NULL;
  80.     double frequency;
  81.     float amplitude;
  82.     long size = 0;
  83.     struct cal_st *new_cal;
  84. /*
  85.  * First, read in file data
  86.  */
  87.     while (fgets (cal_buf, CAL_BUFSIZ, cfile))
  88.     {
  89.     /*
  90.      * Skip over comment lines
  91.      */
  92.     if (cal_buf[0] == '#' || cal_buf[0] == ';' || cal_buf[0] == '!')
  93.     {
  94.         continue;
  95.     }
  96.     /*
  97.      * get frequency value
  98.      */
  99.     token = strtok (cal_buf, " ");
  100.     if (!token)
  101.     {
  102.         error_message (INVALID_CALIBRATION_FILE);
  103.         RAISE_ERROR (NOTHING_SPECIAL);    /* longjmp outa here */
  104.     }
  105.     frequency = strtod (token, &beyondpointer);
  106.     if (*beyondpointer != '\0')
  107.     {
  108.         error_message (INVALID_CALIBRATION_FILE);
  109.         RAISE_ERROR (NOTHING_SPECIAL);    /* longjmp outa here */
  110.     }
  111.     token = strtok (NULL, " ");
  112.     if (!token)
  113.     {
  114.         error_message (INVALID_CALIBRATION_FILE);
  115.         RAISE_ERROR (NOTHING_SPECIAL);    /* longjmp outa here */
  116.     }
  117.     amplitude = (float) strtod (token, &beyondpointer);
  118.     if (*beyondpointer != '\0' && *beyondpointer != '\n')
  119.     {
  120.         error_message (INVALID_CALIBRATION_FILE);
  121.         RAISE_ERROR (NOTHING_SPECIAL);    /* longjmp outa here */
  122.     }
  123.     /*
  124.      * Add new values to arrays
  125.      */
  126.     size++;
  127.     freq_array = grealloc (freq_array, (size * sizeof frequency), 
  128.                    NOTHING_SPECIAL);
  129.     ampl_array = grealloc (ampl_array, (size * sizeof amplitude), 
  130.                    NOTHING_SPECIAL);
  131.     freq_array[size-1] = frequency;
  132.     ampl_array[size-1] = amplitude;
  133.     }
  134.     if (size <= 0)
  135.     {
  136.     error_message (INVALID_CALIBRATION_FILE);
  137.     RAISE_ERROR (NOTHING_SPECIAL);    /* longjmp outa here */
  138.     }
  139.  
  140. /*
  141.  * create new node
  142.  */
  143.     new_cal = gmalloc (sizeof (struct cal_st), NOTHING_SPECIAL);
  144.     new_cal->next = NULL;
  145.     new_cal->filename = cname;
  146.     new_cal->size = size;
  147.     new_cal->db = db_scale;
  148.     new_cal->frequencies = freq_array;
  149.     new_cal->amplitudes = ampl_array;
  150.     new_cal->index = 0;
  151. /*
  152.  * insert new node in list, or at start
  153.  */
  154.     if (clist)
  155.     {
  156.     while (clist->next)
  157.     {
  158.         clist = clist->next;
  159.     }
  160.     clist->next = new_cal;
  161.     }
  162.     else
  163.     {
  164.     *clistp = new_cal;
  165.     }
  166. }
  167.  
  168. void calibration_list__reset (struct cal_st **clistp)
  169. {
  170.     if (*clistp)
  171.     {
  172.     (*clistp)->index = 0;
  173.     calibration_list__reset (&((*clistp)->next));
  174.     }
  175. }
  176.  
  177.  
  178. float calibration_list__apply (struct cal_st **clistp,
  179.                    float value, double frequency)
  180. {
  181.     struct cal_st *clist = *clistp;
  182.     long index;
  183.     BOOLEAN outrange = TRUE;
  184.     double factor;
  185.  
  186. /*
  187.  * Assumption is made that frequencies are applied in an increasing
  188.  * order.  When starting over, call calibration_list__reset
  189.  *
  190.  * Make local copy of 'index'
  191.  */
  192.     index = clist->index;
  193. /*
  194.  * Find applicable index
  195.  */
  196.     while (frequency > clist->frequencies[index])
  197.     {
  198.     outrange = FALSE;
  199.     if (++index >= clist->size)
  200.     {
  201.         index--;
  202.         outrange = TRUE;
  203.         break;
  204.     }
  205.     }
  206.     clist->index = index;  /* update */
  207.     
  208. /*
  209.  * Compute the factor to use
  210.  */
  211.  
  212.     if (frequency == clist->frequencies[index] || outrange)
  213.     {
  214.     factor = clist->amplitudes[index];
  215.     }
  216.     else
  217.     {
  218.     /*
  219.      * linear interpolation is fun
  220.      * (Sorry I don't bother with anything more sophisticated yet)
  221.      */
  222.     double x1 = clist->frequencies[index-1];
  223.     double x2 = clist->frequencies[index];
  224.     double y1 = clist->amplitudes[index-1];
  225.     double y2 = clist->amplitudes[index];
  226.     double slope = (y2 - y1) / (x2 - x1);
  227.     double offset = y1 - (slope * x1);
  228.  
  229.     factor = (slope * frequency) + offset;
  230.     }
  231.  
  232. /*
  233.  * Apply
  234.  */
  235.     if (Db)
  236.     {
  237.     if (clist->db)
  238.     {
  239.         value -= factor;
  240.     }
  241.     else
  242.     {
  243.         double dbfactor;
  244.         if (factor <= 0.0) RAISE_ERROR (NOTHING_SPECIAL);
  245.         dbfactor = 20.0 * log10 (factor);
  246.         value -= dbfactor;
  247.     }
  248.     }
  249.     else
  250.     {
  251.     if (!clist->db)
  252.     {
  253.         if (factor == 0.0) RAISE_ERROR (NOTHING_SPECIAL);
  254.         value /= factor;
  255.     }
  256.     else
  257.     {
  258.         value /= pow (10.0, (factor * 0.05));
  259.     }
  260.     }
  261. /*
  262.  * Do next calibration in list
  263.  */
  264.     if (clist->next)
  265.     {
  266.     value = calibration_list__apply (&clist->next, value, frequency);
  267.     }
  268.  
  269.     return value;
  270. }
  271.     
  272.  
  273.