1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 /*
8 This file contains routines for performing the Daubechies fast wavelet
9 transform analysis of time series data.
10
11 File: Daubechies.c
12 Author: B. Douglas Ward
13 Date: 28 March 2000
14
15 */
16
17
18 /*---------------------------------------------------------------------------*/
19 /*
20 Calculate one iteration of the Daubechies forward FWT in 1-dimension.
21 */
22
Daubechies_forward_pass_1d(int n,float * s)23 void Daubechies_forward_pass_1d (int n, float * s)
24 {
25 int i;
26 int npts;
27 float * a = NULL;
28 float * c = NULL;
29 const float h[4] = { 0.683013, 1.18301, 0.316987, -0.183013 };
30
31 npts = powerof2 (n);
32 a = (float *) malloc (sizeof(float) * npts/2);
33 c = (float *) malloc (sizeof(float) * npts/2);
34
35 for (i = 0; i < npts/2; i++)
36 {
37 a[i] = (h[0]*s[(2*i)%npts] + h[1]*s[(2*i+1)%npts] + h[2]*s[(2*i+2)%npts]
38 + h[3]*s[(2*i+3)%npts]) / 2.0;
39 c[i] = (h[3]*s[(2*i)%npts] - h[2]*s[(2*i+1)%npts] + h[1]*s[(2*i+2)%npts]
40 - h[0]*s[(2*i+3)%npts]) / 2.0;
41 }
42
43 for (i = 0; i < npts/2; i++)
44 {
45 s[i] = a[i];
46 s[i + npts/2] = c[i];
47 }
48
49 free (a); a = NULL;
50 free (c); c = NULL;
51 }
52
53
54 /*---------------------------------------------------------------------------*/
55 /*
56 Calculate the Daubechies forward fast wavelet transform in 1-dimension.
57 */
58
Daubechies_forward_FWT_1d(int n,float * s)59 void Daubechies_forward_FWT_1d (int n, float * s)
60 {
61 int m;
62 int npts;
63
64 npts = powerof2 (n);
65
66 for (m = n-1; m >= 0; m--)
67 {
68 Daubechies_forward_pass_1d (m+1, s);
69 /*
70 ts_print (npts, s);
71 */
72 }
73 }
74
75
76 /*---------------------------------------------------------------------------*/
77 /*
78 Calculate one iteration of the Daubechies inverse FWT in 1-dimension.
79 */
80
Daubechies_inverse_pass_1d(int n,float * s)81 void Daubechies_inverse_pass_1d (int n, float * s)
82 {
83 int i;
84 int npts, nptsd2;
85 float * a = NULL;
86 float * c = NULL;
87 float * r = NULL;
88 const float h[4] = { 0.683013, 1.18301, 0.316987, -0.183013 };
89
90
91 npts = powerof2 (n);
92 nptsd2 = npts/2;
93 a = s;
94 c = s+nptsd2;
95 r = (float *) malloc (sizeof(float) * npts);
96
97
98 for (i = 0; i < nptsd2; i++)
99 {
100 r[2*i] = h[2]*a[(i-1+nptsd2)%nptsd2] + h[1]*c[(i-1+nptsd2)%nptsd2]
101 + h[0]*a[i] + h[3]*c[i];
102
103 r[2*i+1] = h[3]*a[(i-1+nptsd2)%nptsd2] - h[0]*c[(i-1+nptsd2)%nptsd2]
104 + h[1]*a[i] - h[2]*c[i];
105 }
106
107
108 for (i = 0; i < npts; i++)
109 {
110 s[i] = r[i];
111 }
112
113 free (r); r = NULL;
114
115 }
116
117
118 /*---------------------------------------------------------------------------*/
119 /*
120 Calculate the Daubechies inverse fast wavelet transform in 1-dimension.
121 */
122
Daubechies_inverse_FWT_1d(int n,float * s)123 void Daubechies_inverse_FWT_1d (int n, float * s)
124 {
125 int m;
126 int npts;
127
128 npts = powerof2 (n);
129
130 for (m = 1; m <=n; m++)
131 {
132 Daubechies_inverse_pass_1d (m, s);
133 /*
134 ts_print (npts, s);
135 */
136 }
137 }
138
139
140 /*---------------------------------------------------------------------------*/
141 /*
142 Calculate one iteration of the Daubechies forward FWT in 2-dimensions.
143 */
144
Daubechies_forward_pass_2d(int n,float ** s)145 void Daubechies_forward_pass_2d (int n, float ** s)
146 {
147 int i, j;
148 int npts;
149 float * c = NULL;
150
151
152 npts = powerof2 (n);
153
154 for (i = 0; i < npts; i++)
155 {
156 Daubechies_forward_pass_1d (n, s[i]);
157 }
158
159 c = (float *) malloc (sizeof(float) * npts);
160
161 for (j = 0; j < npts; j++)
162 {
163 for (i = 0; i < npts; i++)
164 c[i] = s[i][j];
165 Daubechies_forward_pass_1d (n, c);
166 for (i = 0; i < npts; i++)
167 s[i][j] = c[i];
168 }
169
170 free (c); c = NULL;
171 }
172
173
174 /*---------------------------------------------------------------------------*/
175 /*
176 Calculate the Daubechies forward fast wavelet transform in 2-dimensions.
177 */
178
Daubechies_forward_FWT_2d(int n,float ** s)179 void Daubechies_forward_FWT_2d (int n, float ** s)
180 {
181 int m;
182 int npts;
183
184 npts = powerof2 (n);
185
186 for (m = n-1; m >= 0; m--)
187 {
188 Daubechies_forward_pass_2d (m+1, s);
189 }
190 }
191
192
193 /*---------------------------------------------------------------------------*/
194 /*
195 Calculate one iteration of the Daubechies inverse FWT in 2-dimensions.
196 */
197
Daubechies_inverse_pass_2d(int n,float ** s)198 void Daubechies_inverse_pass_2d (int n, float ** s)
199 {
200 int i, j;
201 int npts;
202 float * c = NULL;
203
204
205 npts = powerof2 (n);
206
207 for (i = 0; i < npts; i++)
208 {
209 Daubechies_inverse_pass_1d (n, s[i]);
210 }
211
212 c = (float *) malloc (sizeof(float) * npts);
213
214 for (j = 0; j < npts; j++)
215 {
216 for (i = 0; i < npts; i++)
217 c[i] = s[i][j];
218 Daubechies_inverse_pass_1d (n, c);
219 for (i = 0; i < npts; i++)
220 s[i][j] = c[i];
221 }
222
223 free (c); c = NULL;
224 }
225
226
227 /*---------------------------------------------------------------------------*/
228 /*
229 Calculate the Daubechies inverse fast wavelet transform in 2-dimensions.
230 */
231
Daubechies_inverse_FWT_2d(int n,float ** s)232 void Daubechies_inverse_FWT_2d (int n, float ** s)
233 {
234 int m;
235 int npts;
236
237 npts = powerof2 (n);
238
239 for (m = 1; m <= n; m++)
240 {
241 Daubechies_inverse_pass_2d (m, s);
242 }
243 }
244
245
246 /*---------------------------------------------------------------------------*/
247
248
249
250
251
252
253
254