几个数字信号处理算法程序

核心提示： 说明:其中COMP为复数类型/*** Name :DFT* Function :Disperse Fuliye Transformation* Params :N -- Total count of sampling points* X -- Input sequence* Return

说明:其中COMP为复数类型

/*****************************************
*
* Name　　 :DFT
*　 Function :Disperse Fuliye Transformation
*　 Params　 :N -- Total count of sampling points
*　　　　　　 X -- Input sequence
*　 Return　 :XN(k)=sum[x(n)*Pow(e,j2*Pi/N)]
*　　　　　　　　　 k,n:0..N-1
*
*
*****************************************/
void DFT(int N,COMP x[],COMP XK[])
{
　　double C = (2*pi)/N;
　　COMP t(0,0),ret(0,0);
　　for(int k=0;k < N;k++)
　　{
　　　　ret = COMP(0,0);
　　　　for(int i=0;i< N;i++)
　　　　{
　　　　　　t = COMP(cos(C*k*i),-sin(C*k*i));
　　　　　　ret += x[i]*t;
　　　　}
　　　　XK[k] = ret;
　　}
　　
}
/*****************************************
*
* Name　　 :FFT
*　 Function :Fast Fuliye Transformation
*　 Params　 :N -- Total count of sampling points
*　　　　　　 X -- Input sequence
*　 Return　 :XN(k)=sum[x(n)*Pow(e,j2*Pi/N)]
*　　　　　　　　　 k,n:0..N-1
*
*
*****************************************/
void FFT(int N,COMP X[],COMP XK[])
{
　　int j=0;
　　COMP U=0,W=0;
　　COMP* A = XK;
　　
　　//Adjust sequence
　　for(int i=0;i< N;i++)
　　{
　　　　if(i==0)
　　　　{
　　　　　　A[0] = X[0];
　　　　}
　　　　else
　　　　{
　　　　　　j=GetInverse(N,j);
　　　　　　A[i] = X[j];
　　　　}
　　}
　　
　　//确定级别数
　　for(int M=0;M< N;M++)
　　{
　　　　if((1<< M)==N)
　　　　　　break;
　　}
　　
　　
　　for(int L=1;L<=M;L++)//1-M级依次确定
　　{
　　　　int LE = (int)pow(2,L);//间隔
　　　　int LE1 = LE/2;//W级数,如W0,W1,W2...
　　　　
　　　　W=COMP(cos(pi/LE1),-sin(pi/LE1));
　　　　U=COMP(1,0);
　　　　for(j=0;j< LE1;j++)//
　　　　{
　　　　　　i=j;
　　　　　　while(i< N)
　　　　　　{
　　　　　　　　int IP = i+LE1;
　　　　　　　　COMP T=A[IP]*U;
　　　　　　　　A[IP]=A[i]-T;//蝶形计算
　　　　　　　　A[i]=A[i]+T;
　　　　　　　　i+=LE;
　　　　　　}
　　　　　　
　　　　　　U=U*W;//不同的W次幂
　　　　}
　　}
}
void initX(int N,COMP x[],float F,float T)
{
　　for(int i=0;i< N;i++)
　　{
　　　　x[i] = COMP(cos(2*pi*F*T*i),0);
　　}
}　　　　　　　　　　　　

3.2 子函数代码实现

/********************************************************************
*　 Name　　: FuncHd
*　　Function: Hd()--Required frequency response function
*
*
*********************************************************************/
COMP FuncHd(double LowLimit,double UpperLimit,COMP x)
{
　　if(x.real()>UpperLimit||x.real() < LowLimit)
　　　　return 0;
　　else
　　　　return 1;
　　
}
void FIR(double LowLimit,double UpperLimit,int N,COMP Hn[])
{
　　int M = 2*N;
　　for(int i=0;i < N;i++)
　　{
　　　　Hn[i] = COMP(0,0);
　　　　for(int k=0;k < M;k++)
　　　　{
　　　　　　COMP C = COMP(cos(2*pi*i*k/(double)M),sin(2*pi*i*k/(double)M));
　　　　　　Hn[i] += C*FuncHd(LowLimit,UpperLimit,COMP(cos(2*pi*k/(double)M),sin(2*pi*k/(double)M)));
　　　　}
　　　　Hn[i] = Hn[i]*COMP(1/(double)M,0);
　　}
}4、结束语

基本算法参考《数字信号处理基础及试验》--王树勋主编。虽然现在DSP算法都有很好C语言实现。但是能够通过自己动手编写代码加深对基础知识的掌握，对自己进行数据采集器件的控制还是有很多益处的。