#include<stdio.h>
#include<math.h>
#include<stdlib.h>
//#include<malloc.h>

float* convolution_HPF(float *, float *, int, int);  //Prototype for convolution of data with HPF coefficient
float* convolution_LPF(float *, float *, int, int);  //Prototype for convolution of data with LPF coefficient
float* downsample(float *, int);                     //Prototype for downsampling the sequence
float mean_abs(float *, int);                        //Prototype for calculating mean of detail coefficeints at each decomposition stage
float linreg(float [], int);

#define L 201                                       //Length of input data
#define N floor(log10(L)/log10(2))                     //No. of Decomposition levels

//Filter coefficients for Daubechies16 wavelet
static float LPF_COEF[32] = {0.0031892209253436892,
    0.03490771432362905,
    0.1650642834886438,
    0.43031272284545874,
    0.6373563320829833,
    0.44029025688580486,
    -0.08975108940236352,
    -0.3270633105274758,
    -0.02791820813292813,
    0.21119069394696974,
    0.027340263752899923,
    -0.13238830556335474,
    -0.006239722752156254,
    0.07592423604445779,
    -0.007588974368642594,
    -0.036888397691556774,
    0.010297659641009963,
    0.013993768859843242,
    -0.006990014563390751,
    -0.0036442796214883506,
    0.00312802338120381,
    0.00040789698084934395,
    -0.0009410217493585433,
    0.00011424152003843815,
    0.00017478724522506327,
    -6.103596621404321e-05,
    -1.394566898819319e-05,
    1.133660866126152e-05,
    -1.0435713423102517e-06,
    -7.363656785441815e-07,
    2.3087840868545578e-07,
    -2.1093396300980412e-08};

static float HPF_COEF[32] = {-2.1093396300980412e-08,
    -2.3087840868545578e-07,
    -7.363656785441815e-07,
    1.0435713423102517e-06,
    1.133660866126152e-05,
    1.394566898819319e-05,
    -6.103596621404321e-05,
    -0.00017478724522506327,
    0.00011424152003843815,
    0.0009410217493585433,
    0.00040789698084934395,
    -0.00312802338120381,
    -0.0036442796214883506,
    0.006990014563390751,
    0.013993768859843242,
    -0.010297659641009963,
    -0.036888397691556774,
    0.007588974368642594,
    0.07592423604445779,
    0.006239722752156254,
    -0.13238830556335474,
    -0.027340263752899923,
    0.21119069394696974,
    0.02791820813292813,
    -0.3270633105274758,
    0.08975108940236352,
    0.44029025688580486,
    -0.6373563320829833,
    0.43031272284545874,
    -0.1650642834886438,
    0.03490771432362905,
    -0.0031892209253436892};

main()
{
	float data[L]={19,
23,
31,
18,
22,
34,
30,
37,
62,
28,
40,
47,
28,
26,
43,
46,
20,
68,
62,
29,
52,
33,
30,
35,
40,
5,
28,
41,
18,
35,
21,
16,
14,
28,
15,
29,
33,
16,
21,
9,
17,
20,
23,
24,
33,
12,
22,
31,
34,
7,
22,
23,
21,
42,
10,
15,
20,
15,
15,
29,
14,
9,
11,
32,
7,
30,
39,
34,
32,
31,
40,
23,
35,
63,
28,
44,
23,
29,
22,
12,
17,
23,
6,
9,
10,
30,
8,
8,
27,
40,
32,
23,
31,
14,
28,
21,
20,
28,
27,
15,
22,
24,
31,
14,
43,
15,
14,
16,
16,
3,
29,
8,
11,
16,
17,
46,
25,
25,
39,
18,
30,
17,
16,
24,
14,
14,
14,
24,
31,
10,
13,
27,
14,
52,
36,
26,
27,
25,
32,
11,
33,
41,
23,
52,
28,
37,
41,
14,
21,
17,
16,
7,
24,
7,
12,
25,
7,
6,
44,
9,
18,
41,
26,
26,
31,
17,
10,
30,
14,
9,
37,
30,
11,
12,
33,
17,
8,
9,
46,
7,
6,
5,
31,
3,
15,
12,
19,
23,
53,
21,
41,
19,
52,
16,
39,
22,
7,
7,
40,
19,
9};
	float *ccout_HPF, *ccout_LPF, *dsout_HPF, *dsout_LPF, *temp1, *temp2, mean[11], slope;
	int i, LENGTH, N_1, k=0, r, l=1, prv_LENGTH;

	N_1 = (int) N;
	printf("N=%d\n", N_1);

	//1st-level decomposition
	r=pow(2,k);
	LENGTH = L/r;    //Length of output sequence of convolution at each step
	if(L<32)
		LENGTH = 32;
	printf("%d\n", LENGTH);

	ccout_HPF = (float *) calloc(sizeof(float), LENGTH);        //ptr to Array storing values of convolution with HPF
	dsout_HPF = (float *) calloc(sizeof(float), (LENGTH/2));   //ptr to Array storing values of downsample output
	ccout_LPF = (float *) calloc(sizeof(float), LENGTH);        //ptr to Array storing values of convolution with LPF
	dsout_LPF = (float *) calloc(sizeof(float), (LENGTH/2));   //ptr to Array storing values of downsample output


	printf("%dst Convolution output with HPF: \n", k+1);
	ccout_HPF = convolution_HPF(data, HPF_COEF, L/r, 32);
	for(i=0; i<LENGTH; i++)
		printf("y[%d] = %f\n", i, *(ccout_HPF+i));

	printf("%dst Downsampled output: \n", k+1);
	dsout_HPF = downsample(ccout_HPF, LENGTH);
	for(i=0; i<LENGTH/2; i++)
		printf("y[%d] = %f\n", i, *(dsout_HPF+i));

	mean[0] = mean_abs(dsout_HPF, LENGTH/2);
	printf("Mean=%f\n", mean[0]);

	printf("%dst Convolution output with LPF: \n", k+1);
	ccout_LPF = convolution_LPF(data, LPF_COEF, L/r, 32);
	for(i=0; i<LENGTH; i++)
		printf("y[%d] = %f\n", i, *(ccout_LPF+i));

	printf("%dst Downsampled output: \n", k+1);
	dsout_LPF = downsample(ccout_LPF, LENGTH);
	//for(i=0; i<LENGTH/2; i++)
		//printf("y[%d] = %f\n", i, *(dsout_LPF+i));

	temp1 = (float *) calloc(sizeof(float), LENGTH/2);      //pointer to Array storing values of approximation at each step
    for(i=0; i<(LENGTH/2); i++)
			*(temp1+i) = *(dsout_LPF+i);
	//for(i=0; i<LENGTH/2; i++)
		//printf("y[%d] = %f\n", i, *(temp1+i));
	prv_LENGTH = LENGTH/2;

	free(ccout_HPF);
	free(ccout_LPF);
	free(dsout_HPF);
	free(dsout_LPF);
	//End of 1st level decomposition

	for(k=1, l=1; k<N_1&&l<N_1; k++, l++)                           //2nd to N_1th stage decomposition
	{
		//k=k+1;
		if((LENGTH/2)<=32)
			LENGTH = 32;
		else
		{
			//r=pow(2,k);
		    LENGTH = LENGTH/2;
		}
		printf("%d\n", LENGTH);

		temp2 = (float *) calloc(sizeof(float), LENGTH);
		ccout_HPF = (float *) calloc(sizeof(float), LENGTH);        //ptr to Array storing values of convolution with HPF
		dsout_HPF = (float *) calloc(sizeof(float), (LENGTH/2));    //ptr to Array storing values of downsample output
		ccout_LPF = (float *) calloc(sizeof(float), LENGTH);        //ptr to Array storing values of convolution with LPF
		dsout_LPF = (float *) calloc(sizeof(float), (LENGTH/2));   //ptr to Array storing values of downsample output

		for(i=0; i<prv_LENGTH; i++)
			*(temp2+i) = *(temp1+i);
		for(i=25; i<LENGTH; i++)
			*(temp2+i) = 0;

		printf("%dth Convolution output with HPF: \n", k+1);
		ccout_HPF = convolution_HPF(temp2, HPF_COEF, LENGTH, 32);
		//for(i=0; i<LENGTH; i++)
			//printf("y[%d] = %f\n", i, *(ccout_HPF+i));

		printf("%dth Downsampled output: \n", k+1);
		dsout_HPF = downsample(ccout_HPF, LENGTH);
		for(i=0; i<LENGTH/2; i++)
			printf("y[%d] = %f\n", i, *(dsout_HPF+i));

		mean[l] = mean_abs(dsout_HPF, LENGTH/2);
		printf("Mean=%f\n", mean[l]);

		printf("%dth Convolution output with LPF: \n", k+1);
		ccout_LPF = convolution_LPF(temp2, LPF_COEF, LENGTH, 32);
		//for(i=0; i<LENGTH; i++)
			//printf("y[%d] = %f\n", i, *(ccout_LPF+i));

		printf("%dth Downsampled output: \n", k+1);
		dsout_LPF = downsample(ccout_LPF, LENGTH);
		//for(i=0; i<LENGTH/2; i++)
			//printf("y[%d] = %f\n", i, *(dsout_LPF+i));

		for(i=0; i<LENGTH/2; i++)
				*(temp1+i) = *(dsout_LPF+i);
		prv_LENGTH = LENGTH/2;

		//printf("\n");
		//for(i=0; i<LENGTH/2; i++)
			//printf("y[%d] = %f\n", i, *(temp1+i));

		free(ccout_HPF);
		free(ccout_LPF);
		free(dsout_HPF);
		free(dsout_LPF);
		free(temp2);

	}

	for(i=LENGTH/2; i<L; i++)
		*(temp1+i) = 0;
	mean[l] = mean_abs(temp1, LENGTH/2);
		printf("Mean=%f\n", mean[l]);

	for(i=0; i<=N; i++)
		printf("Mean = %f\t", mean[i]);

	slope = linreg(mean, N_1);                                                        //Slope calculation using Linear regression method
	printf("%f\n", slope);
	printf("Slope = H + 0.5,\tor H=1-Beta/2 \t H = %f\n", 1-(slope/2));

	/*slope = (log10(mean[10])-log10(mean[0]))/(log10(1024.0)-log10(1.0));            // Slope and H calculation using simple m = ((y2-y1)/(x2-x1))method

	printf("Slope = %f\n", slope);
	printf("Slope = H + 0.5,\tor H=1-Beta/2 \t H = %f\n", 1-(slope/2));*/

	free(temp1);

}

float *convolution_HPF(float *x, float *h, int n, int m)
{
	float *y, *a, *x2;
	int i, j, k;

	if(n-m !=0)
	{
		if(n>m)
		{
			for(i=m; i<n; i++)
			{
				*(h+i) = 0;
				m = n;
			}
		}
		else
		{
			for(i=n; i<m; i++)
			{
				*(x+i) = 0;
				n = m;
			}
		}
	}

	y = (float *)calloc(sizeof(float), (n));
	a = (float *)calloc(sizeof(float), (n));
	x2 = (float *)calloc(sizeof(float), (n));


	*(y+0) = 0;
	*(a+0) = *(h+0);

	for(j=1;j<n;j++)            /*folding h(n) to h(-n)*/
		*(a+j) = *(h+(n-j));

    /*Circular convolution*/
    for(i=0;i<n;i++)
		*(y+0) = *(y+0) + (*(x+i) * (*(a+i)));

	for(k=1;k<n;k++)
	{
		*(y+k) = 0;                    /*circular shift*/

		for(j=1;j<n;j++)
					*(x2+j)=*(a+(j-1));
					*(x2+0)=*(a+(n-1));
		for(i=0;i<n;i++)
		{
					*(a+i)=*(x2+i);
					*(y+k) = *(y+k) + (*(x+i)*(*(x2+i)));
		}
    }

	return y;

	free(y);
	free(a);
	free(x2);
}

float *convolution_LPF(float *x, float *h, int n, int m)
{
	float *y, *a, *x2;
	int i, j, k;

	if(n-m !=0)
	{
		if(n>m)
		{
			for(i=m; i<n; i++)
			{
				*(h+i) = 0;
				m = n;
			}
		}
		else
		{
			for(i=n; i<m; i++)
			{
				*(x+i) = 0;
				n = m;
			}
		}
	}

	y = (float *)calloc(sizeof(float), (n));
	a = (float *)calloc(sizeof(float), (n));
	x2 = (float *)calloc(sizeof(float), (n));


	*(y+0) = 0;
	*(a+0) = *(h+0);

	for(j=1;j<n;j++)            /*folding h(n) to h(-n)*/
		*(a+j) = *(h+(n-j));

    /*Circular convolution*/
    for(i=0;i<n;i++)
		*(y+0) = *(y+0) + (*(x+i) * (*(a+i)));

	for(k=1;k<n;k++)
	{
		*(y+k) = 0;                    /*circular shift*/

		for(j=1;j<n;j++)
					*(x2+j)=*(a+(j-1));
					*(x2+0)=*(a+(n-1));
		for(i=0;i<n;i++)
		{
					*(a+i)=*(x2+i);
					*(y+k) = *(y+k) + (*(x+i)*(*(x2+i)));
		}
    }

	return y;

	free(y);
	free(a);
	free(x2);

}

float * downsample(float *p, int len)
{
	int i, j;
	//float p1[1];
	float *ds;

	ds = (float *) malloc(sizeof(float)*(len/2));

	for(i=0,j=0; i<len/2; i++,j++)
		*(ds+i) = *(p+i+j);

	return ds;

	free(ds);
}

float mean_abs(float *p, int length)
{
	int i=0;
	float sum=0, mean;

	for(i=0; i<length; i++)
	{
		if(*(p+i) >= 0)
			sum = sum + (*(p+i));
		else
			sum = sum +(*(p+i)*(-1));
	}

	mean = (float) sum/length;

	return mean;
}

float linreg(float p[], int N_1)
{
	int i;
	float sum_x=0, sum_y=0, sum_xy=0, sum_x2=0, m, temp=0;

	for(i=0; i<=N_1; i++)
		sum_x = sum_x + log10(pow(2,i));
	//printf("%f\n", sum_x);

	for(i=0; i<=N_1; i++)
		sum_x2 = sum_x2 + (log10(pow(2,i))*log10(pow(2,i)));
	//printf("%f\n", sum_x2);

	for(i=0; i<=N_1; i++)
	{
		temp = p[i];
		if(temp>0 || temp<0)
			sum_y = sum_y + log10 (temp);
	}
	//printf("%f\n", sum_y);

	for(i=0; i<=N_1; i++)
	{
		temp = p[i];
		if(temp>0 ||temp<0)
			sum_xy = sum_xy + (log10(pow(2,i)) * log10(temp));
	}
	//printf("%f\n", sum_xy);

	m = ((N_1+1)*sum_xy - sum_x*sum_y) / ((N_1+1)*sum_x2 - pow(sum_x, 2));

	printf("Slope = %f\t H = 1-(slope)/2 = %f\n", m, 1-(m/2));

	return m;
}