fftview.cpp

数字信号处理试验

		m_iCount = paraSetDlg.count;
		int i;
		for (i=0; i<m_iCount; i++ )
		{
			xreal[i] = CreatAttriReal(i);
			ximg[i] = CreatAttriImg(i);
		}
		CreateDFT(Xreal, Ximg, xreal, ximg,m_iCount);
		CreateAbs(XAbs, Xreal, Ximg, m_iCount);
		CreateArg(XArg, Xreal, Ximg, m_iCount);
		//调用画图对话框和求傅立叶变换
		ShowClineDlg(4, 1);
	}		
}

//显示绘图对话框
void CFFTView::ShowClineDlg(int n, int index2)
{
	int i;
// 	CLineDlg LineDlg;            //初始化
// 	LineDlg.LineCount = m_iCount;
// 	LineDlg.index = n ;
// 	LineDlg.index2 = index2;
// 	LineDlg.xPreal = xreal;
// 	LineDlg.xPimg = ximg;
// 	LineDlg.XPreal = Xreal;
// 	LineDlg.XPimg = Ximg;
// 	LineDlg.XPabs = XAbs;
// 	LineDlg.XParg = XArg;
// 	LineDlg.para1 = m_fPara1;
// 	LineDlg.para2 = m_fPara2;
// 	LineDlg.xMaxreal = mymax(xreal, m_iCount);
// 	LineDlg.xMaximg = mymax(ximg, m_iCount);
// 	LineDlg.XMaxreal = mymax(Xreal, m_iCount);
// 	LineDlg.XMaximg = mymax(Ximg, m_iCount);
// 	LineDlg.XMaxAbs = mymax(XAbs, m_iCount);
// 	LineDlg.XMaxArg = mymax(XArg, m_iCount);
// 	UpdateData(false);
// 	LineDlg.DoModal();	

	CLineDlg *LineDlg = new CLineDlg;
	LineDlg->Create(IDD_DIALOG1, this);
	LineDlg->LineCount = m_iCount;
	LineDlg->index = n;
	LineDlg->index2 = index2;

	LineDlg->xPreal = new double [m_iCount];
	LineDlg->xPimg = new double [m_iCount];
	LineDlg->XPreal = new double [m_iCount];
	LineDlg->XPimg = new double [m_iCount];
	LineDlg->XPabs = new double [m_iCount];
	LineDlg->XParg = new double [m_iCount];
	
	for (i = 0; i<m_iCount; i++)
	{
		LineDlg->xPreal[i] = xreal[i];
	}
	
	for (i = 0; i<m_iCount; i++)
	{
		LineDlg->xPimg [i] = ximg[i];
	}
	
	for (i = 0; i<m_iCount; i++)
	{
		LineDlg->XPreal[i] = Xreal[i];
	}
	
	for (i = 0; i<m_iCount; i++)
	{
		LineDlg->XPimg[i] = Ximg[i];
	}
	
	for (i = 0; i<m_iCount; i++)
	{
		LineDlg->XPabs [i] = XAbs[i];
	}
	
	for (i = 0; i<m_iCount; i++)
	{
		LineDlg->XParg[i] = XArg[i];
	}

	LineDlg->para1 = m_fPara1;
	LineDlg->para2 = m_fPara2;
	LineDlg->xMaxreal = mymax(xreal, m_iCount);
	LineDlg->xMaximg = mymax(ximg, m_iCount);
	LineDlg->XMaxreal = mymax(Xreal, m_iCount);
	LineDlg->XMaximg = mymax(Ximg, m_iCount);
	LineDlg->XMaxAbs = mymax(XAbs, m_iCount);
	LineDlg->XMaxArg = mymax(XArg, m_iCount);
	LineDlg->AddList();
	UpdateData(false);
	LineDlg->ShowWindow(SW_SHOW);
}
//高斯序列
double CFFTView::CreatGaussReal(int num, double para1, double para2)
{
	double ans;
	ans = (double)num - para1;
	ans *= ans;
	ans /= para2;
	ans = (double)exp(-ans);
	return ans;
}

double CFFTView::CreatGaussImg(int num, double para1, double para2)
{
	return 0;
}

//衰减正弦
double CFFTView::CreatAttsinReal(int num, double para1, double para2)  //para1-a, para2-f
{
	double ans;
	ans = (double)(sin(2*PI*para2*num)*exp(-1*para1*num));
	return ans;
}

double CFFTView::CreatAttsinImg(int num, double para1, double para2)
{
	return 0;
}

//三角波序列
double CFFTView::CreatTriReal(int num)
{
	int ans;
	if(num >= 0 && num <=3)
		ans = num + 1;
	else if (num >= 4 && num <= 7)
		ans = 8 - num;
	else 
		ans = 0;
	return (double)ans;
}

double CFFTView::CreatTriImg(int num)
{
	return 0;
}
//反三角序列
double CFFTView::CreatAttriReal(int num)
{
	int ans;
	if(num >= 0 && num <=3)
		ans = 4 - num;
	else if (num >= 4 && num <= 7)
		ans = num - 3;
	else 
		ans = 0;
	return (double)ans;
}

double CFFTView::CreatAttriImg(int num)
{
	return 0;
}

//求最大绝对值
double CFFTView::mymax(double *a, int n)
{
	int i;
 	double *b = a;
 	double num = 0;
 	for (i=0; i<n; i++)
 	{
 		if(fabs(*b) > num)
 			num = fabs(*b);
 		b++;
 	}
 	return num;
}

// 倒位排序函数
// j为正常序列的下标
// nu为长度的log2（长度），长度为16，则为4
// nu在以后要改为长度 count
int CFFTView::ibitr(int j, int count)
{
	int b,j1,i,j2,nu;
	int nu1;
	nu = 0;
	nu1 = count;
	while (nu1 / 2)
	{
		nu1 /= 2;
		nu ++;
	}	
	j1 = j;
	b = 0;
	for (i=1; i<=nu; i++)
	{
		j2 = j1/2;
		b = b*2 + (j1-2*j2);
		j1 = j2;
	}
	return b;
}

// 生成FFT变换的实部和虚部，
// 入口参数：数组名和长度
// Xreal:    存储FFT结果实部的数组名,相当与返回值
// Ximag:    存储FFT结果虚部的数组名,相当于返回值
// xreal:    存储时间序列实部的数组名
// ximga:    存储时间序列虚部的数组名
// n:	     数组的长度
void CFFTView::CreateFFT(double *Xreal, double *Ximag, double *xreal, double *ximag, int n)
{	
	int i,j,k,l,m,deta,p,a;	
	double treal,timag,temp;
	//赋值函数
	for (i = 0; i<n; i++)
	{
		Xreal[i] = xreal[i];
		Ximag[i] = ximag[i];
	}
	//倒位函数，交换数组内容
	for (i = 0; i<n; i++)
	{
		j = ibitr(i, n);
		if (j > i)
		{
			temp = Xreal[i];
			Xreal[i] = Xreal[j];
			Xreal[j] = temp;
			temp = Ximag[i];
			Ximag[i] = Ximag[j];
			Ximag[j] = temp;
		}
	}
	//FFT蝶形运算
	a = 0;
	k = n;
	while (k = k/2)
	{
		deta = (int)pow(2,a);		
		l = 0; 
		m = 0;
		for(i = 0; i < k; i++)
		{			
			for(j = 0; j < deta; j++)   //　l　和　l+deta 的蝶形运算
			{
				p = j * k;
				treal = Xreal[l+deta] * cos(2*PI/n*p)
					+ Ximag[l+deta] * sin(2*PI/n*p);
				timag = Ximag[l+deta] * cos(2*PI/n*p)
					- Xreal[l+deta] * sin(2*PI/n*p);
				Xreal[l+deta] = Xreal[l] - treal;
				Ximag[l+deta] = Ximag[l] - timag;
				Xreal[l] = Xreal[l] + treal;
				Ximag[l] = Ximag[l] + timag;
				l++;
			}
			m += deta * 2;
			l = m;
		}
		a++;
	}
	modify();
}

// 生成DFT变换的实部和虚部，
// 入口参数：数组名和长度
// Xreal:    存储FFT结果实部的数组名,相当与返回值
// Ximag:    存储FFT结果虚部的数组名,相当于返回值
// xreal:    存储时间序列实部的数组名
// ximga:    存储时间序列虚部的数组名
// n:	     数组的长度
void CFFTView::CreateDFT(double *Xreal, double *Ximag, double *xreal, double *ximag, int n)
{
	int k;
	int i;
	for (k=0; k<n; k++)
	{
		Xreal[k] = 0;
		Ximag[k] = 0;
		for (i=0; i<n; i++)
		{
			Xreal[k] += xreal[i]*cos(2*PI*k*i/n) + ximag[i]*sin(2*PI*k*i/n);
			Ximag[k] += ximag[i]*cos(2*PI*k*i/n) - xreal[i]*sin(2*PI*k*i/n);
		}
	}
	modify();
}
// 生成模值
// 入口参数：数组名和长度
// Xabs:     存储模值的数组名，相当与返回值
// Xreal:    存储实部的数组名
// Ximag:    存储虚部的数组名
// n:	     数组的长度
void CFFTView::CreateAbs(double *Xabs, double *Xreal, double *Ximag, int n)
{
	int i;
	for(i=0; i<n; i++)
	{
		Xabs[i] = sqrt(Xreal[i]*Xreal[i] + Ximag[i]*Ximag[i]);
	}
}

// 入口参数：数组名和长度
// Xarg:     存储角度的数组名，相当与返回值
// Xreal:    存储实部的数组名
// Ximag:    存储虚部的数组名
// n:	     数组的长度
void CFFTView::CreateArg(double *Xarg, double *Xreal, double *Ximag, int n)
{
	int i;
	for(i=0; i<n; i++)
	{
		Xarg[i] = atan2(Ximag[i],Xreal[i]);
	}
}


void CFFTView::modify()
{
	int i;
	for (i=0; i<m_iCount; i++)
	{
		if (fabs(Xreal[i]) < 1e-6)
		{
			Xreal[i] = 0.0;
		}
		if (fabs(Ximg[i]) < 1e-6)
		{
			Ximg[i] = 0.0;
		}
	}
}