src/mf.c

   1 #include "a/math.h"
   2 #include "a/mf.h"
   3
   4 a_float a_mf_gauss(a_float x, a_float sigma, a_float c)
   5 {
   6     return a_float_exp(a_float_pow((x - c) / sigma, 2) / -2);
   7 }
   8
   9 a_float a_mf_gauss2(a_float x, a_float sigma1, a_float c1, a_float sigma2, a_float c2)
  10 {
  11     if (x < c1)
  12     {
  13         x = a_mf_gauss(x, sigma1, c1);
  14     }
  15     else if (x > c2)
  16     {
  17         x = a_mf_gauss(x, sigma2, c2);
  18     }
  19     else // c1 <= x <= c2
  20     {
  21         x = 1;
  22     }
  23     return x;
  24 }
  25
  26 a_float a_mf_gbell(a_float x, a_float a, a_float b, a_float c)
  27 {
  28     return 1 / (a_float_pow(a_float_abs((x - c) / a), 2 * b) + 1);
  29 }
  30
  31 a_float a_mf_sig(a_float x, a_float a, a_float c)
  32 {
  33     return 1 / (a_float_exp((c - x) * a) + 1);
  34 }
  35
  36 a_float a_mf_dsig(a_float x, a_float a1, a_float c1, a_float a2, a_float c2)
  37 {
  38     return a_mf_sig(x, a1, c1) - a_mf_sig(x, a2, c2);
  39 }
  40
  41 a_float a_mf_psig(a_float x, a_float a1, a_float c1, a_float a2, a_float c2)
  42 {
  43     return a_mf_sig(x, a1, c1) * a_mf_sig(x, a2, c2);
  44 }
  45
  46 a_float a_mf_trap(a_float x, a_float a, a_float b, a_float c, a_float d)
  47 {
  48     if (x < b)
  49     {
  50         if (x > a) // a < x <= b
  51         {
  52             x = (x - a) / (b - a);
  53         }
  54         else // x <= a
  55         {
  56             x = 0;
  57         }
  58     }
  59     else if (x > c)
  60     {
  61         if (x < d) // c <= x < d
  62         {
  63             x = (d - x) / (d - c);
  64         }
  65         else // d <= x
  66         {
  67             x = 0;
  68         }
  69     }
  70     else // b <= x <= c
  71     {
  72         x = 1;
  73     }
  74     return x;
  75 }
  76
  77 a_float a_mf_tri(a_float x, a_float a, a_float b, a_float c)
  78 {
  79     if (x < b)
  80     {
  81         if (x > a) // a < x <= b
  82         {
  83             x = (x - a) / (b - a);
  84         }
  85         else // x <= a
  86         {
  87             x = 0;
  88         }
  89     }
  90     else
  91     {
  92         if (x < c) // b <= x < c
  93         {
  94             x = (c - x) / (c - b);
  95         }
  96         else // c <= x
  97         {
  98             x = 0;
  99         }
 100     }
 101     return x;
 102 }
 103
 104 a_float a_mf_lins(a_float x, a_float a, a_float b)
 105 {
 106     if (x < a)
 107     {
 108         x = 0;
 109     }
 110     else if (x > b)
 111     {
 112         x = 1;
 113     }
 114     else // a <= x <= b
 115     {
 116         x = (x - a) / (b - a);
 117     }
 118     return x;
 119 }
 120
 121 a_float a_mf_linz(a_float x, a_float a, a_float b)
 122 {
 123     if (x < a)
 124     {
 125         x = 1;
 126     }
 127     else if (x > b)
 128     {
 129         x = 0;
 130     }
 131     else // a <= x <= b
 132     {
 133         x = (b - x) / (b - a);
 134     }
 135     return x;
 136 }
 137
 138 a_float a_mf_s(a_float x, a_float a, a_float b)
 139 {
 140     if (x > (a + b) / 2)
 141     {
 142         if (x < b)
 143         {
 144             x = 1 - 2 * a_float_pow((b - x) / (b - a), 2);
 145         }
 146         else // x >= b
 147         {
 148             x = 1;
 149         }
 150     }
 151     else // x <= (a+b)/2
 152     {
 153         if (x > a)
 154         {
 155             x = 2 * a_float_pow((x - a) / (b - a), 2);
 156         }
 157         else // x <= a
 158         {
 159             x = 0;
 160         }
 161     }
 162     return x;
 163 }
 164
 165 a_float a_mf_z(a_float x, a_float a, a_float b)
 166 {
 167     if (x < (a + b) / 2)
 168     {
 169         if (x > a)
 170         {
 171             x = 1 - 2 * a_float_pow((x - a) / (b - a), 2);
 172         }
 173         else // x <= a
 174         {
 175             x = 1;
 176         }
 177     }
 178     else // x >= (a+b)/2
 179     {
 180         if (x < b)
 181         {
 182             x = 2 * a_float_pow((b - x) / (b - a), 2);
 183         }
 184         else // x >= b
 185         {
 186             x = 0;
 187         }
 188     }
 189     return x;
 190 }
 191
 192 a_float a_mf_pi(a_float x, a_float a, a_float b, a_float c, a_float d)
 193 {
 194     if (x < b)
 195     {
 196         x = a_mf_s(x, a, b);
 197     }
 198     else if (x > c)
 199     {
 200         x = a_mf_z(x, c, d);
 201     }
 202     else
 203     {
 204         x = 1;
 205     }
 206     return x;
 207 }
 208
 209 a_float a_mf(unsigned int e, a_float x, a_float const *a)
 210 {
 211     switch (e)
 212     {
 213     case A_MF_PI:
 214         return a_mf_pi(x, a[0], a[1], a[2], a[3]);
 215     case A_MF_Z:
 216         return a_mf_z(x, a[0], a[1]);
 217     case A_MF_S:
 218         return a_mf_s(x, a[0], a[1]);
 219     case A_MF_LINZ:
 220         return a_mf_linz(x, a[0], a[1]);
 221     case A_MF_LINS:
 222         return a_mf_lins(x, a[0], a[1]);
 223     case A_MF_TRI:
 224         return a_mf_tri(x, a[0], a[1], a[2]);
 225     case A_MF_TRAP:
 226         return a_mf_trap(x, a[0], a[1], a[2], a[3]);
 227     case A_MF_PSIG:
 228         return a_mf_psig(x, a[0], a[1], a[2], a[3]);
 229     case A_MF_DSIG:
 230         return a_mf_dsig(x, a[0], a[1], a[2], a[3]);
 231     case A_MF_SIG:
 232         return a_mf_sig(x, a[0], a[1]);
 233     case A_MF_GBELL:
 234         return a_mf_gbell(x, a[0], a[1], a[2]);
 235     case A_MF_GAUSS2:
 236         return a_mf_gauss2(x, a[0], a[1], a[2], a[3]);
 237     case A_MF_GAUSS:
 238         return a_mf_gauss(x, a[0], a[1]);
 239     case A_MF_NUL:
 240     default:
 241         return 0;
 242     }
 243 }