coif4.cpp 7.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313
  1. /*****************************************************************/
  2. /****** Coifman wavelet and scaling function ( L = 4 ) *******/
  3. /*****************************************************************/
  4. //#include "stdafx.h"
  5. #include <math.h>
  6. #include <stdlib.h>
  7. #include "vector.h"
  8. #include "matrix.h"
  9. #include "myfunc.h"
  10. #include "coif4.h"
  11. /*****************************************************************/
  12. void getcoif4( Vector& sxc, Vector& scc, Vector& wxc, Vector& wlc, int nit ) {
  13. const int nh = int((scc.dim()-1)/11.0);
  14. const int nd = 11*nh+1;
  15. const int ns = 4*nh;
  16. const int nw = 5*nh;
  17. int i, j, k, m;
  18. double h, x, sum;
  19. Vector scal(nd,0.0), wavl(nd,0.0), ftmp(nd,0.0);
  20. /**********************************************************/
  21. h = 1.0/double(nh);
  22. /********************* filter weghts **********************/
  23. double fw[] = { 0.011587596739,-0.029320137980,-0.047639590310,
  24. 0.273021046535, 0.574682393857, 0.294867193696,
  25. -0.054085607092,-0.042026480461, 0.016744410163,
  26. 0.003967883613,-0.001289203356,-0.000509505399 };
  27. /*************** starting of the iteration ****************/
  28. for(i = 0; i <= nh; i++) {
  29. double ih = double (i) * h;
  30. scal[i+ns-nh] = ih;
  31. // for(i = 1; i <= nh; i++)
  32. if ( i > 0 )
  33. scal[i+ns] = 1.0 - ih;
  34. }
  35. /******** iterative procedure for scaling function ********/
  36. for( i = 0; i <= nit; i++ ) {
  37. for( j = -ns; j < nd-ns; j++ ) {
  38. sum = 0.0;
  39. for( k = -4; k <= 7; k++ ) {
  40. m = 2*j-k*nh+ns;
  41. if( m >= 0 && m < nd )
  42. sum += fw[k+4]*scal[m];
  43. }
  44. ftmp[j+ns] = sum*2.0;
  45. }
  46. scal = ftmp;
  47. }
  48. /************************* wavelet ************************/
  49. for( j = -nw; j < nd-nw; j++ ) {
  50. sum = 0.0;
  51. for( k = -6; k <= 5; k++ ) {
  52. m = 2*j-k*nh+ns;
  53. if( m >= 0 && m < nd )
  54. sum += fw[1-k+4]*scal[m]*pow( -1.0, double(k+1) );
  55. }
  56. wavl[j+nw] = sum*2.0;
  57. }
  58. /**************** to print obtained data *****************/
  59. for( i = 0; i < nd; i++ ) {
  60. x = double(i-ns)*h;
  61. sxc[i] = x;
  62. scc[i] = scal[i];
  63. // }
  64. //
  65. // for( i = 0; i < nd; i++ ) {
  66. x = double(i-nw)*h;
  67. wxc[i] = x;
  68. wlc[i] = wavl[i];
  69. }
  70. }
  71. /*****************************************************************/
  72. double scalcoif4( double x, int m, int n, Vector& sc ) {
  73. const int ns = sc.dim()-1;
  74. int i;
  75. double xt, a, b, h, res;
  76. xt = pow( 2.0, double(m) )*x-double(n);
  77. if( xt >= -4.0 && xt <= 7.0 ) {
  78. h = 11.0/double(ns);
  79. i = int((xt+4.0)*double(ns)/11.0);
  80. if( i+1 == ns+1 ) a = -sc[i]/h;
  81. else a = (sc[i+1]-sc[i])/h;
  82. b = sc[i]-a*double(i)*h;
  83. res = a*(xt+4.0)+b;
  84. }
  85. else res = 0.0;
  86. res *= pow( sqrt(2.0), double(m) );
  87. return( res );
  88. }
  89. /*****************************************************************/
  90. double wavlcoif4( double x, int m, int n, Vector& wl ) {
  91. const int ns = wl.dim()-1;
  92. int i;
  93. double xt, a, b, h, res;
  94. xt = pow( 2.0, double(m) )*x-double(n);
  95. if( xt >= -5.0 && xt <= 6.0 ) {
  96. h = 11.0/double(ns);
  97. i = int((xt+5.0)*double(ns)/11.0);
  98. if( i+1 == ns+1 ) a = -wl[i]/h;
  99. else a = (wl[i+1]-wl[i])/h;
  100. b = wl[i]-a*double(i)*h;
  101. res = a*(xt+5.0)+b;
  102. }
  103. else res = 0.0;
  104. res *= pow( sqrt(2.0), double(m) );
  105. return( res );
  106. }
  107. /*****************************************************************/
  108. /****** Coifman scaling functions (L = 4) ******/
  109. /****** on the real line ******/
  110. /*****************************************************************/
  111. void coiflet_line( const int J, // scaling functions level
  112. Vector& qbase ) { // contains scaling functions
  113. // and wavelets
  114. Vector wvlt(qbase), scln(qbase);
  115. int k;
  116. int nit = 40; // number of iterations
  117. int nh = (qbase.dim() - 1) / 11; // number of subdivisions
  118. // of the unit
  119. coifman4( nh, nit, J, scln, wvlt );
  120. qbase = scln;
  121. }
  122. void coifman4( const int nh, const int NIT, const int J, Vector& sbase, Vector& wbase ) {
  123. double h[]= { 0.011587596739,-0.029320137980,-0.047639590310,
  124. 0.273021046535, 0.574682393857, 0.294867193696,
  125. -0.054085607092,-0.042026480461, 0.016744410163,
  126. 0.003967883613,-0.001289203356,-0.000509505399 };
  127. Vector hf(12,h);
  128. Matrix wvlt(11*nh+1,2);
  129. double sum;
  130. int i, it, n;
  131. int nit1, nit2, i0, i0p, i0m, nh4 = 4*nh, nh7 = 7*nh;
  132. /********* initial values for the iteration process *******/
  133. wvlt(nh4,0) = 1.0;
  134. double dnh = (double) nh;
  135. for( i = 1; i <= nh; i++ ) {
  136. wvlt((i + nh4),0) = wvlt((-i + nh4),0) = 1.0 - (double) i / dnh;
  137. }
  138. /********* iteration procedure for scaling function *******/
  139. for( it = 1; it <= NIT; it++ ) {
  140. nit2 = (int) fmod(double(it),2);
  141. nit1 = 1 - nit2;
  142. for( i = -nh4; i <= nh7; i++ ) {
  143. sum = 0.0;
  144. int i2 = 2 * i;
  145. for( n = -4; n <= 7; n++ ) {
  146. i0 = i2 - n*nh;
  147. if( ( i0 <= nh7 ) && ( i0 >= -nh4 ) )
  148. sum+= hf[n+4] * wvlt(nh4+i0,nit1) * 2.0;
  149. }
  150. wvlt(nh4+i,nit2) = sum;
  151. }
  152. }
  153. /********* set scaling function ***************************/
  154. sbase = wvlt.col(nit2) * pow( 2.0, double(J)/2.0 );
  155. /********* set wavelet ************************************/
  156. double vl = pow ( 2.0, double (J) / 2.0);
  157. for( i = -6*nh; i <= 5*nh; i++ ) {
  158. sum = 0.0;
  159. int i2 = 2 * i;
  160. for( n = -6; n <= 5; n++ ) {
  161. i0 = i2 - n*nh;
  162. if( ( i0 <= nh7 ) && ( i0 >= -nh4 ) )
  163. sum+= hf[1-n+4] * wvlt(nh4+i0,nit2) * 2.0 * pow( -1.0, double(n+1) );
  164. }
  165. wbase[6*nh+i] = sum * vl;
  166. }
  167. }
  168. /*****************************************************************/
  169. void coiflet_line_new( const int J,
  170. int& n0,
  171. Vector& qbase ) {
  172. Vector wvlt(qbase), scln(qbase);
  173. int k;
  174. int nit = 40;
  175. int nh = (qbase.dim() - 1) / 11;
  176. n0 = int( power( 2, J ) * nh );
  177. coifman4( nh, nit, J, scln, wvlt );
  178. qbase = scln;
  179. }
  180. /*****************************************************************/
  181. double scalcoif4_new( double x, int& n0, double& ax, double& bx, Vector& sc ) {
  182. int ix1, ix2;
  183. double xa, a, b, ba;
  184. double res = 0.0;
  185. xa = x - ax;
  186. ba = bx - ax;
  187. if( ( xa < ba ) && ( xa > 0.0 ) ) {
  188. ix1 = (int) ( xa * n0 );
  189. ix2 = ix1 + 1;
  190. if( ix2 == sc.dim() ) a = -sc[ix1] * n0;
  191. else a = (sc[ix2] - sc[ix1]) * n0;
  192. b = sc[ix1] - a * (double) ix1 / n0 ;
  193. res = a * xa + b;
  194. }
  195. return( res );
  196. }