home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / fft / real_init.c < prev    next >
Encoding:
Text File  |  2001-08-22  |  4.3 KB  |  184 lines
  1. /* fft/real_init.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. TYPE(gsl_fft_real_wavetable) *
  21. FUNCTION(gsl_fft_real_wavetable,alloc) (size_t n)
  22. {
  23.   int status;
  24.   size_t i;
  25.   size_t n_factors;
  26.   size_t t, product, product_1, q;
  27.   double d_theta;
  28.  
  29.   TYPE(gsl_fft_real_wavetable) * wavetable;
  30.  
  31.   if (n == 0)
  32.     {
  33.       GSL_ERROR_VAL ("length n must be positive integer", GSL_EDOM, 0);
  34.     }
  35.  
  36.   wavetable = (TYPE(gsl_fft_real_wavetable) *) 
  37.     malloc(sizeof(TYPE(gsl_fft_real_wavetable)));
  38.  
  39.   if (wavetable == NULL)
  40.     {
  41.       GSL_ERROR_VAL ("failed to allocate struct", GSL_ENOMEM, 0);
  42.     }
  43.  
  44.   if (n == 1) 
  45.     {
  46.       wavetable->trig = 0;
  47.     }
  48.   else
  49.     {
  50.       wavetable->trig = (TYPE(gsl_complex) *) 
  51.         malloc ((n / 2) * sizeof (TYPE(gsl_complex)));
  52.       
  53.       if (wavetable->trig == NULL)
  54.     {
  55.           /* error in constructor, prevent memory leak */
  56.           
  57.           free(wavetable) ; 
  58.  
  59.       GSL_ERROR_VAL ("failed to allocate trigonometric lookup table", 
  60.                 GSL_ENOMEM, 0);
  61.     }
  62.     }
  63.  
  64.   wavetable->n = n;
  65.  
  66.   status = fft_real_factorize (n, &n_factors, wavetable->factor);
  67.  
  68.   if (status)
  69.     {
  70.       /* error in constructor, prevent memory leak */
  71.       
  72.       free(wavetable->trig);
  73.       free(wavetable) ; 
  74.  
  75.       GSL_ERROR_VAL ("factorization failed", GSL_EFACTOR, 0);
  76.     }
  77.  
  78.   wavetable->nf = n_factors;
  79.  
  80.   d_theta = 2.0 * M_PI / ((double) n);
  81.  
  82.   t = 0;
  83.   product = 1;
  84.   for (i = 0; i < wavetable->nf; i++)
  85.     {
  86.       size_t j;
  87.       const size_t factor = wavetable->factor[i];
  88.       wavetable->twiddle[i] = wavetable->trig + t;
  89.       product_1 = product;    /* product_1 = p_(i-1) */
  90.       product *= factor;
  91.       q = n / product;
  92.  
  93.       for (j = 1; j < factor; j++)
  94.     {
  95.       size_t k;
  96.       size_t m = 0;
  97.       for (k = 1; k < (product_1 + 1) / 2; k++)
  98.         {
  99.           double theta;
  100.           m = m + j * q;
  101.           m = m % n;
  102.           theta = d_theta * m;    /*  d_theta*j*k*q */
  103.           GSL_REAL(wavetable->trig[t]) = cos (theta);
  104.           GSL_IMAG(wavetable->trig[t]) = sin (theta);
  105.  
  106.           t++;
  107.         }
  108.     }
  109.     }
  110.  
  111.   if (t > (n / 2))
  112.     {
  113.       /* error in constructor, prevent memory leak */
  114.  
  115.       free(wavetable->trig);
  116.       free(wavetable) ; 
  117.  
  118.       GSL_ERROR_VAL ("overflowed trigonometric lookup table", 
  119.                         GSL_ESANITY, 0);
  120.     }
  121.  
  122.   return wavetable;
  123. }
  124.  
  125. TYPE(gsl_fft_real_workspace) *
  126. FUNCTION(gsl_fft_real_workspace,alloc) (size_t n)
  127. {
  128.   TYPE(gsl_fft_real_workspace) * workspace;
  129.  
  130.   if (n == 0)
  131.     {
  132.       GSL_ERROR_VAL ("length n must be positive integer", GSL_EDOM, 0);
  133.     }
  134.  
  135.   workspace = (TYPE(gsl_fft_real_workspace) *) 
  136.     malloc(sizeof(TYPE(gsl_fft_real_workspace)));
  137.  
  138.   if (workspace == NULL)
  139.     {
  140.       GSL_ERROR_VAL ("failed to allocate struct", GSL_ENOMEM, 0);
  141.     }
  142.  
  143.   workspace->n = n;
  144.  
  145.   workspace->scratch = (BASE *) malloc (n * sizeof (BASE));
  146.  
  147.   if (workspace->scratch == NULL)
  148.     {
  149.       /* error in constructor, prevent memory leak */
  150.       
  151.       free(workspace) ; 
  152.  
  153.       GSL_ERROR_VAL ("failed to allocate scratch space", GSL_ENOMEM, 0);
  154.     }
  155.  
  156.   return workspace;
  157. }
  158.  
  159.  
  160.  
  161. void
  162. FUNCTION(gsl_fft_real_wavetable,free) (TYPE(gsl_fft_real_wavetable) * wavetable)
  163. {
  164.  
  165.   /* release trigonometric lookup tables */
  166.  
  167.   free (wavetable->trig);
  168.   wavetable->trig = NULL;
  169.  
  170.   free (wavetable) ;
  171. }
  172.  
  173. void
  174. FUNCTION(gsl_fft_real_workspace,free) (TYPE(gsl_fft_real_workspace) * workspace)
  175. {
  176.  
  177.   /* release scratch space */
  178.  
  179.   free (workspace->scratch);
  180.   workspace->scratch = NULL;
  181.  
  182.   free (workspace) ;
  183. }
  184.