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

yuyang911220

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　　 FT
*　 Function isperse 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);
　　}
}