/*--------------------------------------------------------------*
*		Impulse Invariant Method			*
*								*
*	The analog transfer function must be frequency-scaled	*
*	(normalized frequency) before using program		*
*	30.10.92						*
*								*
*---------------------------------------------------------------*
*/

#include 
#include 
/*#include */

void dfilter();
double T;
double a0,a1,a2,b0,b1,b2;
double p1,p2,pr,pi;
double c1,c2,cr,ci;
float A0,A1,B0,B1,B2,temp;

main()
{

/* initialize coefs*/

	A0=0;
	A1=0;
	B0=1;
	B1=0;
	B2=0;
	a0=0;
	a1=0;
	b0=1;
	b1=0;
	b2=0;
	c1=0;
	c2=0;
	p1=0;
	p2=0;
	a2=0;

	printf("impulse invariant discrete filters \n");
	printf("\n");
	printf("enter s-plane coefficients \n");
	printf("enter denominator coeffs: B0,B1,B2 \n");
	scanf("%f%f%f",&B0,&B1,&B2);
	printf("enter numerator coeffs: A0,A1 \n");
	scanf("%f,%f",&A0,&A1);
	T=1;
	dfilter();
	printf("\n");
	printf("press enter to continue\n");
	getchar();
	exit(0);
}

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

void dfilter()
{
	temp=B1*B1 - 4*B0*B2;
	if(B2==0)
	{
		p1=-B0/B1;
		a0=A0/B1;
		b1=-exp(p1*T);
	}
	if(temp>0)
	{
		pr=-B1/(2*B2);
		pi=(pr*pr)-B0/B2;
		pi=sqrt(pi);
		p1=pr+pi;
		p2=pr-pi;
		c1=(A0+A1*p1)/((p1-p2)*B2);
		c2=A1/B2-c1;
		a0=c1+c2;
		a1=-(c1*exp(p2*T)+c2*exp(p1*T));
		b1=-exp(p1*T)-exp(p2*T);
		b2=exp((p1+p2)*T);
	}
	if(temp<0)
	{
		pr=-B1/(2*B2);
		pi=(pr*pr)-B0/B2;
		pi=sqrt(-pi);
		cr=A1/(B2*2);
		ci=-(A0+A1*pr)/(2*pi*B2);
		a0=2*cr;
		a1=-(cr*cos(pi*T)+ci*sin(pi*T))*2*exp(pr*T);
		b1=-2*exp(pr*T)*cos(pi*T);
		b2=exp(2*pr*T);
	}
	printf("discrete filter coeffs: \n");
	printf("a0 a1 a2: \t%f %f %f \n",a0,a1,a2);
	printf("b0 b1 b2: \t%f %f %f \n",b0,b1,b2);
}