wrefine_t.c

VC小波应用程序

/*
 *	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 refinement function
 *		FUNC_ND -- the name of the N-dimensional refinement function
 */
//所有用于细化处理的小波基的尺度函数类型
static void wrefine_step _PROTO((waveletfilter *wfltr, TYPE_ARRAY *a, int n));


//一维离散小波细化函数
/*
函数输入：
TYPE_ARRAY *a：		输入的原始数据，用于细化的小波基；
int nA：		数据数据a的大小，这个数据的大小必须是2的指数值；
int nNew：		细化后小波基的大小，同样必须是2的指数值；
waveletfilter *wfltr：	需要进行细化的小波滤波器；
TYPE_ARRAY *aRefined：	函数的实际输出，细化后的数据。
*/
void FUNC_1D(a, nA, nNew, wfltr, aRefined)
	TYPE_ARRAY *a;	
	int nA;				
	int nNew;		
	waveletfilter *wfltr;	
	TYPE_ARRAY *aRefined;
{
	int iA;
	double scale;

	assert(nA <= nNew);
	for (iA = 0; iA < nA; iA++)
		aRefined[iA] = a[iA];
	for (iA = 2 * nA; iA <= nNew; iA *= 2)
		wrefine_step(wfltr, aRefined, iA);

	scale = sqrt(RATIO(nNew, nA));
	for (iA = 0; iA < nNew; iA++)
		aRefined[iA] *= scale;

	return;
}

#if 0 

//完成N维小波基的细化工作
/*
函数输入：
TYPE_ARRAY *aO：	输入的原始数据，用于细化的小波基；
int nAO[]：		每一维上数据aO的大小，这个数据的大小必须是2的指数值；
int nDim：		小波基的维数以及数组nAO[]的大小；
bool isFwd：		小波基是否完整正变换（如果是，则取真值）
waveletfilter *wfltr：	需要进行细化的小波滤波器；
int nAR[]：		细化后小波基每一维的大小（数值必须是2的指数值）
TYPE_ARRAY *aR：	函数的实际输出，细化后的数据。
*/
void FUNC_ND(aO, nAO, nDim, isFwd, wfltr, nAR, aR)
	TYPE_ARRAY *aO;	
	int nAO[];	
	int nDim;	
	bool isFwd;		
	waveletfilter *wfltr;	
	int nAR[];	
	TYPE_ARRAY *aR;	
{
	int k;
	int iDim, nAOTotal, nAOMax, nARMax;
	int iO1, iO2, iO3, nO, nONext
	int nOPrev = 1;
	int nRPrev = 1;
	int nR, nRNext;
	TYPE_ARRAY *aOTmp, *aRTmp;

	nAOMax = 1;
	nARMax = 1;
	nAOTotal = 1;
	for (iDim = 0; iDim < nDim; iDim++) {
		if (nAO[iDim] > nAOMax)
			nAOMax = nAO[iDim];
		if (nAR[iDim] > nARMax)
			nARMax = nAR[iDim];
		nAOTotal *= nAO[iDim];
	}

	//临时数据空间必须足够大，这样能保证最大维数小波基的细化
	(void) MALLOC_LINTOK(aOTmp, nAOMax, TYPE_ARRAY);
	(void) MALLOC_LINTOK(aRTmp, nARMax, TYPE_ARRAY);

	for (iDim = 0; iDim < nDim; iDim++) {
		nO = nAO[nDim - 1 - iDim];
		nR = nAR[nDim - 1 - iDim];
		nONext = nO * nOPrev;
		nRNext = nR * nRPrev;
		if (nO >= MIN_ORDER) {
			for (iO2 = 0; iO2 < nAOTotal; iO2 += nONext) {
				for (iO1 = 0; iO1 < nOPrev; iO1++) {					
					 //将相关的行、列以及小波基其他信息进行复制
					iO3 = iO1 + iO2;
					for (k = 0; k < nO; k++) {
						aOTmp[k] = aO[iO3];
						iO3 += nOPrev;
					}
					
					//调用一维小波基的细化操作函数
					FUNC_1D(aOTmp, nO, nR, wfltr, aRTmp);

					//将细化的结果从林纾数据空间中粘贴至输出小波基中
					i3 = i1 + i2;
				for (k = 0; k < nR; k++) {
						aR[i3] = aRTmp[k];
						i3 += nRPrev;
					}
				}
			}
		}
		nOPrev = nONext;
		nRPrev = nRNext;
	}

	FREE_LINTOK(aRTmp);
	FREE_LINTOK(aOTmp);

	return;
}
#endif

/* wrefine_step -- perform one iteration of a wavelet refinement */
//完成小波细化操作中一次迭代循环操作
static void wrefine_step(wfltr, a, n)
	waveletfilter *wfltr;
	TYPE_ARRAY *a;
	int n;
{
	int nDiv2 = n / 2;
	int i, iHtilde, j, iA;
	TYPE_ARRAY *aNew;

	(void) MALLOC_LINTOK(aNew, n, TYPE_ARRAY);
	for (i = 0; i < n; i++)
		aNew[i] = 0.0;

	/*
	 *	Wavelet refinement is identical to an inverse wavelet transform (see
	 *	"wxfrm") with the elements of Gtilde set identically to zero.
	 *
	 *	A single step of the reconstruction is summarized by:
	 *
	 *		aNew[0..n-1] = Htilde^t * a[0..n/2-1]
	 *
	 *	where Htilde^t is the transpose of Htilde.
	 */
	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);
			aNew[iA] += wfltr->cHtilde[iHtilde] * a[j];
		}
	}
	for (i = 0; i < n; i++)
		a[i] = aNew[i];

	FREE_LINTOK(aNew);

	return;
}