C语言实现数字滤波器

// 定义PI 参数

#define PI 3.1415926

// 包含文件

#include "stdio.h"
#include "math.h"
#include "stdlib.h"

#include "fft_sub.c"

// 主函数开始
void main()
{

// 定义函数、数组、变量

  void fft();
  float *xr,*xi;
  float *H;
  int i,np,nfft,k;
  float t,dt,df,f,z,fc1,fc2;
  FILE *fpar,*fp1,*fp2,*fp3,*fp4,*fp5;
  char fil1[80],fil2[80],fil3[80],fil4[80],fil5[80];

// 打开参数文件
  fpar=fopen("filter.par","r");
  fscanf(fpar,"%s",fil1);
  fscanf(fpar,"%s",fil2);
  fscanf(fpar,"%s",fil3);
  fscanf(fpar,"%s",fil4);
  fscanf(fpar,"%s",fil5);

  fscanf(fpar,"%d",&np);
  printf("np=%d\n",np);
  fscanf(fpar,"%f",&dt);
  printf("dt=%8.3f ms\n",dt);

// 注意时间毫秒与秒的转换
  dt=(float)(dt/1000.0);
  fscanf(fpar,"%f%f",&fc1,&fc2);
  printf("fc1=%5.1f  fc2=%5.1f\n",fc1,fc2);
  fclose(fpar);

  // （1）calculate fft point，输入数据点数可能不一定是2**k
  k=(int)(log(np)/log(2));
  if(np>pow(2,k))k=k+1;
  nfft=(int)pow(2,k);
  df=(float)(1.0/(nfft*dt));
  printf("nfft=%d  k=%d\n",nfft,k);
  printf("dt=%8.4f  df=%8.4f\n",dt,df);
  // allocate memory
  xr=(float *)calloc(nfft,sizeof(float));
  xi=(float *)calloc(nfft,sizeof(float));
  H=(float *)calloc(nfft,sizeof(float));

  // （2）产生滤波器，generate a filter
  fp1=fopen(fil1,"w");
  for(i=0;i<=nfft/2;i++)
 { f=i*df;
   if(f<=fc2 && f>=fc1)H[i]=1.0;
   else H[i]=0.0;
   fprintf(fp1,"%4d \t %10.4f \t %12.4f\n",i,f,H[i]);
 }

  // 注意滤波器的对称性
  for(i=nfft/2+1;i<nfft;i++)
 { f=i*df;
   H[i]=H[nfft-i];
   fprintf(fp1,"%4d \t %10.4f \t %12.4f\n",i,f,H[i]);
 }
  fclose(fp1);
  printf("The Filter is ok!\n");

  fp2=fopen(fil2,"r");
// 输入数据
  for(i=0;i<np;i++)
 { fscanf(fp2,"%f%f",&t,&z);
   xr[i]=z;
 }
  fclose(fp2);
  printf("The Data is ok!\n");

//（3） fill zero，如果不满足2**k，要补零
  if(np<nfft)
 {for(i=np;i<nfft;i++)
  xr[i]=0;
 }
  for(i=0;i<nfft;i++)
  xi[i]=0.0;
  //（4） calculate fft
  fft(xr,xi,k,1);

  // （5）output amplitude spectrum of the orignal signal
  fp3=fopen(fil3,"w");
  for(i=0;i<nfft;i++)
 { f=i*df;
   fprintf(fp3,"%4d \t %10.4f \t %12.4f\n",i,f,sqrt(xr[i]*xr[i]+xi[i]*xi[i])*dt);
 }
  fclose(fp3);

  // （6）filter the X(f),output Y(f)
  fp4=fopen(fil4,"w");
  for(i=0;i<nfft;i++)
 { xr[i]=xr[i]*H[i];
   xi[i]=xi[i]*H[i];
   f=i*df;
   fprintf(fp4,"%4d \t %10.4f \t %12.4f\n",i,f,sqrt(xr[i]*xr[i]+xi[i]*xi[i])*dt);
 }
  fclose(fp4);

  printf("The |y(f)| is ok!\n");

// （7）do IFFT, get y(t)
  fft(xr,xi,k,-1);

  fp5=fopen(fil5,"w");
  for(i=0;i<=np;i++)
 { t=i*dt;
   fprintf(fp5,"%10.4f \t %12.4f\n",t,xr[i]);
 }
  fclose(fp5);
  printf("The y(t) is ok!\n");

  printf("The Program is finished!\n");
  free(xr);free(xi);free(H);
}