firfilterdlg.cpp

用Visual c++写的有限长单位脉冲响应滤波器

	for(i=1;i<=N;i++)
	{
		temp=data1[i].real*data1[i].real+data1[i].image*data1[i].image;
		if(temp>=maxtemp)
			maxtemp=temp;
	}
	
//	CRect rect3;

//	pwnd->GetClientRect(rect3);  
	int height=rect1.bottom-rect1.top;
	int width=rect1.right-rect1.left;

	double stepcols;
	stepcols=(double)(height-20)/maxtemp;   //140与picture控件的高度有关，也就是每一步长所占的像素

	for(i=1;i<N/2;i++)
	{
		temp1=(data1[i].real*data1[i].real+data1[i].image*data1[i].image);
		temp2=(data1[i+1].real*data1[i+1].real+data1[i+1].image*data1[i+1].image);
		double step;
/*		for(i=1;i<=N;i++)
		{
			step[i]=2*PI*i/(double)N;
		}
*/
		step=width/N;   //380是绘图宽度

		int x1=(int)(5+(i-1)*step);
		int y1=height-(int)(temp1*stepcols); //160是绘图高度

		int x2=(int)(5+i*step);
		int y2=height-(int)(temp2*stepcols);

		PCDC->MoveTo(x1,y1);
		PCDC->LineTo(x2,y2);
	}
	PCDC->SetBkColor(RGB(255,0,255));   //设置的是文字的背景颜色
	PCDC->SetTextColor(RGB(0,0,0));
	PCDC->TextOut(0,0,"正弦输入信号加干扰噪声的频谱");
}

/********************************画滤波后信号的波形************************************************/
void CFIRFilterDlg::DrawOutputWave(int N, Complex data1[])
{
	CWnd *pwnd=GetDlgItem(IDC_OUTPUT_WAVE);
	CDC *PCDC=pwnd->GetDC();

	CRect rect1;
	pwnd->GetClientRect(rect1); 

	CBrush brwhite(RGB(255,255,255)); 
	PCDC->FillRect(rect1,&brwhite);

	/****************画网格的程序段*************************/

	int nwidth=0,nheight=0,i=0,m=0,n=0;

	CPen Dot_pen(PS_DOT,0,RGB(0,0,0));
	CPen Dot_penCenter(PS_DOT,0,RGB(0,0,255));
	CPen Solid_pen(PS_SOLID,0,RGB(0,0,0));
//	CString str;

//	CBrush brwhite(RGB(255,255,255));

//	CRect rect1;

//	pwnd->GetClientRect(rect1);  
	PCDC->SelectObject(&Solid_pen);   //选入设备环境
//	PCDC->FillRect(rect1,&brwhite);


	//画边框
	PCDC->MoveTo(rect1.left,rect1.top);
	PCDC->LineTo(rect1.right,rect1.top);

	PCDC->MoveTo(rect1.right,rect1.top);
	PCDC->LineTo(rect1.right,rect1.bottom);

	PCDC->MoveTo(rect1.left,rect1.bottom);
	PCDC->LineTo(rect1.right,rect1.bottom);

	PCDC->MoveTo(rect1.left,rect1.top);
	PCDC->LineTo(rect1.left,rect1.bottom);

	PCDC->SelectObject(&Dot_penCenter);   //选入设备环境
	//PCDC->FillRect(rect1,&brwhite);

	int CenterLinePos= (rect1.bottom-rect1.top)/2;

	PCDC->MoveTo(rect1.left,CenterLinePos);
	PCDC->LineTo(rect1.right,CenterLinePos);

	//画框内的网格
	//将X和Y坐标分别等分为N1和N2份
	nwidth=(int)((rect1.right-rect1.left)/30);
	nheight=(int)((rect1.bottom-rect1.top)/30);

	PCDC->SelectObject(&Dot_pen);   //选入设备环境

	//画中心线上方的横线
	m=CenterLinePos;
	for(i=1;i<=nheight/2;i++)
	{
		m-=30;
		if(m > rect1.top)
		{
			PCDC->MoveTo(0,m);
	        PCDC->LineTo(rect1.right-rect1.left,m);
		}
	}

	//画中心线下方的横线
	m=CenterLinePos;
	for(i=1;i<=nheight/2;i++)
	{
		m+=30;
		if(m < rect1.bottom)
		{
			PCDC->MoveTo(0,m);
			PCDC->LineTo(rect1.right-rect1.left,m);
		}
	}

	//画竖线
	for(i=1;i<=nwidth;i++)
	{
		n+=30;
		PCDC->MoveTo(n,0);
		PCDC->LineTo(n,rect1.bottom-rect1.top);
	}

	PCDC->SetBkColor(RGB(255,0,255));   //设置的是文字的背景颜色
	PCDC->SetTextColor(RGB(0,0,0));
	PCDC->TextOut(0,0,"滤波后的波形");


	/******************************画刻度的程序段***************************/

	CFont m_ScaleFont;
	CString str;
	m_ScaleFont.CreateFont(15,0,0,0,FW_THIN,FALSE,FALSE,FALSE,DEFAULT_CHARSET,OUT_DEFAULT_PRECIS,CLIP_DEFAULT_PRECIS,DEFAULT_QUALITY,DEFAULT_PITCH,"Arial");

	PCDC->SelectObject(&m_ScaleFont);
	PCDC->SetTextColor(RGB(0,0,0));
	PCDC->SetBkMode(TRANSPARENT);
	PCDC->SetTextAlign(TA_RIGHT|TA_BASELINE);

	//画坐标原点
	PCDC->TextOut(rect1.left-2,(int)((rect1.bottom-rect1.top)/2)+5,"0.0");
	/**************画正半轴********************/


	for(i=1;i<=nheight/2;i++)
	{
		int nPosY1=(rect1.bottom-rect1.top)/2- i*30;
		if(nPosY1 >rect1.top)
		{
			str.Format("%.1f",i*1.0);
		    PCDC->TextOut(rect1.left-2,nPosY1+5,str);
		}
	}


	/************画负半轴*******************/


	for(i=1;i<=nheight/2;i++)
	{
		int nPosY2=(rect1.bottom-rect1.top)/2 + i*30;
		if(nPosY2 < rect1.bottom)
		{
			str.Format("%.1f",-i*1.0);
		    PCDC->TextOut(rect1.left-2,nPosY2+5,str);
		}		
	}


	/************画时间轴的坐标刻度**********/
	PCDC->SetTextAlign(TA_LEFT|TA_TOP);
	PCDC->TextOut(rect1.left-1,rect1.bottom+2,CString("0"));//刻度规定原点

	PCDC->SetTextAlign(TA_CENTER|TA_TOP);

	n=0;
	for(i=1;i<=nwidth;i++)
	{
		n+=30;
		str.Format("%d",i*20);
		PCDC->TextOut(n+2,rect1.bottom+2,str);
	}


	/******************画波形的程序段**********************/

	CPen backpen1(PS_SOLID,0,RGB(0,0,0));
	CPen *pold1=NULL;
	pold1=PCDC->SelectObject(&backpen1);

	int width=rect1.right-rect1.left;
	int height=rect1.bottom-rect1.top;

	double temp1,temp2;
//	int i;
/*	double maxtemp=0.0,temp;

	for(i=1;i<=N;i++)
	{
		temp=(data1[i].real>0)?data1[i].real:(0-data1[i].real);
		if(temp>=maxtemp)
			maxtemp=temp;
	}
*/
//	double stepcols;
//	double stepcols=(double)(height/2-10)/maxtemp;

	double step=(double)(width-10)/N;

	for(i=1;i<N;i++)
	{
		temp1=data1[i].real;
		temp2=data1[i+1].real;
		
		int x1=(int)(step*(i-1));
		int y1=height/2- (int)(temp1*30);

		int x2=(int)(step*i);
		int y2=height/2-(int)(temp2*30);

		if(x1 < rect1.right && x2 < rect1.right && y1 < rect1.bottom && y1 > rect1.top && y2 < rect1.bottom && y2 > rect1.top)
		{
			PCDC->MoveTo(x1,y1);
		    PCDC->LineTo(x2,y2);
		}
		
	}
	
}

/******************************画滤波后信号的频谱的函数*********************************************/
void CFIRFilterDlg::DrawOutputSpectrum(int N,Complex data1[])
{
	CWnd *pwnd=GetDlgItem(IDC_OUTPUTSPECTRUM);
	CDC *PCDC=pwnd->GetDC();

	CRect rect1;
	pwnd->GetClientRect(rect1); 

	CBrush brwhite(RGB(255,255,255)); 
	PCDC->FillRect(rect1,&brwhite);

	/****************画网格的程序段*************************/

	int nwidth=0,nheight=0,i=0,m=0,n=0;

	CPen Dot_pen(PS_DOT,0,RGB(0,0,0));
	CPen Dot_penCenter(PS_DOT,0,RGB(0,0,255));
	CPen Solid_pen(PS_SOLID,0,RGB(0,0,0));
//	CString str;

//	CBrush brwhite(RGB(255,255,255));

//	CRect rect1;

//	pwnd->GetClientRect(rect1);  
	PCDC->SelectObject(&Solid_pen);   //选入设备环境
//	PCDC->FillRect(rect1,&brwhite);


	//画边框
	PCDC->MoveTo(rect1.left,rect1.top);
	PCDC->LineTo(rect1.right,rect1.top);

	PCDC->MoveTo(rect1.right,rect1.top);
	PCDC->LineTo(rect1.right,rect1.bottom);

	PCDC->MoveTo(rect1.left,rect1.bottom);
	PCDC->LineTo(rect1.right,rect1.bottom);

	PCDC->MoveTo(rect1.left,rect1.top);
	PCDC->LineTo(rect1.left,rect1.bottom);

	PCDC->SelectObject(&Dot_penCenter);   //选入设备环境
	//PCDC->FillRect(rect1,&brwhite);

	int CenterLinePos= (rect1.bottom-rect1.top)/2;

	PCDC->MoveTo(rect1.left,CenterLinePos);
	PCDC->LineTo(rect1.right,CenterLinePos);

	//画框内的网格
	//将X和Y坐标分别等分为N1和N2份
	nwidth=(int)((rect1.right-rect1.left)/30);
	nheight=(int)((rect1.bottom-rect1.top)/30);

	PCDC->SelectObject(&Dot_pen);   //选入设备环境

	//画中心线上方的横线
	m=CenterLinePos;
	for(i=1;i<=nheight/2;i++)
	{
		m-=30;
		if(m > rect1.top)
		{
			PCDC->MoveTo(0,m);
	        PCDC->LineTo(rect1.right-rect1.left,m);
		}
	}

	//画中心线下方的横线
	m=CenterLinePos;
	for(i=1;i<=nheight/2;i++)
	{
		m+=30;
		if(m < rect1.bottom)
		{
			PCDC->MoveTo(0,m);
			PCDC->LineTo(rect1.right-rect1.left,m);
		}
	}

	//画竖线
	for(i=1;i<=nwidth;i++)
	{
		n+=30;
		PCDC->MoveTo(n,0);
		PCDC->LineTo(n,rect1.bottom-rect1.top);
	}


	/******************画波形的程序段******************************/

	CPen backpen1(PS_SOLID,1,RGB(0,0,0));
	CPen *pold1=NULL;
	pold1=PCDC->SelectObject(&backpen1);

	double temp,temp1,temp2;
	double maxtemp=0.0;
//	int i;
	for(i=1;i<=N;i++)
	{
		temp=data1[i].real*data1[i].real+data1[i].image*data1[i].image;
		if(temp>=maxtemp)
			maxtemp=temp;
	}
	


	int height=rect1.bottom-rect1.top;
	int width=rect1.right-rect1.left;

	double stepcols;
	stepcols=(double)(height-20)/maxtemp;   //140与picture控件的高度有关，也就是每一步长所占的像素

	for(i=1;i<N/2;i++)
	{
		temp1=(data1[i].real*data1[i].real+data1[i].image*data1[i].image);
		temp2=(data1[i+1].real*data1[i+1].real+data1[i+1].image*data1[i+1].image);
		double step;
/*		for(i=1;i<=N;i++)
		{
			step[i]=2*PI*i/(double)N;
		}
*/
		step=width/N;   //380是绘图宽度

		int x1=(int)(5+(i-1)*step);
		int y1=height-(int)(temp1*stepcols); //160是绘图高度

		int x2=(int)(5+i*step);
		int y2=height-(int)(temp2*stepcols);

		PCDC->MoveTo(x1,y1);
		PCDC->LineTo(x2,y2);
	}

	PCDC->SetBkColor(RGB(255,0,255));   //设置的是文字的背景颜色
	PCDC->SetTextColor(RGB(0,0,0));
	PCDC->TextOut(0,0,"滤波后信号的频谱");


}

/************************************** FFT函数****************************************************/
void CFIRFilterDlg:: FFT(int N, int M, Complex data2[])
{
	int j=1;
	int i=1;
	int K;
	do
	{
		if(i<j)
		{
			Complex T;
			T.real = data2[j].real;
			T.image = data2[j].image;
			data2[j].real=data2[i].real;
			data2[j].image=data2[i].image;
			data2[i].real=T.real;
			data2[i].image=T.image;
		}
		K=N/2;
		while(K<j)
		{
			j-=K;
			K/=2;
		}
		j+=K;
		i++;
	}while(i<N);


	int LE,LE1,IP;
	int L=1;

	do
	{
		LE=(int)pow(2,L);	
		LE1=LE/2;
		Complex U={1.0,0.0};
		Complex W={cos(0-PI/(double)LE1),sin(0-PI/(double)LE1)};
		j=1;
		do
		{
			i=j;

			do
			{
				IP=i+LE1;
				Complex T;
				T.real= data2[IP].real*U.real-data2[IP].image*U.image;
				T.image = data2[IP].real*U.image+data2[IP].image*U.real;
				data2[IP].real = data2[i].real-T.real;
				data2[IP].image = data2[i].image-T.image;
				data2[i].real = data2[i].real+T.real;
				data2[i].image = data2[i].image+T.image;
				i+=LE;
			} while(i<=N);

			double tempreal=U.real;
			double tempimage=U.image;
			U.real = U.real*W.real-W.image*U.image;
			U.image = tempreal*W.image+tempimage*W.real;
			j++;
		} while(j<=LE1);

		L++;

	} while(L<=M);
}

/**********************************滤波器设置函数*********************************************/
void CFIRFilterDlg::FirFunction(Complex A[], Complex B[])   //此函数在MainFunction()被调用
{
	Complex windowdata[513];

	for(int j=0;j<513;j++)
	{
		windowdata[j].real=0;
	}

	CreateFirFilter(windowlength, m_filter_type.GetCurSel(), m_window_type.GetCurSel(), windowdata);
	DrawSYS1(windowlength,windowdata);            //画出系统响应的图形
	
	Complex windowdata1[513];
/*	for(int i=1;i<513;i++)
	{
		windowdata1[i]=windowdata[i];
	}
	int pointtemp=windowpoint;
	int ffttime=0;
	while(pointtemp!=0)
	{
		pointtemp/=2;
		ffttime++;
	}
	ffttime--;
	if(windowpoint>pow(2,ffttime))
	{
		ffttime++;
	}
	for(i=windowpoint+1;i<=pow(2,ffttime);i++)
	{
		windowdata1[i].real=0;
		windowdata1[i].image=0;
	}
	FFT((int)pow(2,ffttime),ffttime,windowdata1);

	DrawSYS2((int)pow(2,ffttime),windowdata1);  //画出系统响应频谱的图形
*/
	
	for(int i=1;i<513;i++)
	{
		windowdata1[i]=windowdata[i];
	}

	DrawSYS3(windowlength,windowdata1);          //画出系统响应幅度频谱的图形


	int max=200+windowlength;
    myconv(A, windowdata, B, max);

}

/*********************************滤波器的创建函数******************************************/
void CFIRFilterDlg::CreateFirFilter(int N, int FilterKind, int WindowKind, Complex inputdata2[])
{
        double W = 0.01*cutofffrequency;
        double W1 = 0.01*lowcutfrequency;
        double W2 = 0.01*highcutfrequency;
        double F1 = 1.0;
        double F2 = 1.0;

        for(int i = 1; i <= N; i++)
        {
          if(WindowKind == 0)//窗的选择,矩形窗
            F2 = 1.0;
          else if(WindowKind == 1)//三角形窗
              if(i <= (N + 1) / 2)
                  F2 = ((double)(i - 1) * (double)2) / (double)N;
               else
				F2 = (double)2 - ((double)(i - 1) * (double)2) / (double)N;
			else if(WindowKind==2) //汉宁窗
				F2=pow(sin(((double)(i-1) * PI) / (double)N), 6);
				//0.6可设定汉宁窗a,这里a=6
			else if(WindowKind==3) //海明窗
				F2 = 0.54 - 0.46 * cos((2*PI * (double)(i-1)) / (double)(N - 1));
				//可设定海明窗a,公式为a-(1-a)*cos(2*PI*i)/(N-1),这里a=0.54
           else if(WindowKind == 4)  //布拉克曼窗
               F2 = (0.42323 - 0.4975499999999999 * cos((2*PI * (double)(i-1)) / (double)(N - 1))) + 0.07922 * cos((4*PI * (double)(i-1)) / (double)(N - 1));
				//这里可设定布拉克曼窗的阶的系数，这里取3阶
				//四阶公式为F2 = ((0.35875000000000001 - 0.48829 * cos((2*PI * (double)i) / (double)N)) + 0.14127999999999999D * Math.cos((4*PI * (double)i) / (double)N)) - 0.011679999999999999 * Math.cos((6*PI * (double)i) / (double)N);

            if(i != (N + 1) / 2)
            {
                if(FilterKind == 0)//低通
                    F1 = sin(W * PI * (double)(i - 1 - (N - 1) / 2)) / (PI * (double)(i - 1 - (N - 1) / 2));
                if(FilterKind == 1)//高通
                    F1 = -sin(W * PI * (double)(i - 1 - (N - 1) / 2)) / (PI * (double)(i - 1 - (N - 1) / 2));
                if(FilterKind == 2)//带通	
                    F1 = sin(W2 * PI * (double)(i - 1 - (N - 1) / 2)) / (PI * (double)(i - 1 - (N - 1) / 2)) - sin(W1 * PI * (double)(i - 1 - (N - 1) / 2)) / (PI * (double)(i - 1 - (N - 1) / 2));
                if(FilterKind == 3)//带阻
                    F1 = -(sin(W2 * PI * (double)(i - 1 - (N - 1) / 2)) / (PI * (double)(i - 1 - (N - 1) / 2))) + sin(W1 * PI * (double)(i - 1 - (N - 1) / 2)) / (PI * (double)(i - 1 - (N - 1) / 2));

            } else
            {
                if(FilterKind == 0)  //低通
                    F1 = W;
                if(FilterKind == 1)  //高通
                    F1 = (double)1 - W;
                if(FilterKind == 2)  //带通  
                    F1 = W2 - W1;
                if(FilterKind == 3)  //带阻
                    F1 = ((double)1 + W1) - W2;
            }
            inputdata2[i].real =F1 * F2;