/*********************************************************************
*利用zhao等人論文原理的，使用遞歸方法的DCT程序                       *
*																	 *
**********************************************************************/
#include"stdio.h"
#include"stdlib.h"
#include"math.h"

#define PI 3.14159265358979

double * DCT(double * f,int N);
double sum(double * a,int k);
//   1 迭代變換函數(shù)
double * Recursion(double * f,int N)
{
	double * t;
    int i;

	t=DCT(f,N);
	t[0]=(t[0]*2/N/1.41421356237310);
	for(i=1;i<N;i++)
		t[i]=(t[i]*2/N);

	return t;
}



//2  迭代求A函數(shù)
double * DCT(double * f,int N)
{
    double *temp1;
	double *temp2;
	int i;

    
	if(N>2)
	{
		temp1=(double *)malloc(N*sizeof(double));
		temp2=(double *)malloc(N*sizeof(double));
		for(i=0;i<=(N/2-1);i++)
		{
			temp1[i]=f[i]+f[N-1-i];

			temp1[i+N/2]=f[i]-f[N-1-i];
			temp1[i+N/2]*=cos(PI*(2*i+1)/2/N);  
		}
		
		f=DCT(temp1,N/2);  //此時(shí)前半部分的temp1已是Pn之前的元素  而f指向temp1的首地址
		temp2=DCT(temp1+N/2,N/2);

		for(i=1;i<=(N/2-1);i++)
		{
			temp2[i]=temp2[i]+temp2[i];
		}

		for(i=0;i<=N/2-1;i++)
		{
		    f[i+N/2]=sum(temp2,i);
		}                   //此時(shí)f記錄的是Pn之前的列向量
    
        for(i=0;i<N;i++)
		{
			if(i%2==0)    temp1[i]=f[i/2];
			else          temp1[i]=f[((i-1+N)/2)];
		}

	}
	else    //N=2時(shí)
	{
         temp1=(double *)malloc(2*sizeof(double));
		 temp1[0]=f[0]+f[1];
		 temp1[1]=cos(PI/4)*(f[0]-f[1]);
		 
	}
	return temp1;

}

//3 處理Bn分解第一個(gè)矩陣的求和函數(shù) 
double sum(double * a,int k)
{
     double sum=0;
	 double sum1=0;
	 double sum2=0;
	 int i;

     if(k%2==0)
	 {
		    sum1=a[0];
			for(i=2;i<=k;i=i+2)
			{
				sum1+=a[i];
				sum2+=a[i-1];
			} 
			sum=sum1-sum2;
	  }
	 else   //k為奇數(shù)
	 {
		 for(i=0;i<=k-1;i=i+2)
		 {
			 sum1+=a[i];
			 sum2+=a[i+1];
		 }
		 sum=sum2-sum1;
	 }
	 return sum;
}