home *** CD-ROM | disk | FTP | other *** search
/ Los Alamos National Laboratory / LANL_CD.ISO / software / compres / src / wavelets.c < prev    next >
C/C++ Source or Header  |  1992-02-11  |  6KB  |  293 lines
  1. /****************************************************************
  2.  
  3. COPYRIGHT (C) 1992 UNIVERSITY OF CALIFORNIA
  4.  
  5. ***************************************************************/
  6.  
  7. #include <math.h>
  8.  
  9. b4_const(modx2sq, rx2, sqrt2, x1)
  10. double *modx2sq, *rx2, *sqrt2, *x1;
  11. {
  12.     double a, b;
  13.  
  14.     *sqrt2 = sqrt(2.0);
  15.     a = (-14.+(63./sqrt(15.)))/1080.;
  16.     a = exp(log(a)/3.);
  17.     b = (14.+(63./sqrt(15.)))/1080.;
  18.     b = -exp(log(b)/3.);
  19.     *x1 = a + b - (1./6.);
  20.     *rx2 = -((a+b)/2.)-(1./6.);
  21.     *modx2sq = (*rx2 * *rx2)+3.*(((a-b)*(a-b)/4.));
  22. }
  23.  
  24.  
  25. /*  Family b4_4. # cos(z/2) factors with Pa, Pb: Na=4, Nb=4.  */
  26.  
  27. b4_4al7(h)
  28. double *h;
  29. {
  30.     double modx2sq, rx2, sqrt2, x1;
  31.  
  32.     b4_const(&modx2sq, &rx2, &sqrt2, &x1);
  33.     h[0] = sqrt2*(6.*x1-1.)/(16.*x1);
  34.     h[1] = sqrt2*(16.*x1-1.)/(64.*x1);
  35.     h[2] = sqrt2*(2.*x1+1.)/(32.*x1);
  36.     h[3] = sqrt2/(64.*x1);
  37. }
  38.  
  39. b4_4ah9(g)
  40. double *g;
  41. {
  42.     int n, nn;
  43.     double *h;
  44.  
  45.     h = (double*)malloc(5*sizeof(double));
  46.     b4_4bl9(h);
  47.     g -= 1;
  48.     nn = -1;
  49.     for(n = 1; n <= 5; n++) {
  50.         g[n]=nn*h[n-1];
  51.         nn *= -1;
  52.     }
  53. }
  54.  
  55. b4_4bl9(h)
  56. double *h;
  57. {
  58.     double modx2sq, rx2, sqrt2, x1;
  59.  
  60.     b4_const(&modx2sq,&rx2,&sqrt2,&x1);
  61.     h[0] = -5.*sqrt2*x1*(48.*modx2sq-16.*rx2+3.)/32.;
  62.     h[1] = -5.*sqrt2*x1*(8.*modx2sq-rx2)/8.;
  63.     h[2] = -5.*sqrt2*x1*(4.*modx2sq+4.*rx2-1.)/16.;
  64.     h[3] = -5.*sqrt2*x1*rx2/8.;
  65.     h[4] = -5.*sqrt2*x1/64.;
  66. }
  67.  
  68.  
  69. b4_4bh7(g)
  70. double *g;
  71. {
  72.     int n, nn;
  73.     double *h;
  74.  
  75.     g--;
  76.     h = (double*)malloc(4*sizeof(double));
  77.     b4_4al7(h);
  78.     nn = -1;
  79.     for(n = 1; n <= 4; n++) {
  80.         g[n]=nn*h[n-1];
  81.         nn *= -1;
  82.     }
  83. }
  84.  
  85. /****************************************************************/
  86.  
  87. b5_const(a0,a1,b0,b1,c,sqrt2)
  88. double *a0, *a1, *b0, *b1, *c, *sqrt2;
  89. {
  90.     double alpha, beta, xgamma, r, s, a;
  91.  
  92.     *sqrt2=sqrt(2.);
  93.     a = cos(acos(-sqrt(15.)/8.)/3.);
  94.     a = (a*8.*sqrt(0.6)/7.)-(9./28.);
  95.     alpha = sqrt(a);
  96.     beta = 27./56.;
  97.     xgamma = -15./56.;
  98.     r = (beta+a-(xgamma/alpha))/2.;
  99.     s = (beta+a+(xgamma/alpha))/2.;
  100.     *a0 = alpha-2.5;
  101.     *a1 = -alpha-2.5;
  102.     *b0 = (25./16.)-(1.25*alpha)+r;
  103.     *b1 = (25./16.)+(1.25*alpha)+s;
  104.     *c = *a0 + *b0 + 1.;
  105. }
  106.  
  107. b5_6al11(h)
  108. double *h;
  109. {
  110.     double a0, a1, b0, b1, c, sqrt2;
  111.  
  112.     b5_const(&a0,&a1,&b0,&b1,&c,&sqrt2);
  113.     h[0] = sqrt2*35.*c*(15.*a1+20.*b1+13.)/512.;
  114.     h[1] = sqrt2*35.*c*(26.*a1+30.*b1+23)/1024.;
  115.     h[2] = sqrt2*35.*c*(4.*a1+3.*b1+4.)/256.;
  116.     h[3] = sqrt2*35.*c*(12.*a1+4.*b1+17.)/2048.;
  117.     h[4] = sqrt2*35.*c*(a1+3.)/1024.;
  118.     h[5] = sqrt2*35.*c/2048.;
  119. }
  120.  
  121. b5_6ah9(g)
  122. double *g;
  123. {
  124.     int n, nn;
  125.     double *h;
  126.  
  127.     h = (double*)malloc(5*sizeof(double));
  128.     b5_6bl9(h);
  129.     g--;
  130.     nn = -1;
  131.     for(n = 1; n <= 5; n++) {
  132.         g[n]=nn*h[n-1];
  133.         nn *= -1;
  134.     }
  135. }
  136.  
  137. b5_6bl9(h)
  138. double *h;
  139. {
  140.     double a0, a1, b0, b1, c, sqrt2;
  141.  
  142.     b5_const(&a0,&a1,&b0,&b1,&c,&sqrt2);
  143.     h[0] = sqrt2*(8.*a0+12.*b0+7.)/(32.*c);
  144.     h[1] = sqrt2*(7.*a0+8.*b0+6.)/(32.*c);
  145.     h[2] = sqrt2*(2.*a0+b0+2.)/(16.*c);
  146.     h[3] = sqrt2*(a0+2.)/(32.*c);
  147.     h[4] = sqrt2/(64.*c);
  148. }
  149.  
  150. b5_6bh11(g)
  151. double *g;
  152. {
  153.     int n, nn;
  154.     double *h;
  155.  
  156.     h = (double*)malloc(6*sizeof(double));
  157.     b5_6al11(h);
  158.     g--;
  159.     nn = -1;
  160.     for(n = 1; n <= 6; n++) {
  161.         g[n]=nn*h[n-1];
  162.         nn *= -1;
  163.     }
  164. }
  165.  
  166. /****************************************************************/
  167.  
  168. b5_4al9(h)
  169. double *h;
  170. {
  171.     double a0, a1, b0, b1, c, sqrt2;
  172.  
  173.     b5_const(&a0,&a1,&b0,&b1,&c,&sqrt2);
  174.     h[0] = sqrt2*35.*c*(8.*a1+12.*b1+7.)/256.;
  175.     h[1] = sqrt2*35.*c*(7.*a1+8.*b1+6.)/256.;
  176.     h[2] = sqrt2*35.*c*(2.*a1+b1+2.)/128.;
  177.     h[3] = sqrt2*35.*c*(a1+2.)/256.;
  178.     h[4] = sqrt2*35.*c/512.;
  179. }
  180.  
  181. b5_4ah11(g)
  182. double *g;
  183. {
  184.     int n, nn;
  185.     double *h;
  186.  
  187.     h = (double*)malloc(6*sizeof(double));
  188.     b5_4bl11(h);
  189.     g--;
  190.     nn = -1;
  191.     for(n = 1; n <= 6; n++) {
  192.         g[n]=nn*h[n-1];
  193.         nn *= -1;
  194.     }
  195. }
  196.  
  197. b5_4bl11(h)
  198. double *h;
  199. {
  200.     double a0, a1, b0, b1, c, sqrt2;
  201.  
  202.     b5_const(&a0,&a1,&b0,&b1,&c,&sqrt2);
  203.     h[0] = sqrt2*(15.*a0+20.*b0+13.)/(64.*c);
  204.     h[1] = sqrt2*(26.*a0+30.*b0+23)/(128.*c);
  205.     h[2] = sqrt2*(4.*a0+3.*b0+4.)/(32.*c);
  206.     h[3] = sqrt2*(12.*a0+4.*b0+17.)/(256.*c);
  207.     h[4] = sqrt2*(a0+3.)/(128.*c);
  208.     h[5] = sqrt2/(256.*c);
  209. }
  210.  
  211. b5_4bh9(g)
  212. double *g;
  213. {
  214.     int n, nn;
  215.     double *h;
  216.  
  217.     h = (double*)malloc(5*sizeof(double));
  218.     b5_4al9(h);
  219.     nn = -1;
  220.     g--;
  221.     for(n = 1; n <= 5; n++) {
  222.         g[n]=nn*h[n-1];
  223.         nn *= -1;
  224.     }
  225. }
  226.  
  227. /****************************************************************/
  228.  
  229. b5_5al10(h)
  230. double *h;
  231. {
  232.     double a0, a1, b0, b1, c, sqrt2;
  233.  
  234.     b5_const(&a0, &a1, &b0, &b1, &c, &sqrt2);
  235.     h--;
  236.  
  237.     h[1] = sqrt2*35.*c*(15.*a1+20.*b1+13.)/512.;
  238.     h[2] = sqrt2*35.*c*(11.*a1+10.*b1+10.)/512.;
  239.     h[3] = sqrt2*35.*c*(5.*a1+2.*b1+6.)/512.;
  240.     h[4] = sqrt2*35.*c*(2.*a1+5.)/1024.;
  241.     h[5] = sqrt2*35.*c/1024.;
  242. }
  243.  
  244. b5_5ah10(g)
  245. double *g;
  246. {
  247.     int n, nn;
  248.     double *h;
  249.  
  250.     h = (double*)malloc(5*sizeof(double));
  251.     b5_5bl10(h);
  252.     h--;
  253.  
  254.     nn = -1;
  255.     g--;
  256.     for(n = 1; n <= 5; n++) {
  257.     g[n]=nn*h[n];
  258.         nn *= -1;
  259.     }
  260. }
  261.  
  262. b5_5bl10(h)
  263. double *h;
  264. {
  265.     double a0, a1, b0, b1, c, sqrt2;
  266.  
  267.     b5_const(&a0, &a1, &b0, &b1, &c, &sqrt2);
  268.     h--;
  269.     h[1] = sqrt2*(15.*a0+20.*b0+13.)/(64.*c);
  270.     h[2] = sqrt2*(11.*a0+10.*b0+10.)/(64.*c);
  271.     h[3] = sqrt2*(5.*a0+2.*b0+6.)/(64.*c);
  272.     h[4] = sqrt2*(2.*a0+5.)/(128.*c);
  273.     h[5] = sqrt2/(128.*c);
  274. }
  275.  
  276. b5_5bh10(g)
  277. double *g;
  278. {
  279.     int n, nn;
  280.     double *h;
  281.  
  282.     h = (double*)malloc(5*sizeof(double));
  283.     b5_5al10(h);
  284.  
  285.     h--;
  286.     g--;
  287.     nn = -1;
  288.     for(n = 1; n <= 5; n++) {
  289.     g[n]=nn*h[n];
  290.         nn *= -1;
  291.     }
  292. }
  293.