single-fft.c

     
   1  //! @file single-fft.c
   2  //! @author J. Marcel van der Veer
   3  //!
   4  //! @section Copyright
   5  //!
   6  //! This file is part of Algol68G - an Algol 68 compiler-interpreter.
   7  //! Copyright 2001-2023 J. Marcel van der Veer [algol68g@xs4all.nl].
   8  //!
   9  //! @section License
  10  //!
  11  //! This program is free software; you can redistribute it and/or modify it 
  12  //! under the terms of the GNU General Public License as published by the 
  13  //! Free Software Foundation; either version 3 of the License, or 
  14  //! (at your option) any later version.
  15  //!
  16  //! This program is distributed in the hope that it will be useful, but 
  17  //! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 
  18  //! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for 
  19  //! more details. You should have received a copy of the GNU General Public 
  20  //! License along with this program. If not, see [http://www.gnu.org/licenses/].
  21  
  22  //! @section Synopsis
  23  //!
  24  //! REAL, COMPLEX fast fourier transform.
  25  
  26  #include "a68g.h"
  27  #include "a68g-genie.h"
  28  #include "a68g-prelude.h"
  29  
  30  #if defined (HAVE_GSL)
  31  
  32  //! @brief Map GSL error handler onto a68g error handler.
  33  
  34  void fft_error_handler (const char *reason, const char *file, int line, int gsl_errno)
  35  {
  36    if (line != 0) {
  37      ASSERT (snprintf (A68 (edit_line), SNPRINTF_SIZE, "%s in line %d of file %s", reason, line, file) >= 0);
  38    } else {
  39      ASSERT (snprintf (A68 (edit_line), SNPRINTF_SIZE, "%s", reason) >= 0);
  40    }
  41    diagnostic (A68_RUNTIME_ERROR, A68 (f_entry), ERROR_FFT, A68 (edit_line), gsl_strerror (gsl_errno));
  42    exit_genie (A68 (f_entry), A68_RUNTIME_ERROR);
  43  }
  44  
  45  //! @brief Detect math errors.
  46  
  47  void fft_test_error (int rc)
  48  {
  49    if (rc != 0) {
  50      fft_error_handler ("math error", "", 0, rc);
  51    }
  52  }
  53  
  54  //! @brief Pop [] REAL on the stack as complex REAL_T [].
  55  
  56  REAL_T *pop_array_real (NODE_T * p, int *len)
  57  {
  58    A68_REF desc;
  59    A68_ARRAY *arr;
  60    A68_TUPLE *tup;
  61    int inc, iindex, k;
  62    BYTE_T *base;
  63    REAL_T *v;
  64    A68 (f_entry) = p;
  65  // Pop arguments.
  66    POP_REF (p, &desc);
  67    CHECK_REF (p, desc, M_ROW_REAL);
  68    GET_DESCRIPTOR (arr, tup, &desc);
  69    *len = ROW_SIZE (tup);
  70    if ((*len) <= 0) {
  71      return NO_REAL;
  72    }
  73    v = (REAL_T *) get_heap_space (2 * (size_t) (*len) * sizeof (REAL_T));
  74    fft_test_error (v == NO_REAL ? GSL_ENOMEM : GSL_SUCCESS);
  75    base = DEREF (BYTE_T, &ARRAY (arr));
  76    iindex = VECTOR_OFFSET (arr, tup);
  77    inc = SPAN (tup) * ELEM_SIZE (arr);
  78    for (k = 0; k < (*len); k++, iindex += inc) {
  79      A68_REAL *x = (A68_REAL *) (base + iindex);
  80      CHECK_INIT (p, INITIALISED (x), M_REAL);
  81      v[2 * k] = VALUE (x);
  82      v[2 * k + 1] = 0.0;
  83    }
  84    return v;
  85  }
  86  
  87  //! @brief Push REAL_T [] on the stack as [] REAL.
  88  
  89  void push_array_real (NODE_T * p, REAL_T * v, int len)
  90  {
  91    A68_REF desc, row;
  92    A68_ARRAY arr;
  93    A68_TUPLE tup;
  94    int inc, iindex, k;
  95    BYTE_T *base;
  96    A68 (f_entry) = p;
  97    NEW_ROW_1D (desc, row, arr, tup, M_ROW_REAL, M_REAL, len);
  98    base = DEREF (BYTE_T, &ARRAY (&arr));
  99    iindex = VECTOR_OFFSET (&arr, &tup);
 100    inc = SPAN (&tup) * ELEM_SIZE (&arr);
 101    for (k = 0; k < len; k++, iindex += inc) {
 102      A68_REAL *x = (A68_REAL *) (base + iindex);
 103      STATUS (x) = INIT_MASK;
 104      VALUE (x) = v[2 * k];
 105      CHECK_REAL (p, VALUE (x));
 106    }
 107    PUSH_REF (p, desc);
 108  }
 109  
 110  //! @brief Pop [] COMPLEX on the stack as REAL_T [].
 111  
 112  REAL_T *pop_array_complex (NODE_T * p, int *len)
 113  {
 114    A68_REF desc;
 115    A68_ARRAY *arr;
 116    A68_TUPLE *tup;
 117    int inc, iindex, k;
 118    BYTE_T *base;
 119    REAL_T *v;
 120    A68 (f_entry) = p;
 121  // Pop arguments.
 122    POP_REF (p, &desc);
 123    CHECK_REF (p, desc, M_ROW_COMPLEX);
 124    GET_DESCRIPTOR (arr, tup, &desc);
 125    *len = ROW_SIZE (tup);
 126    if ((*len) <= 0) {
 127      return NO_REAL;
 128    }
 129    v = (REAL_T *) get_heap_space (2 * (size_t) (*len) * sizeof (REAL_T));
 130    fft_test_error (v == NO_REAL ? GSL_ENOMEM : GSL_SUCCESS);
 131    base = DEREF (BYTE_T, &ARRAY (arr));
 132    iindex = VECTOR_OFFSET (arr, tup);
 133    inc = SPAN (tup) * ELEM_SIZE (arr);
 134    for (k = 0; k < (*len); k++, iindex += inc) {
 135      A68_REAL *re = (A68_REAL *) (base + iindex);
 136      A68_REAL *im = (A68_REAL *) (base + iindex + SIZE (M_REAL));
 137      CHECK_INIT (p, INITIALISED (re), M_COMPLEX);
 138      CHECK_INIT (p, INITIALISED (im), M_COMPLEX);
 139      v[2 * k] = VALUE (re);
 140      v[2 * k + 1] = VALUE (im);
 141    }
 142    return v;
 143  }
 144  
 145  //! @brief Push REAL_T [] on the stack as [] COMPLEX.
 146  
 147  void push_array_complex (NODE_T * p, REAL_T * v, int len)
 148  {
 149    A68_REF desc, row;
 150    A68_ARRAY arr;
 151    A68_TUPLE tup;
 152    int inc, iindex, k;
 153    BYTE_T *base;
 154    A68 (f_entry) = p;
 155    NEW_ROW_1D (desc, row, arr, tup, M_ROW_COMPLEX, M_COMPLEX, len);
 156    base = DEREF (BYTE_T, &ARRAY (&arr));
 157    iindex = VECTOR_OFFSET (&arr, &tup);
 158    inc = SPAN (&tup) * ELEM_SIZE (&arr);
 159    for (k = 0; k < len; k++, iindex += inc) {
 160      A68_REAL *re = (A68_REAL *) (base + iindex);
 161      A68_REAL *im = (A68_REAL *) (base + iindex + SIZE (M_REAL));
 162      STATUS (re) = INIT_MASK;
 163      VALUE (re) = v[2 * k];
 164      STATUS (im) = INIT_MASK;
 165      VALUE (im) = v[2 * k + 1];
 166      CHECK_COMPLEX (p, VALUE (re), VALUE (im));
 167    }
 168    PUSH_REF (p, desc);
 169  }
 170  
 171  //! @brief Push prime factorisation on the stack as [] INT.
 172  
 173  void genie_prime_factors (NODE_T * p)
 174  {
 175    gsl_error_handler_t *save_handler = gsl_set_error_handler (fft_error_handler);
 176    A68_INT n;
 177    A68_REF desc, row;
 178    A68_ARRAY arr;
 179    A68_TUPLE tup;
 180    int len, inc, iindex, k;
 181    BYTE_T *base;
 182    gsl_fft_complex_wavetable *wt;
 183    A68 (f_entry) = p;
 184    POP_OBJECT (p, &n, A68_INT);
 185    CHECK_INIT (p, INITIALISED (&n), M_INT);
 186    wt = gsl_fft_complex_wavetable_alloc ((size_t) (VALUE (&n)));
 187    len = (int) (NF (wt));
 188    NEW_ROW_1D (desc, row, arr, tup, M_ROW_INT, M_INT, len);
 189    base = DEREF (BYTE_T, &ARRAY (&arr));
 190    iindex = VECTOR_OFFSET (&arr, &tup);
 191    inc = SPAN (&tup) * ELEM_SIZE (&arr);
 192    for (k = 0; k < len; k++, iindex += inc) {
 193      A68_INT *x = (A68_INT *) (base + iindex);
 194      STATUS (x) = INIT_MASK;
 195      VALUE (x) = (int) ((FACTOR (wt))[k]);
 196    }
 197    gsl_fft_complex_wavetable_free (wt);
 198    PUSH_REF (p, desc);
 199    (void) gsl_set_error_handler (save_handler);
 200  }
 201  
 202  //! @brief PROC ([] COMPLEX) [] COMPLEX fft complex forward
 203  
 204  void genie_fft_complex_forward (NODE_T * p)
 205  {
 206    gsl_error_handler_t *save_handler = gsl_set_error_handler (fft_error_handler);
 207    int len, rc;
 208    REAL_T *data;
 209    gsl_fft_complex_wavetable *wt;
 210    gsl_fft_complex_workspace *ws;
 211    A68 (f_entry) = p;
 212    data = pop_array_complex (p, &len);
 213    fft_test_error (len == 0 ? GSL_EDOM : GSL_SUCCESS);
 214    wt = gsl_fft_complex_wavetable_alloc ((size_t) len);
 215    ws = gsl_fft_complex_workspace_alloc ((size_t) len);
 216    rc = gsl_fft_complex_forward (data, 1, (size_t) len, wt, ws);
 217    fft_test_error (rc);
 218    push_array_complex (p, data, len);
 219    gsl_fft_complex_wavetable_free (wt);
 220    gsl_fft_complex_workspace_free (ws);
 221    if (data != NO_REAL) {
 222      a68_free (data);
 223    }
 224    (void) gsl_set_error_handler (save_handler);
 225  }
 226  
 227  //! @brief PROC ([] COMPLEX) [] COMPLEX fft complex backward
 228  
 229  void genie_fft_complex_backward (NODE_T * p)
 230  {
 231    gsl_error_handler_t *save_handler = gsl_set_error_handler (fft_error_handler);
 232    int len, rc;
 233    REAL_T *data;
 234    gsl_fft_complex_wavetable *wt;
 235    gsl_fft_complex_workspace *ws;
 236    A68 (f_entry) = p;
 237    data = pop_array_complex (p, &len);
 238    fft_test_error (len == 0 ? GSL_EDOM : GSL_SUCCESS);
 239    wt = gsl_fft_complex_wavetable_alloc ((size_t) len);
 240    ws = gsl_fft_complex_workspace_alloc ((size_t) len);
 241    rc = gsl_fft_complex_backward (data, 1, (size_t) len, wt, ws);
 242    fft_test_error (rc);
 243    push_array_complex (p, data, len);
 244    gsl_fft_complex_wavetable_free (wt);
 245    gsl_fft_complex_workspace_free (ws);
 246    if (data != NO_REAL) {
 247      a68_free (data);
 248    }
 249    (void) gsl_set_error_handler (save_handler);
 250  }
 251  
 252  //! @brief PROC ([] COMPLEX) [] COMPLEX fft complex inverse
 253  
 254  void genie_fft_complex_inverse (NODE_T * p)
 255  {
 256    gsl_error_handler_t *save_handler = gsl_set_error_handler (fft_error_handler);
 257    int len, rc;
 258    REAL_T *data;
 259    gsl_fft_complex_wavetable *wt;
 260    gsl_fft_complex_workspace *ws;
 261    A68 (f_entry) = p;
 262    data = pop_array_complex (p, &len);
 263    fft_test_error (len == 0 ? GSL_EDOM : GSL_SUCCESS);
 264    wt = gsl_fft_complex_wavetable_alloc ((size_t) len);
 265    ws = gsl_fft_complex_workspace_alloc ((size_t) len);
 266    rc = gsl_fft_complex_inverse (data, 1, (size_t) len, wt, ws);
 267    fft_test_error (rc);
 268    push_array_complex (p, data, len);
 269    gsl_fft_complex_wavetable_free (wt);
 270    gsl_fft_complex_workspace_free (ws);
 271    if (data != NO_REAL) {
 272      a68_free (data);
 273    }
 274    (void) gsl_set_error_handler (save_handler);
 275  }
 276  
 277  //! @brief PROC ([] REAL) [] COMPLEX fft forward
 278  
 279  void genie_fft_forward (NODE_T * p)
 280  {
 281    gsl_error_handler_t *save_handler = gsl_set_error_handler (fft_error_handler);
 282    int len, rc;
 283    REAL_T *data;
 284    gsl_fft_complex_wavetable *wt;
 285    gsl_fft_complex_workspace *ws;
 286    A68 (f_entry) = p;
 287    data = pop_array_real (p, &len);
 288    fft_test_error (len == 0 ? GSL_EDOM : GSL_SUCCESS);
 289    wt = gsl_fft_complex_wavetable_alloc ((size_t) len);
 290    ws = gsl_fft_complex_workspace_alloc ((size_t) len);
 291    rc = gsl_fft_complex_forward (data, 1, (size_t) len, wt, ws);
 292    fft_test_error (rc);
 293    push_array_complex (p, data, len);
 294    gsl_fft_complex_wavetable_free (wt);
 295    gsl_fft_complex_workspace_free (ws);
 296    if (data != NO_REAL) {
 297      a68_free (data);
 298    }
 299    (void) gsl_set_error_handler (save_handler);
 300  }
 301  
 302  //! @brief PROC ([] COMPLEX) [] REAL fft backward
 303  
 304  void genie_fft_backward (NODE_T * p)
 305  {
 306    gsl_error_handler_t *save_handler = gsl_set_error_handler (fft_error_handler);
 307    int len, rc;
 308    REAL_T *data;
 309    gsl_fft_complex_wavetable *wt;
 310    gsl_fft_complex_workspace *ws;
 311    A68 (f_entry) = p;
 312    data = pop_array_complex (p, &len);
 313    fft_test_error (len == 0 ? GSL_EDOM : GSL_SUCCESS);
 314    wt = gsl_fft_complex_wavetable_alloc ((size_t) len);
 315    ws = gsl_fft_complex_workspace_alloc ((size_t) len);
 316    rc = gsl_fft_complex_backward (data, 1, (size_t) len, wt, ws);
 317    fft_test_error (rc);
 318    push_array_real (p, data, len);
 319    gsl_fft_complex_wavetable_free (wt);
 320    gsl_fft_complex_workspace_free (ws);
 321    if (data != NO_REAL) {
 322      a68_free (data);
 323    }
 324    (void) gsl_set_error_handler (save_handler);
 325  }
 326  
 327  //! @brief PROC ([] COMPLEX) [] REAL fft inverse
 328  
 329  void genie_fft_inverse (NODE_T * p)
 330  {
 331    gsl_error_handler_t *save_handler = gsl_set_error_handler (fft_error_handler);
 332    int len, rc;
 333    REAL_T *data;
 334    gsl_fft_complex_wavetable *wt;
 335    gsl_fft_complex_workspace *ws;
 336    A68 (f_entry) = p;
 337    data = pop_array_complex (p, &len);
 338    fft_test_error (len == 0 ? GSL_EDOM : GSL_SUCCESS);
 339    wt = gsl_fft_complex_wavetable_alloc ((size_t) len);
 340    ws = gsl_fft_complex_workspace_alloc ((size_t) len);
 341    rc = gsl_fft_complex_inverse (data, 1, (size_t) len, wt, ws);
 342    fft_test_error (rc);
 343    push_array_real (p, data, len);
 344    gsl_fft_complex_wavetable_free (wt);
 345    gsl_fft_complex_workspace_free (ws);
 346    if (data != NO_REAL) {
 347      a68_free (data);
 348    }
 349    (void) gsl_set_error_handler (save_handler);
 350  }
 351  
 352  #endif