00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006    
00007 /*
00008   This file contains routines for performing the Daubechies fast wavelet 
00009   transform analysis of time series data.
00010 
00011   File:    Daubechies.c
00012   Author:  B. Douglas Ward
00013   Date:    28 March 2000
00014 
00015 */
00016 
00017 
00018 /*---------------------------------------------------------------------------*/
00019 /*
00020   Calculate one iteration of the Daubechies forward FWT in 1-dimension.
00021 */
00022 
00023 void Daubechies_forward_pass_1d (int n, float * s)
00024 {
00025   int i;
00026   int npts;
00027   float * a = NULL;
00028   float * c = NULL;
00029   const float h[4] = { 0.683013, 1.18301, 0.316987, -0.183013 };
00030 
00031   npts = powerof2 (n);
00032   a = (float *) malloc (sizeof(float) * npts/2);
00033   c = (float *) malloc (sizeof(float) * npts/2);
00034   
00035   for (i = 0;  i < npts/2;  i++)
00036     {
00037       a[i] = (h[0]*s[(2*i)%npts] + h[1]*s[(2*i+1)%npts] + h[2]*s[(2*i+2)%npts]
00038               + h[3]*s[(2*i+3)%npts]) / 2.0;
00039       c[i] = (h[3]*s[(2*i)%npts] - h[2]*s[(2*i+1)%npts] + h[1]*s[(2*i+2)%npts] 
00040               - h[0]*s[(2*i+3)%npts]) / 2.0;
00041     }
00042 
00043   for (i = 0;  i < npts/2;  i++)
00044     {
00045       s[i] = a[i];
00046       s[i + npts/2] = c[i];
00047     }
00048 
00049   free (a);   a = NULL;
00050   free (c);   c = NULL;
00051 }
00052 
00053 
00054 /*---------------------------------------------------------------------------*/
00055 /*
00056   Calculate the Daubechies forward fast wavelet transform in 1-dimension.
00057 */
00058 
00059 void Daubechies_forward_FWT_1d (int n, float * s)
00060 {
00061   int m;
00062   int npts;
00063 
00064   npts = powerof2 (n);
00065 
00066   for (m = n-1;  m >= 0;  m--)
00067     {
00068       Daubechies_forward_pass_1d (m+1, s);
00069       /*
00070       ts_print (npts, s);
00071       */
00072     }
00073 }
00074 
00075 
00076 /*---------------------------------------------------------------------------*/
00077 /*
00078   Calculate one iteration of the Daubechies inverse FWT in 1-dimension.
00079 */
00080 
00081 void Daubechies_inverse_pass_1d (int n, float * s)
00082 {
00083   int i;
00084   int npts, nptsd2;
00085   float * a = NULL;
00086   float * c = NULL;
00087   float * r = NULL;
00088   const float h[4] = { 0.683013, 1.18301, 0.316987, -0.183013 };
00089 
00090 
00091   npts = powerof2 (n);
00092   nptsd2 = npts/2;
00093   a = s;
00094   c = s+nptsd2;
00095   r = (float *) malloc (sizeof(float) * npts);
00096   
00097 
00098   for (i = 0;  i < nptsd2;  i++)
00099     {
00100       r[2*i]   = h[2]*a[(i-1+nptsd2)%nptsd2] + h[1]*c[(i-1+nptsd2)%nptsd2] 
00101                + h[0]*a[i] + h[3]*c[i];
00102                
00103       r[2*i+1] = h[3]*a[(i-1+nptsd2)%nptsd2] - h[0]*c[(i-1+nptsd2)%nptsd2] 
00104                + h[1]*a[i] - h[2]*c[i];
00105     }
00106 
00107 
00108   for (i = 0;  i < npts;  i++)
00109     {
00110       s[i] = r[i];
00111     }
00112 
00113   free (r);   r = NULL;
00114 
00115 }
00116 
00117 
00118 /*---------------------------------------------------------------------------*/
00119 /*
00120   Calculate the Daubechies inverse fast wavelet transform in 1-dimension.
00121 */
00122 
00123 void Daubechies_inverse_FWT_1d (int n, float * s)
00124 {
00125   int m;
00126   int npts;
00127 
00128   npts = powerof2 (n);
00129 
00130   for (m = 1;  m <=n;  m++)
00131     {
00132       Daubechies_inverse_pass_1d (m, s);
00133       /*
00134       ts_print (npts, s);
00135       */
00136     }
00137 }
00138 
00139 
00140 /*---------------------------------------------------------------------------*/
00141 /*
00142   Calculate one iteration of the Daubechies forward FWT in 2-dimensions.
00143 */
00144 
00145 void Daubechies_forward_pass_2d (int n, float ** s)
00146 {
00147   int i, j;
00148   int npts;
00149   float * c = NULL;
00150 
00151 
00152   npts = powerof2 (n);
00153 
00154   for (i = 0;  i < npts;  i++)
00155     {
00156       Daubechies_forward_pass_1d (n, s[i]);
00157     }
00158 
00159   c = (float *) malloc (sizeof(float) * npts);
00160 
00161   for (j = 0;  j < npts;  j++)
00162     {
00163       for (i = 0;  i < npts;  i++)
00164         c[i] = s[i][j];
00165       Daubechies_forward_pass_1d (n, c);
00166       for (i = 0;  i < npts;  i++)
00167         s[i][j] = c[i];
00168     }
00169 
00170   free (c);   c = NULL;
00171 }
00172 
00173 
00174 /*---------------------------------------------------------------------------*/
00175 /*
00176   Calculate the Daubechies forward fast wavelet transform in 2-dimensions.
00177 */
00178 
00179 void Daubechies_forward_FWT_2d (int n, float ** s)
00180 {
00181   int m;
00182   int npts;
00183 
00184   npts = powerof2 (n);
00185 
00186   for (m = n-1;  m >= 0;  m--)
00187     {
00188       Daubechies_forward_pass_2d (m+1, s);
00189     }
00190 }
00191 
00192 
00193 /*---------------------------------------------------------------------------*/
00194 /*
00195   Calculate one iteration of the Daubechies inverse FWT in 2-dimensions.
00196 */
00197 
00198 void Daubechies_inverse_pass_2d (int n, float ** s)
00199 {
00200   int i, j;
00201   int npts;
00202   float * c = NULL;
00203 
00204 
00205   npts = powerof2 (n);
00206 
00207   for (i = 0;  i < npts;  i++)
00208     {
00209       Daubechies_inverse_pass_1d (n, s[i]);
00210     }
00211 
00212   c = (float *) malloc (sizeof(float) * npts);
00213 
00214   for (j = 0;  j < npts;  j++)
00215     {
00216       for (i = 0;  i < npts;  i++)
00217         c[i] = s[i][j];
00218       Daubechies_inverse_pass_1d (n, c);
00219       for (i = 0;  i < npts;  i++)
00220         s[i][j] = c[i];
00221     }
00222 
00223   free (c);   c = NULL;
00224 }
00225 
00226 
00227 /*---------------------------------------------------------------------------*/
00228 /*
00229   Calculate the Daubechies inverse fast wavelet transform in 2-dimensions.
00230 */
00231 
00232 void Daubechies_inverse_FWT_2d (int n, float ** s)
00233 {
00234   int m;
00235   int npts;
00236 
00237   npts = powerof2 (n);
00238 
00239   for (m = 1;  m <= n;  m++)
00240     {
00241       Daubechies_inverse_pass_2d (m, s);
00242     }
00243 }
00244 
00245 
00246 /*---------------------------------------------------------------------------*/
00247 
00248 
00249 
00250 
00251 
00252 
00253