iird.c

这是数字信号处理方面的一些源码

 (b1+1)->real = 1.0; (b1+1)->imag = 0.0;
 if(*n!=1) {
    for(i=1; i < *n;i++) {
	for(k=0;k<i;k++)   {
	  cs1=(b1+k)->real-(b1+k+1)->real * (b+i)->real;
	  cs2=(b1+k)->imag-(b1+k+1)->real * (b+i)->imag;
	  (b2+k+1)->real = cs1+ (b1+k+1)->imag * (b+i)->imag;
	  (b2+k+1)->imag = cs2- (b1+k+1)->imag * (b+i)->real;
	}
	b2->real = -(b1->real*(b+i)->real - b1->imag*(b+i)->imag);
	b2->imag = -(b1->real*(b+i)->imag + b1->imag*(b+i)->real);
	(b2+i+1)->real = (b1+i)->real;
	(b2+i+1)->imag = (b1+i)->imag;
	for(k=0;k<=i+1;k++)  {
	  (b1+k)->real = (b2+k)->real;
	  (b1+k)->imag = (b2+k)->imag;
	  (b2+k)->real = 0.0;
	  (b2+k)->imag = 0.0;
	}
     }
  }
  for(i=0;i <= *n;i++)
     h[i] = (b1+i)->real / a;

  free(b);
  free(b1);
  free(b2);
  return h;
}
/************************************************************************
bsf - 有理分式变换子程序
      有理分式格式为 H(s)=1/(a0+a1s+a2s^2+...+ans^n), 用 s=F1(Z)/F2(Z)
      来代换 s, 得到输出分式 H(Z);

输入参数:
      double *c     : H(s)的系数指针;
      int ni        : H(s)的阶数;
      double *f1,*f2: F1(Z)的系数和F2(Z)的系数;
      int nf        : S=F1(Z)/F2(Z) 的阶数;
输出参数:
      struct ptr (*a,*b): 输出H(Z)的分子和分母系数, 均为升幂排列;
      int *no       : H(Z)的阶数;
*************************************************************************/
struct rptr* bsf(double *c,int ni,double *f1,double *f2,int nf,struct rptr* ptr,int *no)
{
 int i,j,ii,nd,ne,ng;
 double *d=NULL,*e=NULL,*g=NULL;
 ptr->a=pnpe(f2,nf,ni,ptr->a,no);      /* calculate the numerator coefficients */
 ptr->b=(double *)calloc(*no+2,sizeof(double));
 if(ptr->b==NULL) {
       printf("\nNot enough memory to allocate!");
       exit(0);
    }
 for(i=0;i<= ni;i++) {     /* calculate the denominator coefficients */
    d=pnpe(f1,nf,i,d,&nd);
    e=pnpe(f2,nf,ni-i,e,&ne);
    g=ypmp(d,nd,e,ne,g,&ng);
    for(j=0;j<=ng;j++)
       ptr->b[j] += c[i]*g[j];
    free(d);
    free(e);
    free(g);
 }
 return ptr;
}

/************************************************************************
pnpe - 多项式乘方展开函数, 将 (a0+a1x+a2x^2+...+amx^m)^n 展开为
       b0+b1x+...+bmnx^mn, 输入和结果均为升幂顺序.

输入参数:
       double *a : 输入多项式系数指针;
       int m,n   : 输入多项式的乘方阶数;
输出参数:
       double *b : 输出多项式系数指针;
       int *mn   : 输出多项式阶数;
*************************************************************************/
double *pnpe(double *a,int m,int n,double *b,int *mn)
{
 int i,j,k,nk;
 double *c;
 *mn=m*n;
 c=(double *)calloc(*mn+1,sizeof(double));
 b=(double *)calloc(*mn+1,sizeof(double));
 if(b==NULL) {
       printf("\nNot enough memory to allocate!");
       exit(0);
 }
 if (n==0) {*b=1.00; free(c); return b;}
 else {
    for(i=0;i<=m;i++) b[i]=a[i];
    if(n==1) {free(c); return b;}
    else {
       nk=m+1;
       for(i=1;i<n;i++) {
	  for(j=0;j<=m;j++)
	     for(k=0;k<nk;k++) c[k+j] += a[j]*b[k];
	  nk += m;
	  for(k=0;k<nk;k++) {b[k]=c[k]; c[k]=0;}
       }
    }
    free(c);
 }
 return b;
}
/************************************************************************
ypmp - 多项式相乘函数, 将(a0+a1x+a2x^2+...+amx^m)(b0+b1x+...+bnx^n)
       展开为 c0+c1x+...+c(m+n)x^(m+n), 输入和输出均为升幂顺序.

输入参数:
       double *a : 输入多项式 a 系数指针;
       int m     : 输入多项式 a 阶数;
       double *b : 输入多项式 b 系数指针;
       int n     : 输入多项式 b 阶数;
输出参数:
       double *c : 输出多项式系数指针;
       int *mn   : 输出多项式阶数;
*************************************************************************/
double *ypmp(double *a,int m,double *b,int n,double *c,int *mn)
{
int i,j;
*mn=m+n;
c=(double *)calloc(*mn+1,sizeof(double));
if(c==NULL) {
       printf("\nNot enough memory to allocate!");
       exit(0);
 }
for(i=0;i<=m;i++)
   for(j=0;j<=n;j++)
      c[i+j] += a[i]*b[j];
return c;
}
/************************************************************************
draw_image - 将输入数据的幅度画出图形。该函数可自动调整显示的比例, 使图形
	     充满整个屏幕。

输入参数: double *x    -   输入数据序列的指针;
	  int m        -   输入数据序列的长度;
	  char *title1 -   显示图形的上标题字符串指针;
	  char *xdis1  -   X 坐标左边显示标题字符串指针.
	  char *title2 -   显示图形的左标题字符串指针.
	  char *xdis2  -   X 坐标右边显示标题字符串指针.
	  int dis_type -   显示类型, 0:连线 1:直线.
输出参数: 无
*************************************************************************/
void draw_image(double *x,int m,char *title1,char *title2,
		char *xdis1,char *xdis2,int dis_type)
{
 int gdriver=DETECT, gmode,errorcode;
 int i,scx,scy,y0,signa,signb;
 int style, userpat;
 int start_x=40,start_y=40,end_x=10,end_y=60;
 long tlen;
 double ys,xs,ym;
 char dis[40];
 /*initializes the graphics mode */
 initgraph(&gdriver,&gmode,"");
 errorcode=graphresult();
 if (errorcode != grOk) {
    printf("Graphics error: %s\n",grapherrormsg(errorcode));
    printf("Press any key to halt!\n");
    getch();
    exit(1);
 }
 scx=getmaxx();
 scy=getmaxy();
 ym=1.e-90;
 signa=0;
 signb=0;

 for(i=0;i<m;i++) {
    if ((*(x+i)>0)&&(*(x+i)>ym))  ym = *(x+i);
    if ((*(x+i)<0)&&(- *(x+i)>ym))  ym = - *(x+i);
 }
 for(i=0;i<m;i++)  {
    if (*(x+i)>fabs(ym/20)) signa=1;
    if (*(x+i)<-fabs(ym/20)) signb=1;
 }
 if ((signa==1)&&(signb==1)) ys=(double)((scy - start_y - end_y)>>1)/ym;
 else ys=(double)((scy - start_y - end_y)/ym);
 xs=(double)(scx - start_x - end_x)/m;
 y0=((scy - start_y - end_y)>>1)+start_y;

 /* draw the frame */

 setcolor(LIGHTGREEN);
 rectangle(start_x-1,start_y-20,scx-end_x+1,scy-end_y+20);

 setcolor(DARKGRAY);
 /* select the line style */
 style=DASHED_LINE;
 userpat = 1;
 setlinestyle(style, userpat, 1);
 /* a user defined line pattern */
 /* binary: "0000000000000001"  */
 for(i=0;i<=10;i++)
    line(start_x,start_y+(scy-start_y-end_y)*i/10,scx-end_x,start_y+(scy-start_y-end_y)*i/10);
 for(i=0;i<=10;i++)
    line(start_x+(scx-start_x-end_x)*i/10,start_y,start_x+(scx-start_x-end_x)*i/10,scy-end_y);
 setcolor(GREEN);
 style=SOLID_LINE;
 userpat = 1;
 setlinestyle(style, userpat, 1);
 rectangle(start_x,start_y,scx-end_x,scy-end_y);
 setcolor(YELLOW);
 for(i=0;i<=10;i++)
    line(start_x,start_y+(scy-start_y-end_y)*i/10,start_x+5,start_y+(scy-start_y-end_y)*i/10);
 for(i=0;i<=10;i++)
    line(start_x+(scx-start_x-end_x)*i/10,scy-end_y+15,start_x+(scx-start_x-end_x)*i/10,scy-end_y+20);
 settextstyle(DEFAULT_FONT,HORIZ_DIR,1);
 setcolor(YELLOW);
 if((signa==1)&&(signb==0)) {
    strcpy(dis,"0");
    outtextxy(start_x+2,scy-end_y+4,dis);
    gcvt(ym,5,dis);
    outtextxy(start_x+1,start_y-10,dis);
    outtextxy(start_x-10,scy-end_y+24,xdis1);
    outtextxy(scx-2-strlen(xdis2)*8,scy-end_y+24,xdis2);
 }
 else if((signb==1)&&(signa==0)) {
    strcpy(dis,"0");
    outtextxy(start_x+2,start_y-10,dis);
    gcvt(ym,5,dis);
    outtextxy(start_x+2,scy-end_y+4,"-");
    outtextxy(start_x+10,scy-end_y+4,dis);
    outtextxy(start_x-10,scy-end_y+24,xdis1);
    outtextxy(scx-2-strlen(xdis2)*8,scy-end_y+24,xdis2);
 }
 else {
    line(start_x,y0,scx-end_x,y0);
    strcpy(dis,"0");
    outtextxy(start_x-10,y0,dis);
    gcvt(ym,5,dis);
    outtextxy(start_x+2,start_y-10,dis);
    outtextxy(start_x+2,scy-end_y+4,"-");
    outtextxy(start_x+10,scy-end_y+4,dis);
    outtextxy(start_x-10,scy-end_y+24,xdis1);
    outtextxy(scx-2-strlen(xdis2)*8,scy-end_y+24,xdis2);
 }
 strcpy(dis,"Press any key to continue...");
 setcolor(LIGHTRED);
 outtextxy((scx-28*8)>>1,scy-16,dis);

 settextstyle(DEFAULT_FONT,HORIZ_DIR,2);
 tlen=strlen(title1);
 if ((tlen<<4)<scx) {
    setcolor(LIGHTGREEN);
    outtextxy((start_x+scx-end_x-(tlen<<4))>>1,start_y-40,title1);
 }

 settextstyle(DEFAULT_FONT,VERT_DIR,1);
 tlen=strlen(title2);
 if ((tlen<<4)<scy) {
    setcolor(LIGHTGREEN);
    outtextxy(start_x-20,(scy-end_y-(tlen<<3))>>1,title2);
 }
 /*draw the amplitude image*/
 setcolor(WHITE);
 if((signa==1)&&(signb==0)) y0=scy-end_y;
 else if((signb==1)&&(signa==0)) y0=start_y;
 if (dis_type == 0) {
    for(i=0;i<m-1;i++)
      line(xs*i+start_x,y0-*(x+i)*ys,xs*(i+1)+start_x,y0-*(x+i+1)*ys);
 }
 else if (dis_type == 1) {
    for(i=0;i<=m;i++)
      line(xs*i+start_x,y0-*(x+i)*ys,xs*i+start_x,y0);
 }
 getch();
 closegraph();
}