Make your own free website on Tripod.com

To see the Dsp1.h Header File
To see the dsp21.h Header File
To see the dsp21ex.h Header File
To see the LTI Library Source code
To see the makefile
To see the input for IZT using Power Series Method
To see the input for IZT for Power Series Expansion
To see the input for IZT for Cascade to Paralell Structure Conversion
To see the output for IZT using Power series Method
To see the output for IZT using Power Series Expansion
To see the output for IZT using Cascade to Paralell Structure Conversion
To see the main program
To go back to index
To go back to homepage
/*--------------------------------------------------------------*/
/*								*/
/*	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);
}