分类: C/C++
2008-08-04 09:36:19
在学习数字信号处理算法程序中用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 子函数代码
/***************************************** * * 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语言实现。但是能够通过自己动手编写代码加深对基础知识的掌握,对自己进行数据采集器件的控制还是有很多益处的。
下载本文示例代码