/*--------------------------------------------------------------*/
/*								*/
/*	Program for:						*/
/*	(1) computing the inverse z-transform via the power	*/
/*	series or partial fraction expansion method.		*/
/*								*/
/*	(2) converting a transfer function, H(z), in cascade	*/
/*	form to an equivalent transfer function in paralell,	*/
/*	via partial fraction expansion. The basic building 	*/
/* 	blocks is the second order biquad.			*/
/*								*/
/*								*/
/*	input file:	casacde.dat				*/
/*	output file:	xdata.dat				*/
/*								*/
/*	Manny Ifeachor, June, 1992				*/
/*--------------------------------------------------------------*/

#include "dsp1.h"	/* common structures +include files	*/
#include "dsp21.h"	/* global function + variable declarations*/

/*--------------------------------------------------------------*/
main()
{
extern double A[size], B[size], Ni[size], Di[size];
extern int M,N,N1,M1,nstage,iopt,iir;
extern long npt;
extern FILE *in, *out, *fopen();

	iir=0;
	read__coeffs();	/* go read the system coeffs */
	poly__product1(Di,B,nstage);	/*form B(z)*/
	poly__product1(Ni,A,nstage);	/*form A(z)*/
	M1=M;
	N1=N;
	if(B[M]==0) 		/*system of odd order */
	{
		M1=M-1;
		N1=M1;
	}
	printf("select desired operation\n");
	printf("0	for power series method of IZT\n");
	printf("1	for partial fraction coeffs estimation\n");
	printf("2	for cascade to parallel conversion\n");
	scanf("%d", &iopt);
	switch(iopt)
	{
		case 0:
			printf("enter number of data points required\n");
			scanf("%ld",&npt);
			power__series();
			izt__output();
			break;
		case 1:
			pfraction();
			print__pfcoeffs();
			break;
		case 2:
			pfraction();		/*compute PF coeffs*/
			
			/*compute paralell coeffs for H(z)*/
			
			cascade__parallel();
			print__pfcoeffs();	/*print PF coeffs*/
			
			/* print parallel coeffs for H(z)*/
			
			printpar();
			break;
		default:	break;
		}
		exit(0);
	}

/*----------------------------------------------------------------------------*/

/* function to evaluate the partial fraction coefficients		*/

void pfraction()
{
	int i;
	complex nz[10], dz[10];
	extern double A[size], B[size];
	extern complex pk[10], ck[10];
	extern double B0;
	extern int M,N1,M1,nstage;

	B0=A[M1]/B[M1];			 /*compute constant coeffs */
		
		polar__pole();		/*find positions of poles p1,p2,p3*/
		for(i=1;i<=M1;++i)
		{
						/* evaluate N(z) at p1 */
			nz[i]= poly__polar(A,pk[i],M1);
			 			/* evaluate N(z) at p1 */
			dz[i]=Dz(pk[i],i);	/* evaluate D'(z) at p1 */
			ck[i]=cdiv1(nz[i],dz[i]);  /*obtain c1=N(z)/D'(z), z=p1 */
		}
}
/*------------------------------------------------------------------------*/

/* function to evaluate the denominator polynomial D'(z) */

complex Dz(complex px, int npfc)
{
	int i,j;
	complex dk[10],dz;
	extern complex p[20];

	pole__array(npfc);
				/*evaluate Dk(z)=z(z-p[2]...(z-p[M]),z=px */
	if(M1==2)
	{
		dz=Dzz(px,px,p[1]);
	}
	if(M1>2)
	{
		j=1;
		dk[1]=Dzz(px,px,p[1]);
		for(i=2;i<(M1-1);++i)
		{
			dk[i]=Dzz(dk[i-1],px,p[i]);
			j=i;
		}
		dz=Dzz(dk[j],px,p[j+1]);
	}
	return(dz);
}
/*-------------------------------------------------------------------------*/

/* function to evaluate factors of the form z(z-p) */

complex Dzz(complex pa, complex pb, complex pc)
{
	complex t1,t2,dz;
	t1=cmul1(pa,pb);
	t2=cmul1(pa,pc);
	dz=cadd1(t1,t2,0);
	return(dz);
}

/*-------------------------------------------------------------------------*/

/*function to reorder the poles */

void pole__array(npfc)
{
	int i;

	switch(npfc)
	{
		case 1:
			for (i=1;i<=9;++i)
			{
				p[i]=pk[i+1];
			}
			break;
		case 2:
			p[1]=pk[1];
			for(i=2;i<=9;++i)
			{
				p[i]=pk[i+1];
			}
			break;
		case 3:
			p[1]=pk[1];
			p[2]=pk[2];
			for(i=3;i<=9;++i)
			{
				p[i]=pk[i+1];
			}
			break;
		case 4:
			p[1]=pk[1];
			p[2]=pk[2];
			p[3]=pk[3];
			for(i=4;i<=9;++i)
			{
				p[i]=pk[i+1];
			}
			break;
		case 5:
			for(i=1;i<=4;++i)
			{
				p[i]=pk[i];
			}
			for(i=5;i<=9;++i)
			{
				p[i]=pk[i+1];
			}
			break;
		case 6:
			for(i=1;i<=5;++i)
			{
				p[i]=pk[i];
			}
			for(i=6;i<=9;++i)
			{
				p[i]=pk[i+1];
			}
			break;
		 case 7:
			 for(i=1;i<=6;++i)
			 {
				 p[i]=pk[i];
			 }
			 for(i=7;i<=9;++i)
			 {
				 p[i]=pk[i+1];
			 }
			 break;
		case 8:
			for(i=1;i<=7;++i)
			{
				p[i]=pk[i];
			}
			p[8]=pk[9];
			p[9]=pk[10];
			break;
		case 9:
			for(i=1;i<=8;++i)
			{
				p[i]=pk[i];
			}
			p[9]=pk[10];
			break;
		case10:
			for(i=1;i<=9;++i)
			{
				p[i]=pk[i];
			}
			break;
	}
}

/*-------------------------------------------------------------------------*/

/* function to print the partial fraction coefficients */

void print__pfcoeffs()
{
	int k;
	long i;
	extern int M,M1;
	extern complex ck[10],pk[10];
	extern double B0;
	printf("poles of the z-transform\n");
	printf("\n");
	printf("pk \treal \t\timag \t\tmag \t\tphase\n");
	for(k=1;k<=M1;++k)
	{
		pk[k].angle=pk[k].angle*180/pi;
		i=(long)(pk[k].angle/360);
		pk[k].angle=pk[k].angle-i*360;
		if(pk[k].angle<-180)
			pk[k].angle=pk[k].angle+360;
		if(pk[k].angle>180)
			pk[k].angle=pk[k].angle-360;
		printf("%d \t%f \t%f \t%f \t%f\n",k,pk[k].real,pk[k].imag,pk[k].modulus,pk[k].angle);
	}
	printf("\n");
	printf("\n");
	printf("partial fraction coeffs\n");
	printf("\n");
	printf("B0=%f \n",B0);
	printf("\n");
	printf("Ck \treal \t\timag \t\tmag \t\tphase\n");
	for(k=1;k<=M1;++k)
	{
		ck[k].angle=ck[k].angle*180/pi;
		i=(long)(ck[k].angle/360);
		ck[k].angle=ck[k].angle-i*360;
		if(ck[k].angle<-180)
			ck[k].angle=ck[k].angle+360;
		if(ck[k].angle>180)
			ck[k].angle=ck[k].angle-360;
		printf("%d \t%f \t%f \t%f \t%f\n",k,ck[k].real,ck[k].imag,ck[k].modulus,ck[k].angle);
	}
	printf("press enter to continue \n");
	getchar();
}
/*--------------------------------------------------------------------------*/

/* Function to compute coefficients for parallel realization */

void cascade__parallel()
{
	int i;
	extern double ak[size], bk[size];
	extern complex pk[10], ck[10];
	extern double B0;
	extern int M;

		for(i=0;i< M;++i)
		{
			ncoeff(i);
			
			/* compute coeffs for 2nd order sections */
		}
}
/*---------------------------------------------------------------------------*/
void printpar()
{
	int i,j;
	extern double ak[size];
	extern int M;
	
	printf("stage \tNi(z)\n");
	printf("\n");
	for(i=0;i<=M/2;++i)
	{
		j=2*i;
		printf("%d \t%f \t%f\n",i,ak[j],ak[j+1]);
	}
	j=0;
	printf("\n");
	printf("\n");
	printf("stage \tDi(z)\n");
	printf("\n");
	for(i=0;i<=M/2;++i)
	{
		printf("%d \t%f \t%f \t%f\n",i,Di[j],Di[j+1],Di[j+2]);
	j=j+3;
	}

	printf("press enter to continue \n");
	getchar();
}
/*-------------------------------------------------------------------------*/

/* function to compute numerator coeffs of second order sections from partial fraction coeffs */

void ncoeff(int i)
{
	int j;
	complex temp1,temp2,temp3;
	extern double ak[size];
	extern complex pk[10], ck[10];

	j=2*i;
	temp1=cadd1(ck[j+1],ck[j+2],1);
	ak[j]=temp1.real;
	temp1=cmul1(ck[j+1],pk[j+2]);
	temp2=cmul1(ck[j+2],pk[j+1]);
	temp3=cadd1(temp1,temp2,1);
	ak[j+1]=-temp3.real;
}
/*--------------------------------------------------------------------------*/

void izt__output()
{
	long i;
	extern long npt;
	extern double h[size];

	if((out=fopen("xdata.dat","w"))==NULL)
	{
		printf("cannot open file xdata.dat\n");
		exit(1);
	}
	for(i=0;i<npt;++i)
	{
		fprintf(out,"%15e\n",h[i]);
	}
	fclose(out);
}