/***********************************************************************
 利用雙線性變換法實現數字橢圓濾波器的設計
***********************************************************************/
#include    <math.h>
#include    <dos.h>
#include    <stdlib.h>
#include    <stdio.h>
#include    <string.h>
#include    <conio.h>
#include    <graphics.h>

#define    PI	 ((double)(4.0*atan(1.0)))
/* COMPLEX STRUCTURE */
typedef struct {
    double real, imag;
} COMPLEX;

struct rptr{
   double *a;
   double *b;
};

/***********************************************************************/
void ellipg(double ap,double as,double wp,int *n,int *m,double *h0,double *d0,
       double *a,double *b,double *c);
struct rptr* ebsf(double *a,double *b,double *c,double h0,double d0,int m,
       int ni,double *f1,double *f2,int fn,struct rptr* ptr,int *no);
struct rptr* ebss(double *a,double *b,double *c,double h0,double d0,int m,
       int ni,double *s1,double *s2,int fn,struct rptr* ptr,int *no);
double *pnpe(double *a,int m,int n,double *b,int *mn);
double *ypmp(double *a,int m,double *b,int n,double *c,int *mn);
void elowpass_input(double *wp,double *ap,double *ar,double *f1,double *f2,double *s1,double *s2,int *nf);
void ehighpass_input(double *wp,double *ap,double *ar,double *f1,double *f2,double *s1,double *s2,int *nf);
void ebandpass_input(double *wp,double *ap,double *ar,double *f1,double *f2,double *s1,double *s2,int *nf);
void draw_image(double *x,int m,char *title1,char *title2,
		char *xdis1,char *xdis2,int dis_type);
/***********************************************************************/
void main(void)
{
 double f1[4],f2[4],s1[4],s2[4],*hwdb,jjw,sign;
 double a[30],b[30],c[30],h0,d0;
 int N,m,nf,ns,nz,i,j,k,ftype,type;
 COMPLEX hwdb1,hwdb2;
 double wp,ws,ap,as,jw,jw1,amp1,amp2;
 char title[80],tmp[20];
 struct rptr *ptr=NULL;

 N=500;
 hwdb=(double*)calloc(N+2,sizeof(double));
 if(hwdb==NULL) {
       printf("\nNot enough memory to allocate!");
       exit(0);
 }
 printf("\n1.Lowpass  2.Highpass  3.Bandpass");
 printf("\nPlease select the filter type:");
 ftype=getch();
 putch(ftype);
 ftype=ftype&0x00ff;
 if((ftype!=0x31)&&(ftype!=0x32)&&(ftype!=0x33)) {
    do {
       sound(2000); delay(100); nosound();
       printf("\n1.Lowpass  2.Highpass  3.Bandpass");
       printf("\nPlease select the filter type:");
       ftype=getch();
       putch(ftype);
       ftype=ftype&0x00ff;
    }while((ftype!=0x31)&&(ftype!=0x32)&&(ftype!=0x33));
 }
 switch(ftype) {
    case 0x31: elowpass_input(&wp,&ap,&as,f1,f2,s1,s2,&nf);
	       break;
    case 0x32: ehighpass_input(&wp,&ap,&as,f1,f2,s1,s2,&nf);
	       break;
    case 0x33: ebandpass_input(&wp,&ap,&as,f1,f2,s1,s2,&nf);
	       break;
    default:   elowpass_input(&wp,&ap,&as,f1,f2,s1,s2,&nf);
 }

 ellipg(ap,as,wp,&ns,&m,&h0,&d0,a,b,c);

 printf("Elliptic low-pass filter order = %d",ns);
 printf("\nThe analog filter coefficients of Ha(s):");
 printf("\ni             Ai                  Bi                  Ci");
 for(i=1;i<=m;i++)
    printf("\n%1d   %16.8lf    %16.8lf    %16.8lf",i,a[i],b[i],c[i]);
 printf("\n    H0   = %20.8lf",h0);
 printf("\n    D0   = %20.8lf",d0);
 printf("\nPress any key to calculate the Frequency Response Curve...");
 getch();

 for(i=0;i<N;i++)  {
   jw=pow(5.0*i/N,2);
   jjw=1.0;
   for(j=1;j<=m;j++)  {
     jw1=pow(a[j]-jw,2)/(pow(c[j]-jw,2)+pow(b[j],2)*jw);
     jjw=jjw*jw1;
   }
   if(fmod(ns,2)==0)  jjw=jjw*pow(h0,2);
   else   jjw=jjw*pow(h0,2)/(pow(d0,2)+jw);
   if(jjw<1e-20) jjw=1e-20;
   hwdb[i]=10*log10(jjw);
   if(i%10==9) printf("*");
 }
 strcpy(title,"The Analog Normalized LP (C type) N=");
 itoa(ns,tmp,10);
 strcat(title,tmp);
 strcpy(tmp,"5(rad)");
 draw_image(hwdb,N,title,"The Attenuation (dB)","0",tmp,0);

 ptr=ebss(a,b,c,h0,d0,m,ns,s1,s2,nf,ptr,&nz);
 printf("\nThe analog filter coefficients of H(s):");
 printf("\n(a is numerator coefficient. b is denominator coefficient.)");
 for(i=0;i<=nz;i++){
    printf("\na[%2d]=%20.8lf      b[%2d]=%20.8lf",i,ptr->a[i],i,ptr->b[i]);
    if(i%21==20) getch();
 }
 printf("\n\nPress any key to calculate the filter response of H(s)...");
 getch();
 printf("\nWaitting for calculating...");

 for(k=0;k<N;k++)  {
    jw=2.0*PI*k/N;
    hwdb1.real=ptr->a[0]; hwdb1.imag=ptr->a[1]*jw;
    hwdb2.real=ptr->b[0]; hwdb2.imag=ptr->b[1]*jw;
    sign=1.0;
    for(i=1;i<=(nz-1)/2;i++)  {
       sign=-sign;
       hwdb1.real += sign*ptr->a[2*i]*pow(jw,2.0*i);
       hwdb2.real += sign*ptr->b[2*i]*pow(jw,2.0*i);
       hwdb1.imag += sign*ptr->a[2*i+1]*pow(jw,2.0*i+1);
       hwdb2.imag += sign*ptr->b[2*i+1]*pow(jw,2.0*i+1);
    }
    if(fmod(nz,2.0)==0.0) {
       sign=-sign;
       hwdb1.real += sign*ptr->a[2*i]*pow(jw,2.0*i);
       hwdb2.real += sign*ptr->b[2*i]*pow(jw,2.0*i);
    }
    amp1 = hwdb1.real*hwdb1.real+hwdb1.imag*hwdb1.imag;
    amp2 = hwdb2.real*hwdb2.real+hwdb2.imag*hwdb2.imag;
    if(amp2<1e-60) amp2=1e-60;
    hwdb[k]=amp1/amp2;
    if(hwdb[k]<=1e-20) hwdb[k]=1e-20;
    hwdb[k]=10*log10(hwdb[k]);
    if(k%10==9) printf("*");
 }
 strcpy(title,"Analog Transfer Property (C type) N=");
 itoa(ns,tmp,10);
 strcat(title,tmp);
 strcpy(tmp,"2