几个数字信号处理算法程序
作者：liu_sir
摘要
在学习数字信号处理算法程序中用VC编写的几个通用算法程序。
关键词 离散卷积 FIR
在学习信号处理的过程中,看到书上的大部分算法都是用Fortan或者Basic实现,于是自己试验着用VC实现了一下。
1、卷积计算
离散卷积公式的算法实现
图1 卷积计算界面
1.1 主程序代码(省略了部分不关键代码)
void CInterVolveDlg::CalTheNumByArray()
{
this->UpdateData(TRUE);
FFuncs
funcs[2] = {funch1,funch2}; int
n = this->m_ValueN; double*
x = new double[2*(n+1)];//x(n) double*
y = new double[2*(n+1)];//y(n) double*
h = new double[2*(n+1)];//h(n) //1.init
x(n),h(n),y(n) CButton*
pbtn = (CButton*) this->GetDlgItem(IDC_RADIO1); int
nChoseItem = 0;//函数选择 if(pbtn->GetCheck())
{
nChoseItem
= 0; }
else
{
nChoseItem
= 1; }
for(int
i= 0;i<2*(n+1);i++)
{
if(i< n+1)
{
x[i] = 1;
h[i] = funcs[nChoseItem](i);
}
else
{
x[i] = 0;
h[i] = 0;
}
}
//2.y(i)=SUM(x(m)*h(i-m)) m=0..i
for(i=0;i<2*(n+1);i++)
{
y[i] = Calcy(x,h,i);
}
//显示结果
delete[] x;
delete[] y;
delete[] h;
}
1.2 各个子函数实现
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实现程序界面,具体实现见注释及代码:
图2 DFT与FFT实现界面
2.1 主程序代码
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;
}
2.2 子函数代码
说明:其中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
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语言实现。但是能够通过自己动手编写代码加深对基础知识的掌握，对自己进行数据采集器件的控制还是有很多益处的。
