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

typedef double 
(* FFuncs)(int); //h1(x) double
funch1(int
n) { double
fbase
	= (double)4/(double)5; double fr
	= std::pow(fbase, n); return fr;
	} //h2(x)
double
funch2(int
n) { double
fpi
	= 3.1415927; return 0.5*sin((double)0.5*n);
	} //y(n)
//y(n)=
sum(x(m)*y(n-m))
m=0..n double
Calcy(double x[],double h[],int n) {
double
	yvalue =  0; for(int
	m= 0;m<=n;m++)
	{
		yvalue += x[m]*h[n-m];
	}

	return yvalue;
}

程序界面,具体实现见注释及代码:

图2 DFT与FFT实现界面

void CFFTConversionDlg::OnBnClickedBtncal() 
{
	this->UpdateData(TRUE);
	int
	nN = this->m_NumN; float
	fF = this->m_NumF; float
	fT = this->m_NumT; bool

	bIsTimesof2 = false; 
	for(int i= 0;i<100;i++)
	{
	        if(nN==(2 < < i))
		{
			bIsTimesof2 = true;
			break;
		}
	}
	if(!bIsTimesof2)
	{
		AfxMessageBox("N请输入一个以2为底的幂级数!");
		this->GetDlgItem(IDC_EDTN)->SetFocus();
		return;
	}
	COMP* x = new COMP[nN];//x(n)
	COMP* X = new COMP[nN];//X(k) 
	initX(nN,x,fF,fT);
	CButton* pRadio = (CButton*)this->GetDlgItem(IDC_RADIODFT);


	if(pRadio->GetCheck())
	{
		DFT(nN,x,X);
	}
	else
	{
		FFT(nN,x,X);		
	}

	char buffer[256];
	COMP source = X[nN-1];
	sprintf(buffer,"%f+%fi",source.real(),source.imag());
	CWnd* pwnd = this->GetDlgItem(IDC_EDTRET);
	pwnd->SetWindowText(buffer);

	CListCtrl* pList=(CListCtrl*) this->GetDlgItem(IDC_LIST1);
	CListOper oper;
	oper.FillList(*pList,nN,x,X);

	delete[] x; 
	delete[] X;

}

/*****************************************
*
* 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*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 = 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-T;//蝶形计算
				A=A+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 = 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 = 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 += C*FuncHd(LowLimit,UpperLimit,COMP(cos(2*pi*k/(double)M),sin(2*pi*k/(double)M)));

		}
		Hn = Hn*COMP(1/(double)M,0);
	}

}

4、结束语

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

http://www.vckbase.com/document/viewdoc/?id=1644附：源代码