htswanalysis/MACS/lib/gsl/gsl-1.11/wavelet/daubechies.c

   1 /* wavelet/daubechies.c
   2  *
   3  * Copyright (C) 2004 Ivo Alxneit
   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 3 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /*
  21  * Coefficients for Daubechies wavelets of extremal phase are from
  22  * I. Daubechies, "Orthonormal Bases of Compactly Supported Wavelets",
  23  * Communications on Pure and Applied Mathematics, 41 (1988) 909--996
  24  * (table 1).
  25  * Additional digits have been obtained using the Mathematica package
  26  * Daubechies.m by Tong Chen & Meng Xu available at
  27  * http://www.cwp.mines.edu/wavelets/.
  28  */
  29
  30 #include <config.h>
  31 #include <gsl/gsl_errno.h>
  32 #include <gsl/gsl_wavelet.h>
  33
  34 static const double h_4[4] = { 0.48296291314453414337487159986,
  35   0.83651630373780790557529378092,
  36   0.22414386804201338102597276224,
  37   -0.12940952255126038117444941881
  38 };
  39
  40 static const double g_4[4] = { -0.12940952255126038117444941881,
  41   -0.22414386804201338102597276224,
  42   0.83651630373780790557529378092,
  43   -0.48296291314453414337487159986
  44 };
  45
  46 static const double h_6[6] = { 0.33267055295008261599851158914,
  47   0.80689150931109257649449360409,
  48   0.45987750211849157009515194215,
  49   -0.13501102001025458869638990670,
  50   -0.08544127388202666169281916918,
  51   0.03522629188570953660274066472
  52 };
  53
  54 static const double g_6[6] = { 0.03522629188570953660274066472,
  55   0.08544127388202666169281916918,
  56   -0.13501102001025458869638990670,
  57   -0.45987750211849157009515194215,
  58   0.80689150931109257649449360409,
  59   -0.33267055295008261599851158914
  60 };
  61
  62 static const double h_8[8] = { 0.23037781330889650086329118304,
  63   0.71484657055291564708992195527,
  64   0.63088076792985890788171633830,
  65   -0.02798376941685985421141374718,
  66   -0.18703481171909308407957067279,
  67   0.03084138183556076362721936253,
  68   0.03288301166688519973540751355,
  69   -0.01059740178506903210488320852
  70 };
  71
  72 static const double g_8[8] = { -0.01059740178506903210488320852,
  73   -0.03288301166688519973540751355,
  74   0.03084138183556076362721936253,
  75   0.18703481171909308407957067279,
  76   -0.02798376941685985421141374718,
  77   -0.63088076792985890788171633830,
  78   0.71484657055291564708992195527,
  79   -0.23037781330889650086329118304
  80 };
  81
  82 static const double h_10[10] = { 0.16010239797419291448072374802,
  83   0.60382926979718967054011930653,
  84   0.72430852843777292772807124410,
  85   0.13842814590132073150539714634,
  86   -0.24229488706638203186257137947,
  87   -0.03224486958463837464847975506,
  88   0.07757149384004571352313048939,
  89   -0.00624149021279827427419051911,
  90   -0.01258075199908199946850973993,
  91   0.00333572528547377127799818342
  92 };
  93
  94 static const double g_10[10] = { 0.00333572528547377127799818342,
  95   0.01258075199908199946850973993,
  96   -0.00624149021279827427419051911,
  97   -0.07757149384004571352313048939,
  98   -0.03224486958463837464847975506,
  99   0.24229488706638203186257137947,
 100   0.13842814590132073150539714634,
 101   -0.72430852843777292772807124410,
 102   0.60382926979718967054011930653,
 103   -0.16010239797419291448072374802
 104 };
 105
 106 static const double h_12[12] = { 0.11154074335010946362132391724,
 107   0.49462389039845308567720417688,
 108   0.75113390802109535067893449844,
 109   0.31525035170919762908598965481,
 110   -0.22626469396543982007631450066,
 111   -0.12976686756726193556228960588,
 112   0.09750160558732304910234355254,
 113   0.02752286553030572862554083950,
 114   -0.03158203931748602956507908070,
 115   0.00055384220116149613925191840,
 116   0.00477725751094551063963597525,
 117   -0.00107730108530847956485262161
 118 };
 119
 120 static const double g_12[12] = { -0.00107730108530847956485262161,
 121   -0.00477725751094551063963597525,
 122   0.00055384220116149613925191840,
 123   0.03158203931748602956507908070,
 124   0.02752286553030572862554083950,
 125   -0.09750160558732304910234355254,
 126   -0.12976686756726193556228960588,
 127   0.22626469396543982007631450066,
 128   0.31525035170919762908598965481,
 129   -0.75113390802109535067893449844,
 130   0.49462389039845308567720417688,
 131   -0.11154074335010946362132391724
 132 };
 133
 134 static const double h_14[14] = { 0.07785205408500917901996352196,
 135   0.39653931948191730653900039094,
 136   0.72913209084623511991694307034,
 137   0.46978228740519312247159116097,
 138   -0.14390600392856497540506836221,
 139   -0.22403618499387498263814042023,
 140   0.07130921926683026475087657050,
 141   0.08061260915108307191292248036,
 142   -0.03802993693501441357959206160,
 143   -0.01657454163066688065410767489,
 144   0.01255099855609984061298988603,
 145   0.00042957797292136652113212912,
 146   -0.00180164070404749091526826291,
 147   0.00035371379997452024844629584
 148 };
 149
 150 static const double g_14[14] = { 0.00035371379997452024844629584,
 151   0.00180164070404749091526826291,
 152   0.00042957797292136652113212912,
 153   -0.01255099855609984061298988603,
 154   -0.01657454163066688065410767489,
 155   0.03802993693501441357959206160,
 156   0.08061260915108307191292248036,
 157   -0.07130921926683026475087657050,
 158   -0.22403618499387498263814042023,
 159   0.14390600392856497540506836221,
 160   0.46978228740519312247159116097,
 161   -0.72913209084623511991694307034,
 162   0.39653931948191730653900039094,
 163   -0.07785205408500917901996352196
 164 };
 165
 166 static const double h_16[16] = { 0.05441584224310400995500940520,
 167   0.31287159091429997065916237551,
 168   0.67563073629728980680780076705,
 169   0.58535468365420671277126552005,
 170   -0.01582910525634930566738054788,
 171   -0.28401554296154692651620313237,
 172   0.00047248457391328277036059001,
 173   0.12874742662047845885702928751,
 174   -0.01736930100180754616961614887,
 175   -0.04408825393079475150676372324,
 176   0.01398102791739828164872293057,
 177   0.00874609404740577671638274325,
 178   -0.00487035299345157431042218156,
 179   -0.00039174037337694704629808036,
 180   0.00067544940645056936636954757,
 181   -0.00011747678412476953373062823
 182 };
 183
 184 static const double g_16[16] = { -0.00011747678412476953373062823,
 185   -0.00067544940645056936636954757,
 186   -0.00039174037337694704629808036,
 187   0.00487035299345157431042218156,
 188   0.00874609404740577671638274325,
 189   -0.01398102791739828164872293057,
 190   -0.04408825393079475150676372324,
 191   0.01736930100180754616961614887,
 192   0.12874742662047845885702928751,
 193   -0.00047248457391328277036059001,
 194   -0.28401554296154692651620313237,
 195   0.01582910525634930566738054788,
 196   0.58535468365420671277126552005,
 197   -0.67563073629728980680780076705,
 198   0.31287159091429997065916237551,
 199   -0.05441584224310400995500940520
 200 };
 201
 202 static const double h_18[18] = { 0.03807794736387834658869765888,
 203   0.24383467461259035373204158165,
 204   0.60482312369011111190307686743,
 205   0.65728807805130053807821263905,
 206   0.13319738582500757619095494590,
 207   -0.29327378327917490880640319524,
 208   -0.09684078322297646051350813354,
 209   0.14854074933810638013507271751,
 210   0.03072568147933337921231740072,
 211   -0.06763282906132997367564227483,
 212   0.00025094711483145195758718975,
 213   0.02236166212367909720537378270,
 214   -0.00472320475775139727792570785,
 215   -0.00428150368246342983449679500,
 216   0.00184764688305622647661912949,
 217   0.00023038576352319596720521639,
 218   -0.00025196318894271013697498868,
 219   0.00003934732031627159948068988
 220 };
 221
 222 static const double g_18[18] = { 0.00003934732031627159948068988,
 223   0.00025196318894271013697498868,
 224   0.00023038576352319596720521639,
 225   -0.00184764688305622647661912949,
 226   -0.00428150368246342983449679500,
 227   0.00472320475775139727792570785,
 228   0.02236166212367909720537378270,
 229   -0.00025094711483145195758718975,
 230   -0.06763282906132997367564227483,
 231   -0.03072568147933337921231740072,
 232   0.14854074933810638013507271751,
 233   0.09684078322297646051350813354,
 234   -0.29327378327917490880640319524,
 235   -0.13319738582500757619095494590,
 236   0.65728807805130053807821263905,
 237   -0.60482312369011111190307686743,
 238   0.24383467461259035373204158165,
 239   -0.03807794736387834658869765888
 240 };
 241
 242 static const double h_20[20] = { 0.02667005790055555358661744877,
 243   0.18817680007769148902089297368,
 244   0.52720118893172558648174482796,
 245   0.68845903945360356574187178255,
 246   0.28117234366057746074872699845,
 247   -0.24984642432731537941610189792,
 248   -0.19594627437737704350429925432,
 249   0.12736934033579326008267723320,
 250   0.09305736460357235116035228984,
 251   -0.07139414716639708714533609308,
 252   -0.02945753682187581285828323760,
 253   0.03321267405934100173976365318,
 254   0.00360655356695616965542329142,
 255   -0.01073317548333057504431811411,
 256   0.00139535174705290116578931845,
 257   0.00199240529518505611715874224,
 258   -0.00068585669495971162656137098,
 259   -0.00011646685512928545095148097,
 260   0.00009358867032006959133405013,
 261   -0.00001326420289452124481243668
 262 };
 263
 264 static const double g_20[20] = { -0.00001326420289452124481243668,
 265   -0.00009358867032006959133405013,
 266   -0.00011646685512928545095148097,
 267   0.00068585669495971162656137098,
 268   0.00199240529518505611715874224,
 269   -0.00139535174705290116578931845,
 270   -0.01073317548333057504431811411,
 271   -0.00360655356695616965542329142,
 272   0.03321267405934100173976365318,
 273   0.02945753682187581285828323760,
 274   -0.07139414716639708714533609308,
 275   -0.09305736460357235116035228984,
 276   0.12736934033579326008267723320,
 277   0.19594627437737704350429925432,
 278   -0.24984642432731537941610189792,
 279   -0.28117234366057746074872699845,
 280   0.68845903945360356574187178255,
 281   -0.52720118893172558648174482796,
 282   0.18817680007769148902089297368,
 283   -0.02667005790055555358661744877
 284 };
 285
 286 static int
 287 daubechies_init (const double **h1, const double **g1, const double **h2,
 288                  const double **g2, size_t * nc, size_t * offset,
 289                  size_t member)
 290 {
 291   switch (member)
 292     {
 293     case 4:
 294       *h1 = h_4;
 295       *g1 = g_4;
 296       *h2 = h_4;
 297       *g2 = g_4;
 298       break;
 299
 300     case 6:
 301       *h1 = h_6;
 302       *g1 = g_6;
 303       *h2 = h_6;
 304       *g2 = g_6;
 305       break;
 306
 307     case 8:
 308       *h1 = h_8;
 309       *g1 = g_8;
 310       *h2 = h_8;
 311       *g2 = g_8;
 312       break;
 313
 314     case 10:
 315       *h1 = h_10;
 316       *g1 = g_10;
 317       *h2 = h_10;
 318       *g2 = g_10;
 319       break;
 320
 321     case 12:
 322       *h1 = h_12;
 323       *g1 = g_12;
 324       *h2 = h_12;
 325       *g2 = g_12;
 326       break;
 327
 328     case 14:
 329       *h1 = h_14;
 330       *g1 = g_14;
 331       *h2 = h_14;
 332       *g2 = g_14;
 333       break;
 334
 335     case 16:
 336       *h1 = h_16;
 337       *g1 = g_16;
 338       *h2 = h_16;
 339       *g2 = g_16;
 340       break;
 341
 342     case 18:
 343       *h1 = h_18;
 344       *g1 = g_18;
 345       *h2 = h_18;
 346       *g2 = g_18;
 347       break;
 348
 349     case 20:
 350       *h1 = h_20;
 351       *g1 = g_20;
 352       *h2 = h_20;
 353       *g2 = g_20;
 354       break;
 355
 356     default:
 357       return GSL_FAILURE;
 358     }
 359
 360   *nc = member;
 361   *offset = 0;
 362
 363   return GSL_SUCCESS;
 364 }
 365
 366 static int
 367 daubechies_centered_init (const double **h1, const double **g1,
 368                           const double **h2, const double **g2, size_t * nc,
 369                           size_t * offset, size_t member)
 370 {
 371   switch (member)
 372     {
 373     case 4:
 374       *h1 = h_4;
 375       *g1 = g_4;
 376       *h2 = h_4;
 377       *g2 = g_4;
 378       break;
 379
 380     case 6:
 381       *h1 = h_6;
 382       *g1 = g_6;
 383       *h2 = h_6;
 384       *g2 = g_6;
 385       break;
 386
 387     case 8:
 388       *h1 = h_8;
 389       *g1 = g_8;
 390       *h2 = h_8;
 391       *g2 = g_8;
 392       break;
 393
 394     case 10:
 395       *h1 = h_10;
 396       *g1 = g_10;
 397       *h2 = h_10;
 398       *g2 = g_10;
 399       break;
 400
 401     case 12:
 402       *h1 = h_12;
 403       *g1 = g_12;
 404       *h2 = h_12;
 405       *g2 = g_12;
 406       break;
 407
 408     case 14:
 409       *h1 = h_14;
 410       *g1 = g_14;
 411       *h2 = h_14;
 412       *g2 = g_14;
 413       break;
 414
 415     case 16:
 416       *h1 = h_16;
 417       *g1 = g_16;
 418       *h2 = h_16;
 419       *g2 = g_16;
 420       break;
 421
 422     case 18:
 423       *h1 = h_18;
 424       *g1 = g_18;
 425       *h2 = h_18;
 426       *g2 = g_18;
 427       break;
 428
 429     case 20:
 430       *h1 = h_20;
 431       *g1 = g_20;
 432       *h2 = h_20;
 433       *g2 = g_20;
 434       break;
 435
 436     default:
 437       return GSL_FAILURE;
 438     }
 439
 440   *nc = member;
 441   *offset = (member >> 1);
 442
 443   return GSL_SUCCESS;
 444 }
 445
 446 static const gsl_wavelet_type daubechies_type = {
 447   "daubechies",
 448   &daubechies_init
 449 };
 450
 451 static const gsl_wavelet_type daubechies_centered_type = {
 452   "daubechies-centered",
 453   &daubechies_centered_init
 454 };
 455
 456 const gsl_wavelet_type *gsl_wavelet_daubechies = &daubechies_type;
 457 const gsl_wavelet_type *gsl_wavelet_daubechies_centered =
 458   &daubechies_centered_type;