[人工智能]C语言实现DSP篇——线性卷积、快速傅里叶变换等

线性卷积：

/** Program for convolve **/

#include <math.h>

int N1,N2;		/*输入数组长度*/
int n;			/*输出数组长度*/
int m,i,k;
float x[20];
float h[20];
float y[20];	/*输出数组*/

main()
{	
   N1=10;                           /* x 长度*/
   N2=10;                           /* h 长度*/
   n=N1+N2-1;                      /* 输出 y 的长度*/  

   for(i=0;i<20;i++)               /* 初始化数组 */
   {
    
     x[i]=0;
     h[i]=0;
     y[i]=0;
    
   }

   for(i=0;i<n;i++)		/* 给x数组赋值 */              
   {
      if(i<N1)
      {
         x[i]=i;
  	  }
      
      else
      {
         x[i]=0;
      }
   }

   for(i=0;i<=n;i++)	/* 给h数组赋值 */
   {
      if(i<N2)
      {
         h[i]=1;
      }
      
      else
      {
         h[i]=0;
      }
   }

   for(i=0;i<n;i++)  	/* 计算卷积 */
   {
   	  for (k=0;k<=i;k++)
   	  	y[i]=y[i]+h[k]*x[i-k];    /* */      
   }
}

N点FFT：

/* Program for FFT */

#include <math.h>
#include <stdio.h>


const float pi=3.1415926;
int N;							/* FFT点数 */
float x_re[300], x_im[300];	 	/* 输入信号序列 */
float y_re[300], y_im[300];		/* 输出频谱序列 */
float w_re, w_im;				/* 蝶形因子 */
int m;							/* 蝶形运算的级数,即Log2(N) */
float t_re, t_im, v_re, v_im;  	/* 临时变量 */
int j,i,k,f,n;   				
int a, b, c;


main()
{ 	
	N=64;

	/* 初始化数据空间 */
	for(i=0; i<300; i++)
	{
		x_re[i]=0;
		x_im[i]=0;
	}
	
	/* 设定输入信号序列为单边指数函数 */
	/* 考虑到实际衰减很快，因此可以忽略后面大半部分数值 */
	for(i=0;i<=N;i++)
	{
		 x_re[i]=exp(-i);
		 x_im[i]=0;		
	}
	
	/* 复制到输出数组 */
	for(i=0; i<300; i++)
	{
		y_re[i]=x_re[i];
		y_im[i]=x_im[i];
	}

	/* 用雷德算法对输入信号序列进行倒序重排 */
  	j=0;
  	for(i=0;i<N;i++)
  	{	
  		if(i<j)
  		{ 
  		  t_re=y_re[j]; 
  		  t_im=y_im[j]; 
  		  y_re[j]=y_re[i]; 
  		  y_im[j]=y_im[i];
  		  y_re[i]=t_re;
		  y_im[i]=t_im;
  		}
  		k=N/2;
  		while((k<=j)&(k>0))
  		{ 
  		  j=j-k; 
  		  k=k/2;
  		}
  		j=j+k;
    }
    
    /* 计算蝶形运算的级数log2(N) */
	f=N;
	for(m=1; (f=f/2)!=1; m++);
	
    
    /*** FFT ***/
    for(n=1; n<=m; n++)
    {   
    	a=1;				/* a=2的n次方 */
    	for(i=0;i<n;i++)
    		a=a*2;
    	
    	b=a/2;
    	
    	v_re=1.0; 			/* 蝶形因子 */
    	v_im=0.0;
    	w_re=cos(pi/b);		
    	w_im=-sin(pi/b);
    	
    	for(j=0;j<b;j++)	/* 蝶形运算 */
    	{
    		for(i=j;i<N;i=i+a)
    		{
    			c=i+b;
   
       			t_re=y_re[c]*v_re-y_im[c]*v_im;
    			t_im=y_re[c]*v_im+y_im[c]*v_re;
    			
    			y_re[c]=y_re[i]-t_re;
    			y_im[c]=y_im[i]-t_im;
    			
    			y_re[i]=y_re[i]+t_re;
    			y_im[i]=y_im[i]+t_im;
    		}
    		
    		t_re=v_re*w_re-v_im*w_im;
    		t_im=v_re*w_im+v_im*w_re;
    		
    		v_re=t_re;
    		v_im=t_im;
    	}
    }
    
}

相关估计：

/** (Ex5_2.c)  of (Project Ex5_2 )**/
/** Program for Correlation Algorithm **/

#include "math.h"

float x[50],y[50];   /* 两个输入数组 */
float r[50];         /* 输出的相关函数估计值 */
int mode;	         /* 选择参数：0-无偏相关估计;1-有偏相关估计*/
int m;      		 /* 输出数组长度 */
int n;               /* 参与估计的输入数组长度 */
int k;				 /* 相关系数的阶数 */
int i,j;			 /* 将要用到的循环变量 */
float sum;           /* 将要用到的变量 */


main()
{  /* 参数初始化 */  
   mode=0; 
   m=10;                        
   n=40;
   k=0;
   sum=0;

   for(i=0;i<50;i++)
   {
      x[i]=0;
      y[i]=0;
      r[i]=0;   
   }

   for(i=0;i<n;i++)
   {
      x[i]=1;
	  y[i]=i;
   }

   for(k=0;k<m;k++)
   {
      sum=0;

      for(j=0;j<n-k;j++)  /* 对x[j+k]*y[j]求和，j从1到n-k-1 */
      {
         sum=sum+x[j+k]*y[j];    
      }

      if(mode==0)        /* 无偏估计 */
      {
         r[k]=sum/(float)(n-k);
      }

      else 				/* 有偏估计 */
      {
         r[k]=sum/(float)n;
      }
   }    
}

离散余弦变换：

/* Program for DCT */

#include <math.h>
#include <stdio.h>


const float pi=3.1415926;
int N;							/* FFT点数 */
int N2;							/* N2=2*N */
float x_re[300], x_im[300];	 	/* 输入信号序列 */
float y_re[300], y_im[300];		/* 输出频谱序列 */
float ck[300];					/* 离散余弦变换系数 */
float w_re, w_im;				/* 蝶形因子 */
int m;							/* 蝶形运算的级数,即Log2(N) */
float t_re, t_im, v_re, v_im;  	/* 临时变量 */
int j,i,k,f,n;   				
int a, b, c;


main()
{ 	
	N=8;
	N2=2*N;

	/* 初始化数据空间 */
	for(i=0; i<300; i++)
	{
		x_re[i]=0;
		x_im[i]=0;
		ck[i]=0;
	}
	
	/* 设定输入信号序列为单边指数函数 */
	/* 考虑到实际衰减很快，因此可以忽略后面大半部分数值 */
	for(i=0;i<N;i++)
	{
		 x_re[i]=exp(-i);
		 x_im[i]=0;		
	}
	
	/* 对信号序列进行延拓 */
	for(i=N;i<N2;i++);
	
	/* 复制到输出数组 */
	for(i=0; i<300; i++)
	{
		y_re[i]=x_re[i];
		y_im[i]=x_im[i];
	}

	/* 用雷德算法对输入信号序列进行倒序重排 */
  	j=0;
  	for(i=0;i<N2;i++)
  	{	
  		if(i<j)
  		{ 
  		  t_re=y_re[j]; 
  		  t_im=y_im[j]; 
  		  y_re[j]=y_re[i]; 
  		  y_im[j]=y_im[i];
  		  y_re[i]=t_re;
		  y_im[i]=t_im;
  		}
  		k=N2/2;
  		while((k<=j)&(k>0))
  		{ 
  		  j=j-k; 
  		  k=k/2;
  		}
  		j=j+k;
    }
    
    /* 计算蝶形运算的级数log2(N) */
	f=N2;
	for(m=1; (f=f/2)!=1; m++);
	
    
    /*** FFT ***/
    for(n=1; n<=m; n++)
    {   
    	a=1;				/* a=2的n次方 */
    	for(i=0;i<n;i++)
    		a=a*2;
    	
    	b=a/2;
    	
    	v_re=1.0; 			/* 蝶形因子 */
    	v_im=0.0;
    	w_re=cos(pi/b);		
    	w_im=-sin(pi/b);
    	
    	for(j=0;j<b;j++)	/* 蝶形运算 */
    	{
    		for(i=j;i<N2;i=i+a)
    		{
    			c=i+b;
   
       			t_re=y_re[c]*v_re-y_im[c]*v_im;
    			t_im=y_re[c]*v_im+y_im[c]*v_re;
    			
    			y_re[c]=y_re[i]-t_re;
    			y_im[c]=y_im[i]-t_im;
    			
    			y_re[i]=y_re[i]+t_re;
    			y_im[i]=y_im[i]+t_im;
    		}
    		
    		t_re=v_re*w_re-v_im*w_im;
    		t_im=v_re*w_im+v_im*w_re;
    		
    		v_re=t_re;
    		v_im=t_im;
    	}
    }
    
    /* 提取实部作为离散余弦变换的系数 */
    
   	for(i=0;i<N;i++)
   	{
   		ck[0]=ck[0]+x_re[i];
   	}
   	ck[0]=ck[0]/sqrt((float)N);
   	
    for(i=1;i<N;i++)
    {
    	ck[i]=cos((i*pi)/N2)*y_re[i]-sin(-(i*pi)/N2)*y_im[i];
    	ck[i]=sqrt(2.0/N)*ck[i];
    }
}

IIR滤波器：

/* Program for IIR */

#include <math.h>
#include <stdlib.h>

#define pi 3.1415925

float fp,fr,fs;		/* 通带截止频率、阻带截止频率、抽样频率 */
float ap, ar;		/* 容限(dB) */
float ptr_a[50],ptr_b[50]; /* 系统函数H(z)的分母系数和分子系数 */
float hwdb[50];				/* 幅频响应值 */
float b_real[50],b_imag[50];
float b1_real[50],b1_imag[50],b2_real[50],b2_imag[50];
float d[50],e[50],g[50];
float f1[2],f2[2],h[50];


/* 函数定义 */

void bcg(float ap,float as,float wp,float ws,int *n,float *h)                     /*求H(s)分母系数*/
{
   int i,k;
   float a,p,wc,cs1,cs2;
   float c;

   c=(pow(10.0,0.1*as)-1.0)/(pow(10.0,0.1*ap)-1.0);
   *n=(int)(fabs(log10(c)/log10(ws/wp)/2.0)+0.99999);                                     /*求N*/
   wc=wp;
   a=pow(wc,(double)(*n));

   for(i=0;i<*n;i++)                           /*求极点*/
   {
      p=pi*(0.5+(2.0*i+1.0)/(2.0*(*n)));
      b_real[i]=wc*cos(p);
      b_imag[i]=wc*sin(p);
   }     

   b1_real[0]=-(b_real[0]);
   b1_imag[0]=-(b_imag[0]);
   b1_real[1]=1.0;
   b1_imag[1]=0.0;

   if(*n!=1)
   {
      for(i=1;i<*n;i++)
      {
         for(k=0;k<i;k++)
         {
            cs1=b1_real[k]-b1_real[k+1]*b_real[i];
            cs2=b1_imag[k]-b1_real[k+1]*b_imag[i];
            b2_real[k+1]=cs1+b1_imag[k+1]*b_imag[i];
            b2_imag[k+1]=cs2-b1_imag[k+1]*b_real[i];
         }

         b2_real[0]=-(b1_real[0]*b_real[i]-b1_imag[0]*b_imag[i]);
         b2_imag[0]=-(b1_real[0]*b_imag[i]+b1_imag[0]*b_real[i]);
         b2_real[i+1]=b1_real[i];
         b2_imag[i+1]=b1_imag[i];

         for(k=0;k<=i+1;k++)
         {
            b1_real[k]=b2_real[k];
            b1_imag[k]=b2_imag[k];
            b2_real[k]=0.0;
            b2_imag[k]=0.0;
         }
      }
   }

   for(i=0;i<=*n;i++)
      h[i]=b1_real[i]/a;
}

void pnpe(float *a,int m,int n,float *b,int *mn)                               
{	
	/*多项式相乘的展开系数*/
   
   int i,j,k,nk;
   float c[50];

   *mn=m*n;

   for(i=0;i<*mn+1;i++)
   {
      c[i]=0;
      b[i]=0;
   }

   if(n==0)
      b[0]=1.00;
   else
   {
      for(i=0;i<=m;i++)
         b[i]=a[i];
      if(n==1)
      {
         for(i=0;i<=m;i++)
            b[i]=a[i];
      }
      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;
            }
         }
      }
   }
}



void ypmp(float *a,int m,float *b,int n,float *c,int *mn)
{
   int i,j;

   *mn=m+n;

   for(i=0;i<*mn+1;i++)
      c[i]=0;

   for(i=0;i<=m;i++)
   {
      for(j=0;j<=n;j++)
         c[i+j]+=a[i]*b[j];
   }
}


void bsf(float *c,int ni,float *f1,float *f2,int nf,float *ptr_a,float *ptr_b,int *no)
{
   
   int i,j,nd,ne,ng;

   /*计算分子系数*/
   pnpe(f2,nf,ni,ptr_a,no);                                            

   for(i=0;i<*no+2;i++)
      ptr_b[i]=0;

   /*计算分母系数*/
   for(i=0;i<=ni;i++)                                                           
   {
      pnpe(f1,nf,i,d,&nd);
      pnpe(f2,nf,ni-i,e,&ne);
      ypmp(d,nd,e,ne,g,&ng);

      for(j=0;j<=ng;j++)
         ptr_b[j]+=c[i]*g[j];
   }
}


//*** 主函数 ***//
main()
{
   int N,nf,ns,nz,i,k;
   float hwdb1_real,hwdb1_imag;	/* 系统函数的分子部分 */
   float hwdb2_real,hwdb2_imag; /* 系统函数的分母部分 */
   float wp,ws,jw,amp1,amp2;

   N=50;

   /* 初始化数据空间 */
   for(i=0;i<50;i++)
   {
      b_real[i]=0;
      b_imag[i]=0;
      b1_real[i]=0;
      b1_imag[i]=0;
      b2_real[i]=0;
      b2_imag[i]=0;
      d[i]=0;
      e[i]=0;
      g[i]=0;
      ptr_a[i]=0;
      ptr_b[i]=0;
      hwdb[i]=0;
   }

	/* 滤波器设计参量 */
   fp=100;	/* Hz */
   fr=300;
   fs=1000;
   ap=3;	/* db */
   ar=20;
   wp=tan(pi*fp/fs);	/* 双线性变换 */
   ws=tan(pi*fr/fs);
   *f1=-1.0;
   *(f1+1)=1.0;
   *f2=1.0;
   *(f2+1)=1.0;
   nf=1;
   
   /* 求系统函数 H(s) */
   bcg(ap,ar,wp,ws,&ns,h);            
   
   /* 由H(s)转换成H(z) */
   bsf(h,ns,f1,f2,nf,ptr_b,ptr_a,&nz);           

   /* 用exp(jw)代替z */
   for(k=0;k<N;k++)
    {
      jw=k*pi/N;
      hwdb1_real=0;  hwdb1_imag=0;
      hwdb2_real=0;  hwdb2_imag=0;

      for(i=0;i<=nz;i++)
      {
		/* 分子 */      
         hwdb1_real+=ptr_b[i]*cos(i*jw);
         hwdb1_imag+=ptr_b[i]*sin(i*jw);
        /* 分母 */
         hwdb2_real+=ptr_a[i]*cos(i*jw);
         hwdb2_imag+=ptr_a[i]*sin(i*jw);
      }
      
      /* 求幅频响应 */
      amp1=(pow(hwdb1_real,2)+pow(hwdb1_imag,2));
      amp2=(pow(hwdb2_real,2)+pow(hwdb2_imag,2));          
      hwdb[k]=10*log10(amp1/amp2);
      if(hwdb[k]<-200)  hwdb[k]=-200;        
   }
}

LMS自适应滤波：

/* Program for LMS */

#include <math.h>

#define pi 3.1415926

float x[500];		/* 输入信号序列 */
float y[500];		/* 输出信号序列 */
float d[500];		/* 期望输出序列 */
float e[500];		/* 输出误差序列, e(n)=d(n)-y(n) */
float w[50];		/* 横向滤波器系数 */
int L;				/* 实际信号长度 */
int N;				/* 横向滤波器级数 */
int i, j;
float u;			/* 梯度算法步长 */


main()
{ 
   	L=500;
   	N=5;
	u=0.00001;

   	/* 初始化数据空间 */
   	for(i=0;i<500;i++)
   	{ 
   		x[i]=0;
      	d[i]=0;
      	e[i]=0;
      	y[i]=0;
   	}                                                                   

   	/* 设定输入信号为正弦信号
	   期望输出设为比实际输入超前两个单位时间的信号 */
   	for(i=0;i<L;i++)
   	{ 
      	x[i]=(float)10*sin(pi*i/20);
      	d[i]=x[i-2];
   	}

   	/* 初始化滤波器系数 */
   	for(i=0;i<N;i++)
    {  	
    	w[i]=0;
    }

	/* 梯度算法 */
   	for(i=N;i<L;i++)
   	{        
      /*sum=0;*/
      	for(j=0;j<N;j++)
         	y[i]=y[i]+w[j]*x[i-j];
      /*y[i]=sum;*/
      	e[i]=d[i]-y[i];		/* 计算误差 */
      	
      	for(j=0;j<N;j++)	/* 根据误差修正滤波器系数 */
         	w[j]=w[j]+2*u*e[i]*x[i-j];
   }   
}

窗函数法设计FIR滤波器：

/** Program to design an FIR **/

#include "math.h"

float hd[51];	/* 理想的冲击响应系数序列 */
float h[51];	/* 加窗后的冲击响应系数序列 */
float w[51];	/* 窗函数的离散值序列 */
float db[300];  /* 幅度响应数组（对数值）�*/
float pi;		/* 圆周率pi */
float wc;		/* (截止频率/pi) */
int n;			/* 窗函数序列的长度，也是加窗后的响应序列的长度 */
int m;			/* 选择参数(1-5):选择不同的窗函数 */
int l;			/* l>>n */
float im,re;
float a,b,p,wf,d;
int k,i;

main()
{
   	  /* 参数初始化 */
      for(i=0;i<51;i++)
      {		
      	hd[i]=0;
      	h[i]=0;
      	w[i]=0;
      }

	  for(i=0;i<300;i++)
      {		
      	db[i]=0;
      }
            
      m=5;
      n=21;    /* n<=51, 为奇数 */
      wc=0.2;  /* wc=0.10-0.90 */
      l=300;   
      a=(n-1)/2;
      pi=4.0*atan(1.0);

      for(i=0;i<n;i++)
      {	/* 理想矩形函数（频域）的时域形式 */
         if(i==a)
            hd[i]=wc;
         else
         {
            b=i-a;
            hd[i]=sin(pi*b*wc)/(pi*b);
         }
      }

      switch(m)  /* 根据m的值选择不同的窗函数 */
      {
		/* m=1时用矩形窗 */
         case 1:    			
            for(i=0;i<n;i++)
            w[i]=1.0;  			
         break;
         
		/* m=2时用三角窗(巴特莱特窗) */
         case 2:   				
            for(i=0;i<n;i++)
            {
               if(i>=a)
               w[i]=2.0-2.0*i/(n-1);
               else
               w[i]=2.0*i/(n-1);
            }
         break;


		/* m=3时用汉宁窗 */
         case 3:   			
            for(i=0;i<n;i++)
              w[i]=0.5*(1.0-cos(2.0*pi*i/(n-1)));
            break;


		/* m=4时用汉明窗 */
         case 4:
            for(i=0;i<n;i++)
              w[i]=0.54-0.46*cos(2.0*pi*(float)i/(n-1));
            break;


		/* m=5时用布莱克曼窗 */
         case 5:
            for(i=0;i<n;i++)
              w[i]=0.42-0.5*cos(2.0*pi*i/(n-1))+0.08*cos(4.0*pi*i/(n-1));
         break;
      }


      for(i=0;i<n;i++)
         h[i]=hd[i]*w[i];  /* 对理想冲击响应加窗（时域相乘） */

      p=pi/l;

      
      /* 对加窗后的响应序列进行傅里叶变换，观察其幅频响应 */
      for(k=0;k<=l-1;k++)
      {
         wf=(pi*k)/l;
         re=0.0;
         im=0.0;
         
         for(i=0;i<n;i++)
         {
            re=re+h[i]*cos((float)i*wf);
            im=im+h[i]*sin((float)i*wf);
         }
         
         d=sqrt(pow(re,2)+pow(im,2)); /*求模(幅度)*/
         db[k]=20.0*log10(d);         /*转换为对数表示形式*/
      }
}