WVLT.rar_visual c++ 小波

/* * This is a prototype file, not meant to be compiled directly. * The including source file must define the following: * * TYPE_ARRAY -- the base (scalar) type of all argument arrays * FUNC_1D -- the name of the 1-dimensional transform function * FUNC_ND -- the name of the N-dimensional transform function */ /* this is the minimum-order wavelet filter allowed */ #define MIN_ORDER 2 /* multidimensional transforms support a maximum of this many dimensions */ #define MXN_D 32 static void wfltr_convolve _PROTO((waveletfilter *wfltr, bool isFwd, TYPE_ARRAY *aIn, int incA, int n, TYPE_ARRAY *aXf)); static void wxfrm_1d_varstep _PROTO((TYPE_ARRAY *a, int incA, int nA, bool isFwd, waveletfilter *wfltr, TYPE_ARRAY *aXf)); static void wxfrm_nd_nonstd _PROTO((TYPE_ARRAY *a, int nA[], int nD, bool isFwd, waveletfilter *wfltr)); static void wxfrm_nd_std _PROTO((TYPE_ARRAY *a, int nA[], int nD, bool isFwd, waveletfilter *wfltr)); /* * To avoid lots of unnecessary malloc/free calls, especially in * multidimensional transforms, we use a "static global" buffer that, at the * highest level, will be allocated once to the maximum possible size and * freed when done. */ static TYPE_ARRAY *aTmp1D = NULL; /* wfltr_convolve -- perform one convolution of a (general) wavelet transform */ static void wfltr_convolve(wfltr, isFwd, aIn, incA, n, aXf) waveletfilter *wfltr; /* in: convolving filter */ bool isFwd; /* in: TRUE <=> forward transform */ TYPE_ARRAY *aIn; /* in: input data */ int incA; /* in: spacing of elements in aIn[] and aXf[] */ int n; /* in: size of aIn and aXf */ TYPE_ARRAY *aXf; /* out: output data (OK if == aIn) */ { int nDiv2 = n / 2; int iA, iHtilde, iGtilde, j, jH, jG, i; double sum; bool flip; /* * According to Daubechies: * * H is the analyzing smoothing filter. * G is the analyzing detail filter. * Htilde is the reconstruction smoothing filter. * Gtilde is the reconstruction detail filter. */ /* the reconstruction detail filter is the mirror of the analysis smoothing filter */ int nGtilde = wfltr->nH; /* the analysis detail filter is the mirror of the reconstruction smoothing filter */ int nG = wfltr->nHtilde; if (isFwd) { /* * A single step of the analysis is summarized by: * aTmp1D[0..n/2-1] = H * a (smooth component) * aTmp1D[n/2..n-1] = G * a (detail component) */ for (i = 0; i < nDiv2; i++) { /* * aTmp1D[0..nDiv2-1] contains the smooth components. */ sum = 0.0; for (jH = 0; jH < wfltr->nH; jH++) { /* * Each row of H is offset by 2 from the previous one. * * We assume our data is periodic, so we wrap the aIn[] values * if necessary. If we have more samples than we have H * coefficients, we also wrap the H coefficients. */ iA = MOD(2 * i + jH - wfltr->offH, n); sum += wfltr->cH[jH] * aIn[incA * iA]; } aTmp1D[i] = sum; /* * aTmp1D[nDiv2..n-1] contains the detail components. */ sum = 0.0; flip = TRUE; for (jG = 0; jG < nG; jG++) { /* * We construct the G coefficients on-the-fly from the * Htilde coefficients. * * Like H, each row of G is offset by 2 from the previous * one. As with H, we also allow the coefficents of G to * wrap. * * Again as before, the aIn[] values may wrap. */ iA = MOD(2 * i + jG - wfltr->offG, n); if (flip) sum -= wfltr->cHtilde[nG - 1 - jG] * aIn[incA * iA]; else sum += wfltr->cHtilde[nG - 1 - jG] * aIn[incA * iA]; flip = !flip; } aTmp1D[nDiv2 + i] = sum; } } else { /* * The inverse transform is a little trickier to do efficiently. * A single step of the reconstruction is summarized by: * * aTmp1D = Htilde^t * aIn[incA * (0..n/2-1)] * + Gtilde^t * aIn[incA * (n/2..n-1)] * * where x^t is the transpose of x. */ for (i = 0; i < n; i++) aTmp1D[i] = 0.0; /* necessary */ for (j = 0; j < nDiv2; j++) { for (iHtilde = 0; iHtilde < wfltr->nHtilde; iHtilde++) { /* * Each row of Htilde is offset by 2 from the previous one. */ iA = MOD(2 * j + iHtilde - wfltr->offHtilde, n); aTmp1D[iA] += wfltr->cHtilde[iHtilde] * aIn[incA * j]; } flip = TRUE; for (iGtilde = 0; iGtilde < nGtilde; iGtilde++) { /* * As with Htilde, we also allow the coefficents of Gtilde, * which is the mirror of H, to wrap. * * We assume our data is periodic, so we wrap the aIn[] values * if necessary. If we have more samples than we have Gtilde * coefficients, we also wrap the Gtilde coefficients. */ iA = MOD(2 * j + iGtilde - wfltr->offGtilde, n); if (flip) { aTmp1D[iA] -= wfltr->cH[nGtilde - 1 - iGtilde] * aIn[incA * (j + nDiv2)]; } else { aTmp1D[iA] += wfltr->cH[nGtilde - 1 - iGtilde] * aIn[incA * (j + nDiv2)]; } flip = !flip; } } } for (i = 0; i < n; i++) aXf[incA * i] = aTmp1D[i]; return; } /* wxfrm_1d_varstep -- 1-dimensional discrete wavelet transform with variable step size */ static void wxfrm_1d_varstep(aIn, incA, nA, isFwd, wfltr, aXf) TYPE_ARRAY *aIn; /* in: original array */ int incA; /* in: spacing of elements in aIn[] and aXf[] */ int nA; /* in: size of a (must be power of 2) */ bool isFwd; /* in: TRUE <=> forward transform */ waveletfilter *wfltr; /* in: wavelet filter to use */ TYPE_ARRAY *aXf; /* out: transformed array */ { int iA; if (nA < MIN_ORDER) return; if (isFwd) { /* * Start at largest hierarchy, work towards smallest. */ wfltr_convolve(wfltr, isFwd, aIn, incA, nA, aXf); for (iA = nA / 2; iA >= MIN_ORDER; iA /= 2) wfltr_convolve(wfltr, isFwd, aXf, incA, iA, aXf); } else { /* * Start at smallest hierarchy, work towards largest. */ if (aXf != aIn) { for (iA = 0; iA < nA; iA++) aXf[incA * iA] = aIn[incA * iA]; /* required for inverse */ } for (iA = MIN_ORDER; iA <= nA; iA *= 2) wfltr_convolve(wfltr, isFwd, aXf, incA, iA, aXf); } return; } /* wxfrm_[fd]a1d -- 1-dimensional discrete wavelet transform */ void FUNC_1D(a, nA, isFwd, wfltr, aXf) TYPE_ARRAY *a; /* in: original array */ int nA; /* in: size of a (must be power of 2) */ bool isFwd; /* in: TRUE <=> forward transform */ waveletfilter *wfltr; /* in: wavelet filter to use */ TYPE_ARRAY *aXf; /* out: transformed array */ { assert(aTmp1D == NULL); (void) MALLOC_LINTOK(aTmp1D, nA, TYPE_ARRAY); wxfrm_1d_varstep(a, 1, nA, isFwd, wfltr, aXf); FREE_LINTOK(aTmp1D); aTmp1D = NULL; return; } /* wxfrm_[fd]and -- n-dimensional discrete wavelet transform */ void FUNC_ND(aIn, nA, nD, isFwd, isStd, wfltr, aXf) TYPE_ARRAY *aIn; /* in: original data */ int nA[]; /* in: size of each dimension of a (must be powers of 2) */ int nD; /* in: number of dimensions (size of nA[]) */ bool isFwd; /* in: TRUE <=> forward transform */ bool isStd; /* in: TRUE <=> standard basis */ waveletfilter *wfltr; /* in: wavelet filter to use */ TYPE_ARRAY *aXf; /* out: transformed data (ok if == a) */ { int k, d, nDMax, nATot; int lg2nA; if (nD < 1 || MXN_D < nD) { (void) fprintf(stderr, "number of dimensions must be between 1 and %d -- exiting\n", MXN_D); exit(1); } for (d = 0; d < nD; d++) { for (lg2nA = 0; (1 << lg2nA) < nA[d]; lg2nA++) continue; if ((1 << lg2nA) != nA[d]) { (void) fprintf(stderr, "size of dimension #%d (= %d) is not a power of 2 -- exiting\n", d, nA[d]); exit(1); } } nDMax = 1; nATot = 1; for (d = 0; d < nD; d++) { if (nA[d] > nDMax) nDMax = nA[d]; nATot *= nA[d]; } /* * Make the "static global" temp array aTmp1D, big enough to hold the * longest dimension. */ assert(aTmp1D == NULL); (void) MALLOC_