T536_DTU/app_public/wavelib/wavelib.c

/*******************************************************************************
版权所有：	珠海欧力配网自动化股份有限公司
文件名称：	wavelib.c
文件版本：	01.00
创建作者:	ygl
创建日期：	2022-03-09
功能说明：	modwt小波计算
其它说明：	支持：DB,sym、coif、haar
注意事项:只求近似值，不考虑数据重构
*/
#ifdef  __KERNEL__

#include "head.h"
#define  M_SQRT2     1.41421356237309504880

#else
#include "head.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define _USE_MATH_DEFINES
#include <math.h>
#include "wavelib.h"

#define rt_malloc    malloc
// #define rt_printf    printf
#define rt_free      free

#endif

hilbert_para coef1;	//四组系数
hilbert_para coef2;
hilbert_para coef3;
hilbert_para coef4;
struct file *g_fpuo = NULL,*g_fpio = NULL;

static const float db1[2] = {
    7.071067811865475244008443621048490392848359376884740365883398e-01,
    7.071067811865475244008443621048490392848359376884740365883398e-01
};

static const float db2[4] = {
    4.829629131445341433748715998644486838169524195042022752011715e-01,
    8.365163037378079055752937809168732034593703883484392934953414e-01,
    2.241438680420133810259727622404003554678835181842717613871683e-01,
    -1.294095225512603811744494188120241641745344506599652569070016e-01
};


static const float db3[6] = {
    3.326705529500826159985115891390056300129233992450683597084705e-01,
    8.068915093110925764944936040887134905192973949948236181650920e-01,
    4.598775021184915700951519421476167208081101774314923066433867e-01,
    -1.350110200102545886963899066993744805622198452237811919756862e-01,
    -8.544127388202666169281916918177331153619763898808662976351748e-02,
    3.522629188570953660274066471551002932775838791743161039893406e-02
};


static const float db4[8] = {
    2.303778133088965008632911830440708500016152482483092977910968e-01,
    7.148465705529156470899219552739926037076084010993081758450110e-01,
    6.308807679298589078817163383006152202032229226771951174057473e-01,
    -2.798376941685985421141374718007538541198732022449175284003358e-02,
    -1.870348117190930840795706727890814195845441743745800912057770e-01,
    3.084138183556076362721936253495905017031482172003403341821219e-02,
    3.288301166688519973540751354924438866454194113754971259727278e-02,
    -1.059740178506903210488320852402722918109996490637641983484974e-02
};


// static const float db5[10] = {
//     1.601023979741929144807237480204207336505441246250578327725699e-01,
//     6.038292697971896705401193065250621075074221631016986987969283e-01,
//     7.243085284377729277280712441022186407687562182320073725767335e-01,
//     1.384281459013207315053971463390246973141057911739561022694652e-01,
//     -2.422948870663820318625713794746163619914908080626185983913726e-01,
//     -3.224486958463837464847975506213492831356498416379847225434268e-02,
//     7.757149384004571352313048938860181980623099452012527983210146e-02,
//     -6.241490212798274274190519112920192970763557165687607323417435e-03,
//     -1.258075199908199946850973993177579294920459162609785020169232e-02,
//     3.335725285473771277998183415817355747636524742305315099706428e-03
// };


// All coif coefficents have to be multiplied by sqrt(2)
static const float coif1[6] = {
    -5.142972847076845595317549230122688830344559947132656813651045e-02,
    2.389297284707684559531754923012268883034455994713265681365104e-01,
    6.028594569415369119063509846024537766068911989426531362730209e-01,
    2.721405430584630880936490153975462233931088010573468637269790e-01,
    -5.142972847076845595317549230122688830344559947132656813651045e-02,
    -1.107027152923154404682450769877311169655440052867343186348954e-02
};


static const float coif2[12] = {
    1.158759673871686817889714882853120395708315073355502818875931e-02,
    -2.932013798346856448679594524397843054053420947418409889774786e-02,
    -4.763959031100813225872995081511549408622753909592460525840745e-02,
    2.730210465347666137982239328923516270034828327990699588033501e-01,
    5.746823938568638472459483149751499367740786490481481391460366e-01,
    2.948671936956191896750637208703777973914107635455611537640778e-01,
    -5.408560709171142997443672832006888537570221990444706777525838e-02,
    -4.202648046077160694657530752545884878978719268926222513485613e-02,
    1.674441016327950635146257083249391698866289538037299820224006e-02,
    3.967883612962012109043447090269950094081810916481648252817197e-03,
    -1.289203356140659543141355500990678257894936161704492503370186e-03,
    -5.095053991076441489598480835620951586540050976664367876412655e-04
};

static const float coif3[18] = {
    -2.682418670922068664584689955153722375535836177157637134187840e-03,
    5.503126707831385107969640263617469178794666057252906037981936e-03,
    1.658356047917034608134280439996549525220639437145367606178002e-02,
    -4.650776447872697640390293095170192691113917841041002855534619e-02,
    -4.322076356021191118175840907244577856782537221435748296465882e-02,
    2.865033352736474630249006862976158896891076238443844211133873e-01,
    5.612852568703300445990941995240077241406247774064453800050914e-01,
    3.029835717728241602862575774374668529867757043461413348549577e-01,
    -5.077014075488886159516471867138370972545857441670871832472707e-02,
    -5.819625076158553022607041679522801089624825903982541419721721e-02,
    2.443409432116695639462954438418928805487699080947974989338820e-02,
    1.122924096203786563399489540091488781245346096838814728167341e-02,
    -6.369601011048822977293753932627342482077585617391852852955559e-03,
    -1.820458915566242322836631665832145136570132777862391313328351e-03,
    7.902051009575939937150950543290226440287715441826917281929124e-04,
    3.296651737931830308416338897758022998655744276957481989605186e-04,
    -5.019277455327664998007173088097694083956570594580641192332170e-05,
    -2.446573425530813115445387662881902303945941576472342106918209e-05
};


static const float coif4[24] = {
    6.309612114309468490753696608619526520153127603444406835368201e-04,
    -1.152225143769973488683007937016166047881572156705066038094891e-03,
    -5.194525163470323267558201363327294331811309729430512113592118e-03,
    1.136246148326482276463392678363118465908960082105224676102131e-02,
    1.886723856956305960822813160712701905823879297781452350370094e-02,
    -5.746424190192718517290527411385172124443396690932404284859269e-02,
    -3.965265296244913762718094206756579981738035770770645437919302e-02,
    2.936674050161006858761278962798582650835466243678172528509866e-01,
    5.531264550395492870333469741987846570947502710783248169642137e-01,
    3.071573096678856987248881030393884808414165269795297009902001e-01,
    -4.711273752389572084912399351781012121935994396763702238263689e-02,
    -6.803811467802056988332974920928626798429778679560269769187728e-02,
    2.781363695846951303169163645831936314699164412528991864702607e-02,
    1.773583142270308388403079552822372238681544967313003044695583e-02,
    -1.075631615508724933047071603601897536695959225169888787867102e-02,
    -4.001010844950535391911552472397083276670126595827549403173754e-03,
    2.652664913530499860820143301690017184933302935238430721089152e-03,
    8.955939276952843603555618778866181384528643960440369133096025e-04,
    -4.165001950941708741516836418852536615951250588002878691463468e-04,
    -1.838296167136253805617482342622910940008368723403836355183423e-04,
    4.408022661597206973006038672236031501663774161685451815597956e-05,
    2.208284691230832960893331999804142845136324572860276715790883e-05,
    -2.304919162676504406778986897925054839632903355820414483306851e-06,
    -1.262179179994622253884862172782890488140153502131112374520603e-06
};



// static const float coif5[30] = {
//     -1.499645228345950331670593167919531667975440598691604525531231e-04,
//     2.535527523580334712936363872191554706055603482812691726895588e-04,
//     1.540286725995222360335148244676269541414659303531250711822333e-03,
//     -2.941078164035693185044038586065593320891475311414770624555173e-03,
//     -7.164112349410053294382279572472252500899544810929605832362178e-03,
//     1.655218330649288840540841623080651353621667424921282557975513e-02,
//     1.991901719798432056056857854066125809443504706772520641876273e-02,
//     -6.499837825472324963374262221660858232544804226063450042795603e-02,
//     -3.680255347446873527191823500872992242220223547780834450868002e-02,
//     2.980959014587191795511466861338063554509597132272839414668911e-01,
//     5.475082713540367154128337935687830970431964302909253422329131e-01,
//     3.097002590784203529311533316221254677074498876376965941549923e-01,
//     -4.386731482823615640442730013366750193381707273908757638050452e-02,
//     -7.464442013283971243472663968192859562973186442054433655531762e-02,
//     2.919469277528073666095772398605275751022315529465178441510318e-02,
//     2.310457227706684192610065243663928370022983285246219996141160e-02,
//     -1.397129268638200558584119246355879336305763752871371182932059e-02,
//     -6.476749751505861835547590642967453082384538848552165075614441e-03,
//     4.781116799130657606400088024549264921093190305150784065791191e-03,
//     1.719383484385504023022397097446276782318002683055773803854075e-03,
//     -1.174947934413537690027670037110105795928147523549002426409332e-03,
//     -4.508222400696236312231932151038336110220594834213702970043431e-04,
//     2.134457975086291667348984871136041914777578046177470626552867e-04,
//     9.924691139873533169989496559631669037970741600337089424730635e-05,
//     -2.914684388622130824599478843558087403539428940986384077972155e-05,
//     -1.504031798197685905639227292876711236513927746903476131955063e-05,
//     2.616809660013118152124234488302931243021794024318439103773996e-06,
//     1.457502921355163070577152619048168436286350537937563166257584e-06,
//     -1.148199649902979726237655584441763456854312591680755421569962e-07,
//     -6.791060677322355511541065559242475254516249773485524025251102e-08
// };

static const float sym2[4] = {
    0.48296291314469025, 0.83651630373746899,
    0.22414386804185735, -0.12940952255092145 };


static const float sym3[6] = {
    0.33267055295095688, 0.80689150931333875,
    0.45987750211933132, -0.13501102001039084,
    -0.085441273882241486, 0.035226291882100656 };

static const float sym4[8] = {
    0.032223100604042702, -0.012603967262037833,
    -0.099219543576847216, 0.29785779560527736,
    0.80373875180591614, 0.49761866763201545,
    -0.02963552764599851, -0.075765714789273325 };

// static const float sym5[10] = {
//     0.019538882735286728, -0.021101834024758855,
//     -0.17532808990845047, 0.016602105764522319,
//     0.63397896345821192, 0.72340769040242059,
//     0.1993975339773936, -0.039134249302383094,
//     0.029519490925774643, 0.027333068345077982 };

/* 要求:
* (1)C0^2 + C1^2 + C2^2 + C3^2 = 1
* (2)C1 * C3 + C0 * C2 = 0
*/

/* 矩阵变换需要用到的参数 */
static Filter DBn = {
	.Len = 4,
	.C0 = (const float) 0.482962913144690,
	.C1 = (const float) 0.836516303737469,
	.C2 = (const float) 0.224143868041857,
	.C3 = (const float) -0.129409522550921,
};


/* 按分解矩阵进行分解 */
static void transform(float *input_Data, float *output_Data, int length)
{
    int i=0;
    int count=0;

	for(i=0; i<=(length-3); i=i+1)
    {
        output_Data[count++] = DBn.C0 * input_Data[i] + DBn.C1 * input_Data[i+1] + DBn.C2 * input_Data[i+2] + DBn.C3 * input_Data[i+3];
    }
    output_Data[count++] = DBn.C2 * input_Data[0] + DBn.C3 * input_Data[1] + DBn.C0 * input_Data[length-2] + DBn.C1 * input_Data[length-1];
    output_Data[count] = DBn.C1 * input_Data[0] - DBn.C0 * input_Data[1] + DBn.C3 * input_Data[length-2] - DBn.C2 * input_Data[length-1];
}


/******************************************************************************
*	功能：利用矩阵进行分解
*	@srcData:	源信号;
*	@resData:	保存结果;
*	@srcLen:	源信号长度;
*******************************************************************************/
void DWT_Matrix_Transform(float *srcData, s16 *resData, int srcLen)
{
	
	int i=0;
    float out1Data[srcLen];

	memset(out1Data, 0, srcLen);
	transform(srcData, out1Data, srcLen);

    for (i = 0; i < srcLen; i +=1)			
    {
		resData[i] = (s16)rt_round(out1Data[i]);	//此处有风险，小心处理
		// *((s16 *)&resData[i]) = (s16)(out1Data[i]);	//此处有风险，小心处理
    }

	return ;
}

/******************************************************************************
函数名称：	hilbert_para_init
函数版本：	01.01
创建作者：	ygl
创建日期：	2021.11.17
函数说明：	hilbert初始化。
参数说明：	无
返回值：  	长度
修改记录：	
*******************************************************************************/
void _hilbert_para_init(void)	//初始化系数
{
    memset(&coef1, 0, sizeof(coef1));
    memset(&coef2, 0, sizeof(coef2));
    memset(&coef3, 0, sizeof(coef3));
    memset(&coef4, 0, sizeof(coef4));			
}
/******************************************************************************
函数名称：	filter_Init
函数版本：	01.01
创建作者：	ygl
创建日期：	2021.11.17
函数说明：	卷积核初始化。
参数说明：	无
返回值：  	长度
修改记录：	
*******************************************************************************/
void filter_Init(void)
{
	#if 0
	g_fpuo = rt_file_open("/app/uo.txt", O_CREAT|O_RDWR,0);		//近似系数存储文件
	if(IS_ERR(g_fpuo))
	{
		rt_printf("LINE %d write file error!\n",__LINE__);
		rt_printf("write file error!\n");
		// return -1;
	}
	g_fpio = rt_file_open("/app/io.txt", O_CREAT|O_RDWR,0);		//近似系数存储文件
	if(IS_ERR(g_fpio))
	{
		rt_printf("LINE %d write file error!\n",__LINE__);
		rt_printf("write file error!\n");
		// return -1;
	}
	#endif
	_hilbert_para_init();	
}

float hilbert_filter(float *channel1, float *channel2, int len)
{
		int i=0;
		const float c1=2.0216;
		const float c2=1.2562;
		const float c3=0.2276;
		const float c4=1.5219;
		const float c5=0.6205;
		const float c6=0.0463;
		
		// const float c1=2.1441;
		// const float c2=1.4380;
		// const float c3=0.2904;
		// const float c4=1.6444;
		// const float c5=0.7421;
		// const float c6=0.0658;	
		
        float channel_f1 = 0;
        float channel_f2 = 0;
        float channel_f3 = 0;
        float channel_f4 = 0;
		float Qval = 0;

		memset(&coef1,0,sizeof(coef1));
		memset(&coef2,0,sizeof(coef2));
		memset(&coef3,0,sizeof(coef3));
		memset(&coef4,0,sizeof(coef4));
		
		for(i=0;i<len;i++)
		{
		    channel_f1 = c1*(coef1.y2)-c2*(coef1.y4)+c3*(coef1.y6)+c3*(coef1.x1)-c2*(coef1.x3)+c1*(coef1.x5)-coef1.x7;	
                 
		    coef1.y6 = coef1.y5;
		    coef1.y5 = coef1.y4;
		    coef1.y4 = coef1.y3;
		    coef1.y3 = coef1.y2;
		    coef1.y2 = coef1.y1;
		    coef1.y1 = channel_f1;
		    coef1.x7 = coef1.x6;
		    coef1.x6 = coef1.x5;
		    coef1.x5 = coef1.x4;
		    coef1.x4 = coef1.x3;
		    coef1.x3 = coef1.x2;
		    coef1.x2 = coef1.x1;
		    coef1.x1 = channel1[i];  
		    
		    
		    channel_f2 = c4*(coef2.y2)-c5*(coef2.y4)+c6*(coef2.y6)+c6*(channel1[i])-c5*(coef2.x2)+c4*(coef2.x4)-(coef2.x6);     
		    coef2.y6 = coef2.y5;
		    coef2.y5 = coef2.y4;
		    coef2.y4 = coef2.y3;
		    coef2.y3 = coef2.y2;
		    coef2.y2 = coef2.y1;
		    coef2.y1 = channel_f2;		   
		    coef2.x6 = coef2.x5;
		    coef2.x5 = coef2.x4;
		    coef2.x4 = coef2.x3;
		    coef2.x3 = coef2.x2;
		    coef2.x2 = coef2.x1;
		    coef2.x1 = channel1[i];    
		    
		    channel_f3 = c1*(coef3.y2)-c2*(coef3.y4)+c3*(coef3.y6)+c3*(coef3.x1)-c2*(coef3.x3)+c1*(coef3.x5)-coef3.x7;	
                    coef3.y6 = coef3.y5;
		    coef3.y5 = coef3.y4;
		    coef3.y4 = coef3.y3;
		    coef3.y3 = coef3.y2;
		    coef3.y2 = coef3.y1;
		    coef3.y1 = channel_f3;
		    coef3.x7 = coef3.x6;
		    coef3.x6 = coef3.x5;
		    coef3.x5 = coef3.x4;
		    coef3.x4 = coef3.x3;
		    coef3.x3 = coef3.x2;
		    coef3.x2 = coef3.x1;
		    coef3.x1 = channel2[i];  
		    
		    
		    channel_f4 = c4*(coef4.y2)-c5*(coef4.y4)+c6*(coef4.y6)+c6*(channel2[i])-c5*(coef4.x2)+c4*(coef4.x4)-(coef4.x6);     
		    coef4.y6 = coef4.y5;
		    coef4.y5 = coef4.y4;
		    coef4.y4 = coef4.y3;
		    coef4.y3 = coef4.y2;
		    coef4.y2 = coef4.y1;
		    coef4.y1 = channel_f4;		   
		    coef4.x6 = coef4.x5;
		    coef4.x5 = coef4.x4;
		    coef4.x4 = coef4.x3;
		    coef4.x3 = coef4.x2;
		    coef4.x2 = coef4.x1;
		    coef4.x1 = channel2[i];  
		    
		    Qval += (channel_f1*channel_f4 - channel_f2*channel_f3)/2;
                    
		}
		Qval =(float) Qval/len; 
		// rt_printf("len = %d,qval = %f\r\n",len,Qval);
		return Qval;
}


int filtlength(const char* name) {
    int len = strlen(name);
    int i = 0;
    char *new_str = NULL;
    int N = 0;
	if (!strcmp(name,"haar") || !strcmp(name,"db1")) {
		return 2;
	}
    else if (len > 2 && strstr(name, "db") != NULL)
    {
        new_str = (char*)rt_malloc(sizeof(char)*(len-2 + 1));
        for (i = 2; i < len + 1; i++)
            new_str[i - 2] = name[i];

        N = atoi(new_str);
        rt_free(new_str);
        if (N>38)
        {
            rt_printf("\n Filter Not in Database \n");
            return -1;
        }

        return N * 2;
    }

    else if (len > 4 && strstr(name, "coif") != NULL)
    {
        new_str = (char*)rt_malloc(sizeof(char)*(len - 4 + 1));
        for (i = 4; i < len + 1; i++)
            new_str[i - 4] = name[i];

        N = atoi(new_str);
        rt_free(new_str);
        if (N>17)
        {
            rt_printf("\n Filter Not in Database \n");
            return -1;
        }

        return N * 6;
    }
    else if (len > 3 && strstr(name, "sym") != NULL)
    {
        new_str = (char*)rt_malloc(sizeof(char)*(len - 3 + 1));
        for (i = 3; i < len + 1; i++)
            new_str[i - 3] = name[i];

        N = atoi(new_str);
        rt_free(new_str);
        if (N>20 || N < 2)
        {
            rt_printf("\n Filter Not in Database \n");
            return -1;
        }

        return N * 2;
    }
	else {
		rt_printf("\n Filter Not in Database \n");
		return -1;
	}

}

void copy_reverse(const float *in, int N,float *out)
{
    int count = 0;
    for (count = 0; count < N; count++)
        out[count] = in[N - count - 1];
}

void qmf_even(const float *in, int N,float *out)
{
    int count = 0;
    for (count = 0; count < N; count++)
    {
        out[count] = in[N - count - 1];
        if (count % 2 != 0)
        {
            out[count] = -1 * out[count];
        }
    }
}

void qmf_wrev(const float *in, int N, float *out)
{
    float *sigOutTemp;
    sigOutTemp = (float*)rt_malloc(N*sizeof(float));

    qmf_even(in, N, sigOutTemp);
    copy_reverse(sigOutTemp, N, out);

    rt_free(sigOutTemp);
    return;
}

void wave_copy(const float *in, int N, float *out)
{
    int count = 0;
    for (count = 0; count < N; count++)
        out[count] = in[count];
}

int filtcoef(const char* name, float *lp1, float *hp1, float *lp2, float *hp2) {
    int i = 0; 
    int N = filtlength(name);
	if (!strcmp(name,"haar") || !strcmp(name,"db1")) {
        copy_reverse(db1, N, lp1);
        qmf_wrev(db1, N, hp1);
        wave_copy(db1, N, lp2);
        qmf_even(db1, N, hp2);

		return N;
	}
	else if (!strcmp(name,"db2")){
        copy_reverse(db2, N, lp1);
        qmf_wrev(db2, N, hp1);
        wave_copy(db2, N, lp2);
        qmf_even(db2, N, hp2);

        return N;
	}
	else if (!strcmp(name,"db3")){
        copy_reverse(db3, N, lp1);
        qmf_wrev(db3, N, hp1);
        wave_copy(db3, N, lp2);
        qmf_even(db3, N, hp2);

        return N;
	}
	else if (!strcmp(name,"db4")){
        copy_reverse(db4, N, lp1);
        qmf_wrev(db4, N, hp1);
        wave_copy(db4, N, lp2);
        qmf_even(db4, N, hp2);

        return N;
	}

	else if (!strcmp(name,"coif1")){
        float *coeffTemp;
        coeffTemp = (float*)rt_malloc(N*sizeof(float));

        wave_copy(coif1, N, coeffTemp);
        for (i = 0; i < N; ++i) {
            coeffTemp[i] *= M_SQRT2;
        }

        copy_reverse(coeffTemp, N, lp1);
        qmf_wrev(coeffTemp, N, hp1);
        wave_copy(coeffTemp, N, lp2);
        qmf_even(coeffTemp, N, hp2);
        rt_free(coeffTemp);

        return N;
	}
	else if (!strcmp(name,"coif2")){
        float *coeffTemp;
        coeffTemp = (float*)rt_malloc(N*sizeof(float));

        wave_copy(coif2, N, coeffTemp);
        for (i = 0; i < N; ++i) {
            coeffTemp[i] *= M_SQRT2;
        }

        copy_reverse(coeffTemp, N, lp1);
        qmf_wrev(coeffTemp, N, hp1);
        wave_copy(coeffTemp, N, lp2);
        qmf_even(coeffTemp, N, hp2);
        rt_free(coeffTemp);

        return N;
	}
	else if (!strcmp(name,"coif3")){
        float *coeffTemp;
        coeffTemp = (float*)rt_malloc(N*sizeof(float));

        wave_copy(coif3, N, coeffTemp);
        for (i = 0; i < N; ++i) {
            coeffTemp[i] *= M_SQRT2;
        }

        copy_reverse(coeffTemp, N, lp1);
        qmf_wrev(coeffTemp, N, hp1);
        wave_copy(coeffTemp, N, lp2);
        qmf_even(coeffTemp, N, hp2);
        rt_free(coeffTemp);

        return N;
	}
	else if (!strcmp(name,"coif4")){
        float *coeffTemp;
        coeffTemp = (float*)rt_malloc(N*sizeof(float));

        wave_copy(coif4, N, coeffTemp);
        for (i = 0; i < N; ++i) {
            coeffTemp[i] *= M_SQRT2;
        }

        copy_reverse(coeffTemp, N, lp1);
        qmf_wrev(coeffTemp, N, hp1);
        wave_copy(coeffTemp, N, lp2);
        qmf_even(coeffTemp, N, hp2);
        rt_free(coeffTemp);

        return N;
	}

	else if (!strcmp(name,"sym2") || !strcmp(name,"sym1")){
        copy_reverse(sym2, N, lp1);
        qmf_wrev(sym2, N, hp1);
        wave_copy(sym2, N, lp2);
        qmf_even(sym2, N, hp2);
        return N;
	}

	else if (!strcmp(name,"sym3")){
        copy_reverse(sym3, N, lp1);
        qmf_wrev(sym3, N, hp1);
        wave_copy(sym3, N, lp2);
        qmf_even(sym3, N, hp2);
        return N;
	}

	else if (!strcmp(name,"sym4")){
        copy_reverse(sym4, N, lp1);
        qmf_wrev(sym4, N, hp1);
        wave_copy(sym4, N, lp2);
        qmf_even(sym4, N, hp2);
        return N;
	}

	else {
		rt_printf("\n Filter Not in Database \n");
		return -1;
	}

	return 0;
}


wave_object wave_init(const char* wname) {
	wave_object obj = NULL;
	int retval;
	retval = 0;

	if (wname != NULL) {
		retval = filtlength(wname);
	}

	obj = (wave_object)rt_malloc(sizeof(struct wave_set) + sizeof(float)* 4 * retval);

	obj->filtlength = retval;
	obj->lpd_len = obj->hpd_len = obj->lpr_len = obj->hpr_len = obj->filtlength;
	strcpy(obj->wname, wname);
	if (wname != NULL) {
		filtcoef(wname,obj->params,obj->params+retval,obj->params+2*retval,obj->params+3*retval);
	}
	obj->lpd = &obj->params[0];
	obj->hpd = &obj->params[retval];
	obj->lpr = &obj->params[2 * retval];
	obj->hpr = &obj->params[3 * retval];
	return obj;
}

wt_object wt_init(wave_object wave,const char* method, int siglength,int J) {
	int i;
	wt_object obj = NULL;

	if (J > 100) {
		rt_printf("\n The Decomposition Iterations Cannot Exceed 100. Exiting \n");
		return NULL;
	}

	if (!strcmp(method, "modwt") || !strcmp(method, "MODWT")) {

		if (!strstr(wave->wname,"haar")) {
			if (!strstr(wave->wname,"db")) {
				if (!strstr(wave->wname, "sym")) {
					if (!strstr(wave->wname, "coif")) {
						rt_printf("\n MODWT is only implemented for orthogonal wavelet families - db, sym and coif \n");
						return NULL;
					}
				}
			}
		}

		// obj = (wt_object)rt_malloc(sizeof(struct wt_set) + sizeof(float)* (siglength * 2 * (J + 1)));
		// obj->outlength = siglength * (J + 1); // Default
		obj = (wt_object)rt_malloc(sizeof(struct wt_set) + sizeof(float)* (siglength * 2));
		obj->outlength = siglength ; // 
		strcpy(obj->ext, "per"); // Default
	}

	obj->wave = wave;
	obj->siglength = siglength;
	obj->modwtsiglength = siglength;
	obj->J = J;
	// obj->MaxIter = MaxIter;
	strcpy(obj->method, method);

	if (siglength % 2 == 0) {
		obj->even = 1;
	}
	else {
		obj->even = 0;
	}

	obj->cobj = NULL;

	strcpy(obj->cmethod, "direct"); // Default
	obj->cfftset = 0;
	obj->lenlength = J + 2;
	obj->output = &obj->params[0];

	if (!strcmp(method, "modwt") || !strcmp(method, "MODWT")) {
		for (i = 0; i < siglength * 2; ++i) {
			obj->params[i] = 0.0;
		}
	}

	return obj;
}

static void modwt_per(wt_object wt, int M, float *inp, float *cA, int len_cA) {
	int l, i, t, len_avg;
	float s;
	float *filt;
	len_avg = wt->wave->lpd_len;

	filt = (float*)rt_malloc(sizeof(float)* len_avg);
	s = M_SQRT2;//sqrt(2.0);
	for (i = 0; i < len_avg; ++i) {
		filt[i] = wt->wave->lpd[i] / s;
	}

	for (i = 0; i < len_cA; ++i) {
		t = i;
		cA[i] = filt[0] * inp[t];

		for (l = 1; l < len_avg; l++) {
			t -= M;
			while (t >= len_cA) {
				t -= len_cA;
			}
			while (t < 0) {
				// t += len_cA; //ygl 传入波形数据，有可能首尾不对称。此处首核卷积改为第一点开始    
                t += M;
			}

			cA[i] += filt[l] * inp[t];
		}
	}
	rt_free(filt);
}

static void modwt_direct(wt_object wt, const float *inp) {
	int i, J, temp_len, iter, M;
	float  *cA;


	if (strcmp(wt->ext, "per")) {
		rt_printf("MODWT direct method only uses periodic extension per. \n");
		rt_printf(" Use MODWT fft method for symmetric extension sym \n");
		return;
	}

	temp_len = wt->siglength;
	J = wt->J;
	wt->length[0] = wt->length[J] = temp_len;
	wt->outlength = wt->length[J + 1] =  temp_len;	//只考虑近似值
	M = 1;
	for (iter = 1; iter < J; ++iter) {
		M = 2 * M;
		wt->length[iter] = temp_len;
	}


	cA = (float*)rt_malloc(sizeof(float)* temp_len);
	M = 1;

	for (i = 0; i < temp_len; ++i) {
		wt->params[i] = inp[i];
	}

	for (iter = 0; iter < J; ++iter) {
		if (iter > 0) {
			M = 2 * M;
		}

		modwt_per(wt, M, wt->params, cA, temp_len);


		for (i = 0; i < temp_len; ++i) {
			wt->params[i] = cA[i];
		}

	}

	rt_free(cA);

}

void modwt(wt_object wt, const float *inp) {
	if (!strcmp(wt->cmethod, "direct")) {
		modwt_direct(wt, inp);
	}
	else {
		rt_printf("Error- Available Choices for this method are - direct and fft \n");
		return;
	}
}

void wave_free(wave_object object) {
	rt_free(object);
}

void wt_free(wt_object object) {
	rt_free(object);
}


/*------------------------------ 内部函数 -------------------------------------
内部函数以下划线‘_’开头，不需要检查参数的合法性.
*/


/*------------------------------ 测试函数 -------------------------------------
一个实体文件必须带一个本模块的测试函数来进行单元测试，如果的确不方便在本模块中
进行单元测试，必须在此注明实际的测试位置（例如在哪个实体文件中使用哪个测试函数）.
*/

// #ifndef __KERNEL__ 
#if 0
/**************************************************************************
函数名称: main()
函数版本：1.00 
作者：    ygl
创建日期：2022.03.10
函数功能说明：小波测试函数，在vscode当前文件下测试。使用D:\\MinGW\\bin\\gcc.exe，需要配置vscode相关"launch.json"和"task.json“
            支持：haar  db1~5   sym2~5  coif1~5 共四类小波函数
输入参数：
输出参数：
返回值： 
***************************************************************************/
int main() {
	wave_object obj;
	wt_object wt;
	float *inp, *out, *diff;
	int N, i, J;

	FILE *ifp;
    FILE *fpAppCoef = NULL;
	float temp[10000];

	char *name = "haar";
	obj = wave_init(name);
	// wave_summary(obj);

	// ifp = fopen("signal.txt", "r");
    ifp = fopen("myDWTData4.txt", "r");
    fseek(ifp,0,SEEK_END);
	N = ftell(ifp)/7;
	fseek(ifp,0,SEEK_SET);
	i = 0;
	if (!ifp) {
		printf("Cannot Open File");
		exit(100);
	}
	while (!feof(ifp)) {
        fscanf(ifp, "%f \n", &temp[i]);
		i++;
	}

    //***********************************把注释取消掉，可以将变换结果存于文件****************************************
	if(NULL == (fpAppCoef = fopen("DWT_AppData.txt", "w")))			//近似系数存储文件
	{
		printf("write file error!\n");
		return -1;
	}
    //***********************************************************************************************************/
	// N = 177;
	
	fclose(ifp);

	inp = (float*)malloc(sizeof(float)* N);
	out = (float*)malloc(sizeof(float)* N);
	diff = (float*)malloc(sizeof(float)* N);
	//wmean = mean(temp, N);

	for (i = 0; i < N; ++i) {
		inp[i] = temp[i];
		//printf("%g \n",inp[i]);
	}
	J = 4;

	wt = wt_init(obj, "modwt", N, J);// Initialize the wavelet transform object
	
	modwt(wt, inp);// Perform MODWT
	//MODWT output can be accessed using wt->output vector. Use wt_summary to find out how to extract appx and detail coefficients
	
	for (i = 0; i < wt->outlength; ++i) {
		printf("%g ",wt->output[i]);
        fprintf(fpAppCoef, "%f\t", wt->output[i]);		//将数据流写入文件
	}

	// imodwt(wt, out);// Perform ISWT (if needed)
	// Test Reconstruction

	// for (i = 0; i < wt->siglength; ++i) {
	// 	diff[i] = out[i] - inp[i];
	// }


    fclose(fpAppCoef);
	wave_free(obj);
	wt_free(wt);

	rt_free(inp);
	rt_free(out);
	rt_free(diff);
	return 0;
}

#endif
/*------------------------------ 文件结束 -------------------------------------
*/