/*------------------------------------------------------------------------------*
 * File Name: wavelet_utils.h												    *
 * Creation: Sandy 4/24/06														*
 * Purpose: support utility functions for wavelet functions               		*
 * Copyright (c) OriginLab Corp.2005											*
 * All Rights Reserved															*
 * 																				*
 * Modification Log:															*
 *	Folger 07/28/08 BETTER_OPREATION_OUTPUT_COLUMN_COMMENTS						*
 *------------------------------------------------------------------------------*/

#ifndef _WAVELET_UTILS_H
#define _WAVELET_UTILS_H

#include  <nag8\OC_nag_ex.h>

static int g_nWaveletType[] = {
	Nag_Haar,
	Nag_DB2, Nag_DB3, Nag_DB4, Nag_DB5, Nag_DB6, Nag_DB7, Nag_DB8, Nag_DB9, Nag_DB10,
	Nag_Bior1_1, Nag_Bior1_3, Nag_Bior1_5, Nag_Bior2_2, Nag_Bior2_4, Nag_Bior2_6, Nag_Bior2_8, Nag_Bior3_1, Nag_Bior3_3, Nag_Bior3_5, Nag_Bior3_7
};

static int g_nCwtWaveletType[] = {
	Nag_Morlet, Nag_DGauss, Nag_MexHat
};



static int g_nWaveModes[] = {Nag_PerWavExt, Nag_ZeroWavExt};


//Category last updated by Cloud on 2008-6-6
//Commetnts last updated by Sandy on 2006-4-25
/** >Analysis>Signal Processing
	Remarks:
		Compute the full length decomposing a signal with definite level. 
	Paramaters:
		nSignalLen = [input] Length of the signal.
		nLevel     = [input] the level want to decompose

	Example1:
		void get_enough_length_to_be_decompose_ex1()
		{
		
			vector vSignal;
			vSignal.Data(0,129,1);
			
			int nSignalLen=vSignal.GetSize();
			int nLen = get_enough_length_to_be_decompose(nSignalLen, 7);
			//nLen = 256;
		}
	Return:
		Return the enough length of signal to decompose with multi-levels.
*/
uint get_enough_length_to_be_decompose(int nSignalLen, int nLevel);


//Category last updated by Cloud on 2008-6-6
//Commetnts last updated by Sandy on 2006-4-25
/** >Analysis>Signal Processing
	Remarks:
		Compute the reconstructed signal's size ahead of implement the invert transform. 
	Paramaters:
		mode 		   = [input] The Extend mode.
		nSizeofCoeffs  = [input] size of coeffs
		filter		   = [input] wavelet filter

	Example1:Put the two coeffs into 1st and 2nd col of an actived worksheet is needed.
	         
		void get_reconstruct_size_ex1()
		{
			Worksheet wks = Project.ActiveLayer();
			if(!wks)
				return;
			
			Dataset ds1(wks, 0), ds2(wks, 1);
			vector vca, vcd;
			vca = ds1, vcd = ds2;
			Nag_WavFilt filter;

			int i, mshift, n, nc, ny;
			
			if(!create_wavelet_filter(g_nWaveletType[0], filter))
				return;
			
			Nag_WavModeType mode = nWaveModes[0];
			
			
			nc = vca.GetSize();
			
			ny = get_reconstruct_size(mode, nc, filter);
			
			
			vector vy(ny);
			
			NagError fail;
			fail.code = NE_NOERROR;
			nag_idwt(mode, nc, vca, vcd, &filter, ny, vy, &fail);
			if (fail.code != NE_NOERROR)
				return;			
			
			nag_free_wavfilter(&filter);
			
			int nCol = wks.AddCol();
			Dataset ds3(wks, nCol);
			ds3.SetSize(vy.GetSize());
			ds3 = vy;
		}
	Return:
		Return the size of reconstructed signal
*/
uint get_reconstruct_size(Nag_WavModeType mode, int nSizeofCoeffs, Nag_WavFilt& filter);


//Category last updated by Cloud on 2008-6-6
//Commetnts last updated by Sandy on 2006-4-25
/** >Analysis>Signal Processing
	Remarks:
		create a wavelet filter with the input wavelet name. 
	Paramaters:
		wavname		   = [input] The name of wavelet.
		filter		   = [output] wavelet filter

	Example1:        
		void create_wavelet_filter_ex1()
		{
			Nag_WavFilt filter;
			NagError fail;
			//INIT_FAIL
			fail.code = NE_NOERROR;
			fail.print = 0;
			
			if(!create_wavelet_filter(g_nWaveletType[0], filter))
			 	return;
			
			//using filter
			
			nag_free_wavfilter(&filter);
		}
		
	Return:
		Return true if success
*/
bool create_wavelet_filter(Nag_WavType wavname, Nag_WavFilt &filter);

//Category last updated by Cloud on 2008-6-6
//Commetnts last updated by Sandy on 2006-4-25
/**>Analysis>Signal Processing
	Remarks:
		Compute the coefficients's size ahead of implement the dwt transform. 
	Paramaters:
		mode 		   = [input] The Extend mode.
		nSizeofCoeffs  = [input] size of signal
		filter		   = [input] wavelet filter

	Example1:Put the signal into 1st col of an actived worksheet is needed.
	         
		void get_coeffs_size_ex1()
		{
			Worksheet wks = Project.ActiveLayer();
			if(!wks)
				return;
			
			Dataset ds1(wks, 0);
			vector vs;
			vs = ds1;
			
			Nag_WavFilt filter;

			int i, n, nc;
			
			if(!create_wavelet_filter(g_nWaveletType[0], filter))
				return;

			Nag_WavModeType mode = nWaveModes[0];
			n = vs.GetSize();
			nc = get_coeffs_size(mode, n, filter);
			vs.SetSize(nc*2);
			
			vector ca(nc), cd(nc);
						
			NagError fail;
			fail.code = NE_NOERROR;
			nag_dwt(mode, vs.GetSize(), vs, &filter, nc, ca, cd, &fail);
			if (fail.code != NE_NOERROR)
			{
				nag_free_wavfilter(&filter);
				return;
			}
			
			nag_free_wavfilter(&filter);
			
			int nCol1 = wks.AddCol();
			int nCol2 = wks.AddCol();
			Dataset ds2(wks, nCol1);
			Dataset ds3(wks, nCol2);
			ds2.SetSize(ca.GetSize());
			ds2 = ca;
			ds3.SetSize(cd.GetSize());
			ds3 = cd;
		}
		
	Return:
		Return the size of coefficients
*/
uint get_coeffs_size(Nag_WavModeType mode, int nSizeofSingal, Nag_WavFilt& filter);


//Category last updated by Cloud on 2008-6-6
//Commetnts last updated by Sandy on 2006-4-25
/**>Analysis>Signal Processing
	Remarks:
		Multilevel 1-D wavelet decomposition. 
	Paramaters:
		nLevel         = [input] The level to decompose signal at 
		mode		   = [input] The Extend mode.
		vSignal        = [input] The signal vector.
		filter		   = [input] Wavelet filter.
		vCoeff		   = [output] output all coefficients in a decomposition vector, each coefficient's length refer to vL.
		vL			   = [output] bookkeeping vector of vCoeff.

	Example1:  Put the signal into 1st col of an actived worksheet is needed.      
		void multi_scale_dwt_ex1()
		{
			Worksheet wks = Project.ActiveLayer();
			if(!wks)
				return;
			
			Dataset ds1(wks, 0);
			vector vs;
			vs = ds1;
			
			Nag_WavFilt filter;
			
			Nag_WavModeType mode = nWaveModes[0];

			int i, n, nc, nLevel;
			
			int nSignalLen=vs.GetSize();
			
			nLevel=min(7,(int)(log10(nSignalLen)/log10(2)));
			
			n = get_enough_length_to_be_decompose(nSignalLen, nLevel);
			
			vs.SetSize(n);
			
			if(!create_wavelet_filter(g_nWaveletType[0], filter))
				return;
			
			vector vC,vL;	
			if(!multi_scale_dwt(nLevel,mode,vs,filter,vC,vL))
			{
				nag_free_wavfilter(&filter);
				return;

			}
			nag_free_wavfilter(&filter);
			
			int nCol1 = wks.AddCol();
			int nCol2 = wks.AddCol();
			Dataset ds2(wks, nCol1);
			Dataset ds3(wks, nCol2);
			ds2.SetSize(vC.GetSize());
			ds2 = vC;
			ds3.SetSize(vL.GetSize());
			ds3 = vL;
			
		}
		
	Return:
		Return true if success
*/
bool multi_scale_dwt(int nLevel, Nag_WavModeType mode, vector& vSignal, Nag_WavFilt filter, vector& vCoeff, vector& vL);


//Category last updated by Cloud on 2008-6-6
//Commetnts last updated by Sandy on 2006-4-25
/**>Analysis>Signal Processing
	Remarks:
		Multilevel 1-D wavelet reconstruct
	Paramaters:
		nLevel         = [input] The level to reconstruct signal 
		mode		   = [input] The Extend mode.
		vSignal        = [output] The recontrusted signal vector.
		filter		   = [input] Wavelet filter.
		vCoeff		   = [input] a decomposition vector, each coefficient's length refer to vL.
		vL			   = [input] bookkeeping vector of vCoeff.

	Example1:   Put the signal into 1st col of an actived worksheet is needed. 
	            input threshold value of which level(only this level) need to cut.     
		void multi_scale_idwt_ex1(int level, double threshold)
		{
			Worksheet wks = Project.ActiveLayer();
			if(!wks)
				return;
			
			Dataset ds1(wks, 0);
			vector vx;
			vx = ds1;
			
			Nag_WavFilt filter;

			Nag_WavModeType mode = nWaveModes[0];
			
			int i, n, nc, nLevel;
			int nSignalLen=vx.GetSize();
			
			nLevel=min(7,(int)(log10(nSignalLen)/log10(2)));
			
			n = get_enough_length_to_be_decompose(nSignalLen, nLevel);
			
			vx.SetSize(n);
			

			if(!create_wavelet_filter(g_nWaveletType[0], filter))
				return;
				
			vector vC,vL;	
			if(!multi_scale_dwt(nLevel,mode,vx,filter,vC,vL))
			{
				nag_free_wavfilter(&filter);
				return;
			}
			
			//get detail coefficients of certain level
			int nWholeLen=vL.GetSize();
			int len=vL[nWholeLen-level];
			vector cd;
			cd.SetSize(len);
			int nl=vL[0];
			for(i=0;i<(nLevel-level);i++)
			{
				nl+=vL[i+1];
			}
			vC.GetSubVector(cd, nl, nl+len-1);
			
			//discard the small data of cd with theshold 
			cd.Replace(threshold,0,MATREPL_USE_ABSOLUTE_VALUE_IN_TEST|MATREPL_TEST_LESSTHAN);
			vC.SetSubVector(cd, nl);
					
			vector vy;
			vy.SetSize(vx.GetSize());
			if(!multi_scale_idwt(nLevel,mode,vy,filter,vC,vL))
			{
				nag_free_wavfilter(&filter);
				return;
			}
			
			vy.SetSize(nSignalLen);

			nag_free_wavfilter(&filter);
			
			int nCol = wks.AddCol();
			Dataset ds3(wks, nCol);

			ds3.SetSize(vy.GetSize());
			ds3 = vy;
		}
		
	Return:
		Return true if success
*/
bool multi_scale_idwt(int nLevel, Nag_WavModeType mode, vector& vSignal, Nag_WavFilt filter, vector& vCoeff, vector& vL);

//------ Folger 07/28/08 BETTER_OPREATION_OUTPUT_COLUMN_COMMENTS
string	get_wavelet_type_string(int nType);
//------