expression.cpp

这是我学习数字信号处理时编写的程序！功能就是显示一些常用的函数的图形！并可以做FFT计算！

// Expression.cpp: implementation of the CExpression class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "Zhangjg_dsp.h"
#include "Expression.h"
#include "Zhangjg_dspDoc.h"
#include "math.h"
#define PI 3.1415926
#include <complex>
using namespace std;
#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CExpression::CExpression()
{

}

CExpression::~CExpression()
{

}

void CExpression::AttenuationSin(int n, double a, double f, double *X)
{
       for(int i=0;i<=n-1;i++)
	   X[i]=exp(-a*i)*sin(2*PI*f*i);
}

void CExpression::Gauss( int p, int q,double *X,int n)
{   
       for(int i=0;i<=n-1;i++)
		X[i]=exp(-(i-p)*(i-p)/q);
}
void  CExpression::Fft(complex<double> *FD,int M)      //FD use for put in and puting out;
{
	long N;       //length;
	// 循环变量
   
	int i; //级数i
    long n;//时域n
	
	// 中间变量
	long    P;
	
	// 角度
	double	angle;
	
	complex<double> *W,*X2,*X1,*X;
	
	// 计算付立叶变换点数
	N = 1 << M;
	complex<double> Temple;//中间变量；
	X1 = new complex<double>[N];    //初始数据；	
	// 将时域点写入X1
	memcpy(X1, FD, sizeof(complex<double>) * N);

	if(N==4)                 //be careful when N ==4;the mistake is here;so it is nasisery
	{
		W  = new complex<double>[N *N];
			// 计算加权系数,即旋转因子；
	for(n = 0; n < N *N; n++)
	{
		angle = -n * PI * 2 / N;
		W[n] = complex<double> (cos(angle), sin(angle));
	}
		complex<double> Temple1;//中间变量；
		for(int k=0;k<=3;k++)
		{

			Temple=complex<double>(0,0);
			for(int n=0;n<=3;n++)
			{
				int p=k*n;
				Temple1=X1[n]*W[p];
				Temple=Temple+Temple1;
			}
			FD[k]=Temple;
		}
		    // 释放内存
	delete W;
	delete X1;
	}
	else
	{
			// 分配运算所需存储器
	W  = new complex<double>[N / 2];

	X2 = new complex<double>[N];    //中间计算时用到的
			// 计算加权系数,即旋转因子；
	for(n = 0; n < N / 2; n++)
	{
		angle = -n * PI * 2 / N;
		W[n] = complex<double> (cos(angle), sin(angle));
	}
		//对时域点进行排序；
	for(long I=1,J=0;I<=N-1;I++)
	{

		for(i=1;i<=M;i++)
		{
			
			if(J<1<<(M-i))
			{
				J=J+(1<<(M-i));
				if(J>I)
				{
				Temple=X1[I];
				X1[I]=X1[J];
                X1[J]=Temple;
				}
				 break;
			}
			else 
			{
				J=J-(1<<(M-i));
		     	continue;
			}
		
		}
	}


	for(i=1;i<=M;i++)   //级数运算；
	{
		long B=1<<(i-1);     //输入数据的间距；
		for(n=0;n<=N-1;)
		{
           for(long J=0;J<=(1<<(i-1))-1;J++)
		   {

			   P=J*(1<<(M-i));
			   X2[n]=X1[n]+X1[n+B]*W[P];
			   X2[n+B]=X1[n]-X1[n+B]*W[P];
			 n++;
		   }
		   n=n+B;
		} 
		X  = X1;
		X1 = X2;
		X2 = X;
	}
    memcpy(FD, X1, sizeof(complex<double>) * N);
    // 释放内存
	delete W;
	delete X1;
	delete X2;
	}
}
void  CExpression::Dft(complex<double> *FD,int M)  
{	
	// 计算付立叶变换点数
	long N = 1 << M;
	complex<double> *X;
	complex<double> Temple;//中间变量；
	X= new complex<double>[N];    //初始数据；	
    for(int k=0;k<=N-1;k++)
	{
		Temple=complex<double>(0,0);
		X[k]=complex<double>(0,0);
		for(int n=0;n<=N-1;n++)
		{
			Temple=FD[n]*exp(complex<double>(0,-2*PI*k*n/N));
			X[k]+=Temple;
		}
	}
    memcpy(FD, X, sizeof(complex<double>) * N);
    // 释放内存
	delete X;
}
void CExpression::Rectangle(int n,double *X)
	{
        for(int i=0;i<=n-1;i++)
		X[i]=1;
	}

void CExpression::Triangle(int n, double *X)
{   
	int t=n/8;
for(int i=0;i<=t-1;i++)
{    
	for(int k=0,j=8*i;j<8*(i+1),k<8;j++,k++)
		if(k<=3)
		X[j]=k+1;
		else
		X[j]=8-k;
}
}

void CExpression::ReTriangle(int n, double *X)
{
 int t=n/8;
for(int i=0;i<=t-1;i++)
{    
	for(int k=0,j=8*i;j<8*(i+1),k<8;j++,k++)
		if(k<=3)
		X[j]=4-k;
		else
		X[j]=k-3;
}
}
void CExpression::DrawCoordinate(int stepx, int stepy,int pointx, int pointy,CDC* pDC, CZhangjg_dspDoc* pDoc)
{
		//CZhangjg_dspDoc *pDoc=GetDocument();//获取活动文档的指针；

	    pDC->MoveTo(pointx,pointy+stepy*5);//y轴的长
		pDC->LineTo(pointx,pointy-stepy*5);
     	pDC->LineTo(pointx-4,pointy-stepy*5+10);//画箭头
		pDC->MoveTo(pointx,pointy-stepy*5);
		pDC->LineTo(pointx+4,pointy-stepy*5+10);
                                       
		    pDC->MoveTo(pointx-stepx*pDoc->m_N-100,pointy);
		    pDC->LineTo(pointx+stepx*pDoc->m_N-100,pointy);
			pDC->LineTo(pointx+stepx*pDoc->m_N-110,pointy-7);
            pDC->MoveTo(pointx+stepx*pDoc->m_N-100,pointy);
			pDC->LineTo(pointx+stepx*pDoc->m_N-110,pointy+7);//x轴的长，画箭头 need to be inproment add parameter;
		
			CString str;
			for(int i=-3;i<32;i++)
			{
				if(i!=0)
				{
	   str.Format("%d", i);
	   pDC->TextOut(pointx+stepx*i-2,pointy+2,str);//x轴的细节
				}

			}

            for(i=-3;i<=4;i++)
			{
				if(i!=0)
				{
				str.Format("%d",i);
               pDC->TextOut(pointx-15,pointy-stepy*i-6,str);
				}
			}

			       pDC->TextOut(pointx+5,pointy-stepy*4,"x(n)");
				   pDC->TextOut(pointx+stepy*20-2,pointy+2,"t");
				   pDC->TextOut(pointx-stepx*10,pointy-stepy*2,"数字信号处理!!");
}