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 Haar fast wavelet transform
00009   analysis of time series data.
00010 
00011   File:    Haar.c
00012   Author:  B. Douglas Ward
00013   Date:    28 March 2000
00014 
00015 */
00016 
00017 /*---------------------------------------------------------------------------*/
00018 /*
00019   Calculate Haar in-place forward fast wavelet transform in 1-dimension.
00020 */
00021 
00022 void Haar_ip_FFWT_1d (int n, float * s)
00023 {
00024   float a;
00025   float c;
00026   int i;
00027   int j;
00028   int k;
00029   int l;
00030   int m;
00031 
00032 
00033   i = 1;
00034   j = 2;
00035   m = powerof2 (n);
00036 
00037   for (l = n-1;  l >= 0;  l--)
00038     {
00039       printf ("l = %d \n", l);
00040       m /= 2;
00041 
00042       for (k = 0;  k < m;  k++)
00043         {
00044           a = (s[j*k] + s[j*k+i]) / 2.0;
00045           c = (s[j*k] - s[j*k+i]) / 2.0;
00046           s[j*k] = a;
00047           s[j*k+i] = c;
00048         }
00049       
00050       i *= 2;
00051       j *= 2;
00052       /*
00053       ts_print (powerof2(n), s);
00054       */
00055     }
00056 }
00057 
00058 
00059 /*---------------------------------------------------------------------------*/
00060 /*
00061   Calculate Haar in-place inverse fast wavelet transform in 1-dimension.
00062 */
00063 
00064 void Haar_ip_IFWT_1d (int n, float * s)
00065 {
00066   float a0;
00067   float a1;
00068   int i;
00069   int j;
00070   int k;
00071   int l;
00072   int m;
00073 
00074 
00075   i = powerof2 (n-1);
00076   j = 2*i;
00077   m = 1;
00078 
00079   for (l = 1;  l <= n;  l++)
00080     {
00081       printf ("l = %d \n", l);
00082 
00083       for (k = 0;  k < m;  k++)
00084         {
00085           a0 = s[j*k] + s[j*k+i];
00086           a1 = s[j*k] - s[j*k+i];
00087           s[j*k] = a0;
00088           s[j*k+i] = a1;
00089         }
00090       
00091       i /= 2;
00092       j /= 2;
00093       m *= 2;
00094       /*
00095       ts_print (powerof2(n), s);
00096       */
00097     }
00098 }
00099 
00100 
00101 /*---------------------------------------------------------------------------*/
00102 /*
00103   Calculate one iteration of the Haar forward FWT in 1-dimension.
00104 */
00105 
00106 void Haar_forward_pass_1d (int n, float * s)
00107 {
00108   int i;
00109   int npts;
00110   float * a = NULL;
00111   float * c = NULL;
00112 
00113   npts = powerof2 (n);
00114   a = (float *) malloc (sizeof(float) * npts/2);
00115   c = (float *) malloc (sizeof(float) * npts/2);
00116   
00117   for (i = 0;  i < npts/2;  i++)
00118     {
00119       a[i] = (s[2*i] + s[2*i+1]) / 2.0;
00120       c[i] = (s[2*i] - s[2*i+1]) / 2.0;
00121     }
00122 
00123   for (i = 0;  i < npts/2;  i++)
00124     {
00125       s[i] = a[i];
00126       s[i + npts/2] = c[i];
00127     }
00128 
00129   free (a);   a = NULL;
00130   free (c);   c = NULL;
00131 }
00132 
00133 
00134 /*---------------------------------------------------------------------------*/
00135 /*
00136   Calculate the Haar forward fast wavelet transform in 1-dimension.
00137 */
00138 
00139 void Haar_forward_FWT_1d (int n, float * s)
00140 {
00141   int m;
00142   int npts;
00143 
00144   npts = powerof2 (n);
00145 
00146   for (m = n-1;  m >= 0;  m--)
00147     {
00148       Haar_forward_pass_1d (m+1, s);
00149       /*
00150       ts_print (npts, s);
00151       */
00152     }
00153 }
00154 
00155 
00156 /*---------------------------------------------------------------------------*/
00157 /*
00158   Calculate one iteration of the Haar inverse FWT in 1-dimension.
00159 */
00160 
00161 void Haar_inverse_pass_1d (int n, float * s)
00162 {
00163   int i;
00164   int npts;
00165   float * r = NULL;
00166 
00167   npts = powerof2 (n);
00168   r = (float *) malloc (sizeof(float) * npts);
00169   
00170   for (i = 0;  i < npts/2;  i++)
00171     {
00172       r[2*i]   = s[i] + s[i + npts/2];
00173       r[2*i+1] = s[i] - s[i + npts/2]; 
00174     }
00175 
00176   for (i = 0;  i < npts;  i++)
00177     {
00178       s[i] = r[i];
00179     }
00180 
00181   free (r);   r = NULL;
00182 }
00183 
00184 
00185 /*---------------------------------------------------------------------------*/
00186 /*
00187   Calculate the Haar inverse fast wavelet transform in 1-dimension.
00188 */
00189 
00190 void Haar_inverse_FWT_1d (int n, float * s)
00191 {
00192   int m;
00193   int npts;
00194 
00195   npts = powerof2 (n);
00196 
00197   for (m = 1;  m <=n;  m++)
00198     {
00199       Haar_inverse_pass_1d (m, s);
00200       /*
00201       ts_print (npts, s);
00202       */
00203     }
00204 }
00205 
00206 
00207 /*---------------------------------------------------------------------------*/
00208 /*
00209   Calculate one iteration of the Haar forward FWT in 2-dimensions.
00210 */
00211 
00212 void Haar_forward_pass_2d (int n, float ** s)
00213 {
00214   int i, j;
00215   int npts;
00216   float * c = NULL;
00217 
00218 
00219   npts = powerof2 (n);
00220 
00221   for (i = 0;  i < npts;  i++)
00222     {
00223       Haar_forward_pass_1d (n, s[i]);
00224     }
00225 
00226   c = (float *) malloc (sizeof(float) * npts);
00227 
00228   for (j = 0;  j < npts;  j++)
00229     {
00230       for (i = 0;  i < npts;  i++)
00231         c[i] = s[i][j];
00232       Haar_forward_pass_1d (n, c);
00233       for (i = 0;  i < npts;  i++)
00234         s[i][j] = c[i];
00235     }
00236 
00237   free (c);   c = NULL;
00238 }
00239 
00240 
00241 /*---------------------------------------------------------------------------*/
00242 /*
00243   Calculate the Haar forward fast wavelet transform in 2-dimensions.
00244 */
00245 
00246 void Haar_forward_FWT_2d (int n, float ** s)
00247 {
00248   int m;
00249   int npts;
00250 
00251   npts = powerof2 (n);
00252 
00253   for (m = n-1;  m >= 0;  m--)
00254     {
00255       Haar_forward_pass_2d (m+1, s);
00256     }
00257 }
00258 
00259 
00260 /*---------------------------------------------------------------------------*/
00261 /*
00262   Calculate one iteration of the Haar inverse FWT in 2-dimensions.
00263 */
00264 
00265 void Haar_inverse_pass_2d (int n, float ** s)
00266 {
00267   int i, j;
00268   int npts;
00269   float * c = NULL;
00270 
00271 
00272   npts = powerof2 (n);
00273 
00274   for (i = 0;  i < npts;  i++)
00275     {
00276       Haar_inverse_pass_1d (n, s[i]);
00277     }
00278 
00279   c = (float *) malloc (sizeof(float) * npts);
00280 
00281   for (j = 0;  j < npts;  j++)
00282     {
00283       for (i = 0;  i < npts;  i++)
00284         c[i] = s[i][j];
00285       Haar_inverse_pass_1d (n, c);
00286       for (i = 0;  i < npts;  i++)
00287         s[i][j] = c[i];
00288     }
00289 
00290   free (c);   c = NULL;
00291 }
00292 
00293 
00294 /*---------------------------------------------------------------------------*/
00295 /*
00296   Calculate the Haar inverse fast wavelet transform in 2-dimensions.
00297 */
00298 
00299 void Haar_inverse_FWT_2d (int n, float ** s)
00300 {
00301   int m;
00302   int npts;
00303 
00304   npts = powerof2 (n);
00305 
00306   for (m = 1;  m <= n;  m++)
00307     {
00308       Haar_inverse_pass_2d (m, s);
00309     }
00310 }
00311 
00312 
00313 /*---------------------------------------------------------------------------*/
00314 
00315 
00316