摘要
在学习数字信号处理算法程序中用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 < 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语言实现。但是能够通过自己动手编写代码加深对基础知识的掌握，对自己进行数据采集器件的控制还是有很多益处的。
http://www.vckbase.com/document/viewdoc/?id=1644
附：源代码