gsl_wavelet__daubechies.c 12 KB

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