ID #5734

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

摘要

在学习数字信号处理算法程序中用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语言实现。但是能够通过自己动手编写代码加深对基础知识的掌握,对自己进行数据采集器件的控制还是有很多益处的。

本文配套源码


2011-07-01 18:29
阅读:
I'm VC , Just U know Y
本站部分文章来源于互联网,版权归原作者所有。

延伸阅读: