1 /* Resampling library

     2  * Copyright (C) <2001> David A. Schleef <ds@schleef.org>

     3  *

     4  * This library is free software; you can redistribute it and/or

     5  * modify it under the terms of the GNU Library General Public

     6  * License as published by the Free Software Foundation; either

     7  * version 2 of the License, or any later version.

     8  *

     9  * This library is distributed in the hope that it will be useful,

    10  * but WITHOUT ANY WARRANTY; without even the implied warranty of

    11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU

    12  * Library General Public License for more details.

    13  *

    14  * You should have received a copy of the GNU Library General Public

    15  * License along with this library; if not, write to the

    16  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,

    17  * Boston, MA 02111-1307, USA.

    18  */

    19 

    20 #ifdef HAVE_CONFIG_H

    21 #ifndef __SYMBIAN32__

    22 #include <config.h>

    23 #else

    24 #include "config.h"

    25 #endif

    26 #endif

    27 

    28 #include <string.h>

    29 #include <math.h>

    30 #include <stdio.h>

    31 #include <stdlib.h>

    32 

    33 #include "functable.h"

    34 #include "debug.h"

    35 

    36 

    37 #ifdef __SYMBIAN32__

    38 EXPORT_C

    39 #endif

    40 

    41 void

    42 functable_func_sinc (double *fx, double *dfx, double x, void *closure)

    43 {

    44   if (x == 0) {

    45     *fx = 1;

    46     *dfx = 0;

    47     return;

    48   }

    49 

    50   *fx = sin (x) / x;

    51   *dfx = (cos (x) - sin (x) / x) / x;

    52 }

    53 #ifdef __SYMBIAN32__

    54 EXPORT_C

    55 #endif

    56 

    57 

    58 void

    59 functable_func_boxcar (double *fx, double *dfx, double x, void *closure)

    60 {

    61   double width = *(double *) closure;

    62 

    63   if (x < width && x > -width) {

    64     *fx = 1;

    65   } else {

    66     *fx = 0;

    67   }

    68   *dfx = 0;

    69 }

    70 #ifdef __SYMBIAN32__

    71 EXPORT_C

    72 #endif

    73 

    74 

    75 void

    76 functable_func_hanning (double *fx, double *dfx, double x, void *closure)

    77 {

    78   double width = *(double *) closure;

    79 

    80   if (x < width && x > -width) {

    81     x /= width;

    82     *fx = (1 - x * x) * (1 - x * x);

    83     *dfx = -2 * 2 * x / width * (1 - x * x);

    84   } else {

    85     *fx = 0;

    86     *dfx = 0;

    87   }

    88 }

    89 #ifdef __SYMBIAN32__

    90 EXPORT_C

    91 #endif

    92 

    93 

    94 

    95 Functable *

    96 functable_new (void)

    97 {

    98   Functable *ft;

    99 

   100   ft = malloc (sizeof (Functable));

   101   memset (ft, 0, sizeof (Functable));

   102 

   103   return ft;

   104 }

   105 #ifdef __SYMBIAN32__

   106 EXPORT_C

   107 #endif

   108 

   109 

   110 void

   111 functable_free (Functable * ft)

   112 {

   113   free (ft);

   114 }

   115 #ifdef __SYMBIAN32__

   116 EXPORT_C

   117 #endif

   118 

   119 

   120 void

   121 functable_set_length (Functable * t, int length)

   122 {

   123   t->length = length;

   124 }

   125 #ifdef __SYMBIAN32__

   126 EXPORT_C

   127 #endif

   128 

   129 

   130 void

   131 functable_set_offset (Functable * t, double offset)

   132 {

   133   t->offset = offset;

   134 }

   135 #ifdef __SYMBIAN32__

   136 EXPORT_C

   137 #endif

   138 

   139 

   140 void

   141 functable_set_multiplier (Functable * t, double multiplier)

   142 {

   143   t->multiplier = multiplier;

   144 }

   145 #ifdef __SYMBIAN32__

   146 EXPORT_C

   147 #endif

   148 

   149 

   150 void

   151 functable_calculate (Functable * t, FunctableFunc func, void *closure)

   152 {

   153   int i;

   154   double x;

   155 

   156   if (t->fx)

   157     free (t->fx);

   158   if (t->dfx)

   159     free (t->dfx);

   160 

   161   t->fx = malloc (sizeof (double) * (t->length + 1));

   162   t->dfx = malloc (sizeof (double) * (t->length + 1));

   163 

   164   t->inv_multiplier = 1.0 / t->multiplier;

   165 

   166   for (i = 0; i < t->length + 1; i++) {

   167     x = t->offset + t->multiplier * i;

   168 

   169     func (&t->fx[i], &t->dfx[i], x, closure);

   170   }

   171 }

   172 #ifdef __SYMBIAN32__

   173 EXPORT_C

   174 #endif

   175 

   176 

   177 void

   178 functable_calculate_multiply (Functable * t, FunctableFunc func, void *closure)

   179 {

   180   int i;

   181   double x;

   182 

   183   for (i = 0; i < t->length + 1; i++) {

   184     double afx, adfx, bfx, bdfx;

   185 

   186     afx = t->fx[i];

   187     adfx = t->dfx[i];

   188     x = t->offset + t->multiplier * i;

   189     func (&bfx, &bdfx, x, closure);

   190     t->fx[i] = afx * bfx;

   191     t->dfx[i] = afx * bdfx + adfx * bfx;

   192   }

   193 

   194 }

   195 #ifdef __SYMBIAN32__

   196 EXPORT_C

   197 #endif

   198 

   199 

   200 double

   201 functable_evaluate (Functable * t, double x)

   202 {

   203   int i;

   204   double f0, f1, w0, w1;

   205   double x2, x3;

   206   double w;

   207 

   208   if (x < t->offset || x > (t->offset + t->length * t->multiplier)) {

   209     RESAMPLE_DEBUG ("x out of range %g", x);

   210   }

   211 

   212   x -= t->offset;

   213   x *= t->inv_multiplier;

   214   i = floor (x);

   215   x -= i;

   216 

   217   x2 = x * x;

   218   x3 = x2 * x;

   219 

   220   f1 = 3 * x2 - 2 * x3;

   221   f0 = 1 - f1;

   222   w0 = (x - 2 * x2 + x3) * t->multiplier;

   223   w1 = (-x2 + x3) * t->multiplier;

   224 

   225   w = t->fx[i] * f0 + t->fx[i + 1] * f1 + t->dfx[i] * w0 + t->dfx[i + 1] * w1;

   226 

   227   /*w = t->fx[i] * (1-x) + t->fx[i+1] * x; */

   228 

   229   return w;

   230 }

   231 #ifdef __SYMBIAN32__

   232 EXPORT_C

   233 #endif

   234 

   235 

   236 

   237 double

   238 functable_fir (Functable * t, double x, int n, double *data, int len)

   239 {

   240   int i, j;

   241   double f0, f1, w0, w1;

   242   double x2, x3;

   243   double w;

   244   double sum;

   245 

   246   x -= t->offset;

   247   x /= t->multiplier;

   248   i = floor (x);

   249   x -= i;

   250 

   251   x2 = x * x;

   252   x3 = x2 * x;

   253 

   254   f1 = 3 * x2 - 2 * x3;

   255   f0 = 1 - f1;

   256   w0 = (x - 2 * x2 + x3) * t->multiplier;

   257   w1 = (-x2 + x3) * t->multiplier;

   258 

   259   sum = 0;

   260   for (j = 0; j < len; j++) {

   261     w = t->fx[i] * f0 + t->fx[i + 1] * f1 + t->dfx[i] * w0 + t->dfx[i + 1] * w1;

   262     sum += data[j * 2] * w;

   263     i += n;

   264   }

   265 

   266   return sum;

   267 }

   268 #ifdef __SYMBIAN32__

   269 EXPORT_C

   270 #endif

   271 

   272 

   273 void

   274 functable_fir2 (Functable * t, double *r0, double *r1, double x,

   275     int n, double *data, int len)

   276 {

   277   int i, j;

   278   double f0, f1, w0, w1;

   279   double x2, x3;

   280   double w;

   281   double sum0, sum1;

   282   double floor_x;

   283 

   284   x -= t->offset;

   285   x *= t->inv_multiplier;

   286   floor_x = floor (x);

   287   i = floor_x;

   288   x -= floor_x;

   289 

   290   x2 = x * x;

   291   x3 = x2 * x;

   292 

   293   f1 = 3 * x2 - 2 * x3;

   294   f0 = 1 - f1;

   295   w0 = (x - 2 * x2 + x3) * t->multiplier;

   296   w1 = (-x2 + x3) * t->multiplier;

   297 

   298   sum0 = 0;

   299   sum1 = 0;

   300   for (j = 0; j < len; j++) {

   301     w = t->fx[i] * f0 + t->fx[i + 1] * f1 + t->dfx[i] * w0 + t->dfx[i + 1] * w1;

   302     sum0 += data[j * 2] * w;

   303     sum1 += data[j * 2 + 1] * w;

   304     i += n;

   305   }

   306 

   307   *r0 = sum0;

   308   *r1 = sum1;

   309 }