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
/*------------------------------ 文件结束 -------------------------------------
 */