wavefilterdlg.cpp

IIR数字滤波

	PCDC->LineTo(310,75);
	//y(n)图像绘图区,X轴（10,75,310,75）,绘图高度为140，中心高度100,绘图宽度为300左右
	//-------------------------------

	pwnd=GetDlgItem(IDC_OUT_PIC2);
	PCDC=pwnd->GetDC();
	CRect rc4(0,0,320,200);
	PCDC->FillRect(rc4,&brwhite);
	PCDC->MoveTo(10,195);
	PCDC->LineTo(310,195);
	PCDC->SelectObject(&pen1);
	//y(n)图像绘图区,X轴（10,195,310,195）,绘图高度为190，绘图宽度为300左右

}

void CWavefilterDlg::DrawInPic1(int N, Complex data1[])
{

	CWnd *pwnd=GetDlgItem(IDC_IN_PIC1);
	CDC *PCDC=pwnd->GetDC();

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

	double temp;
	int i;
	for(i=1;i<=N;i++)
	{
		temp=data1[i].real;
		double step=300.0/N;
		int x=(int)(10+(int)(step*(i-1)));
		int y=75- (int)(temp*65);
		PCDC->MoveTo(x,75);
		PCDC->LineTo(x,y);
	}
}

void CWavefilterDlg::DrawInPic2(int N, Complex data1[])
{
	CWnd *pwnd=GetDlgItem(IDC_IN_PIC2);
	CDC *PCDC=pwnd->GetDC();

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

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

	for(i=1;i<=N;i++)
	{
		temp=sqrt(data1[i].real*data1[i].real+data1[i].image*data1[i].image);
		double step;
		step=300.0/N;
		int x=(int)(10+(i-1)*step);
		int y=195- (int)(temp*stepcols);
		PCDC->MoveTo(x,195);
		PCDC->LineTo(x,y);
	}
}

void CWavefilterDlg::DrawOutPic1(int N, Complex data1[])
{
	CWnd *pwnd=GetDlgItem(IDC_OUT_PIC1);
	CDC *PCDC=pwnd->GetDC();

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

	double temp;
	int i;
	double maxtemp=0.0;
	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;
	stepcols=70.0/maxtemp;
	for(i=0;i<N;i++)
	{
		temp=data1[i].real;
		double step=300.0/N;
		int x=(int)(10+(int)(step*(i-1)));
		int y=75- (int)(temp*stepcols);
		PCDC->MoveTo(x,75);
		PCDC->LineTo(x,y);
	}
}

void CWavefilterDlg::DrawOutPic2(int N, Complex data1[])
{
	CWnd *pwnd=GetDlgItem(IDC_OUT_PIC2);
	CDC *PCDC=pwnd->GetDC();

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

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

	for(i=1;i<=N;i++)
	{
		temp=sqrt(data1[i].real*data1[i].real+data1[i].image*data1[i].image);
		double step;
		step=300.0/N;
		int x=(int)(10+(i-1)*step);
		int y=195- (int)(temp*stepcols);
		PCDC->MoveTo(x,195);
		PCDC->LineTo(x,y);
	}
}

void CWavefilterDlg::MainFunction()
{
	int i=1;
	if(m_wave.GetCurSel()==1)
		m_zhankong_msg.Format("%d%%",zhankong);
	m_point_msg.Format("%d",pointnum);
	m_period_msg.Format("%2gT",period/10.0);
	m_phase_msg.Format("%1gPI",phase*0.01);
	m_pointperperiod_msg.Format("每周期点数：%d",(int)pointnum*10/period);
	UpdateData(false);
	DrawCoordinate();
	Complex *x=new Complex[1024];
	Complex *X=new Complex[pointnum+1];
	int ftype=m_wave.GetCurSel();
	InitData(ftype,x);

	for(i=0;i<=pointnum;i++)
		X[i]=x[i];

	DrawInPic1(pointnum,x);



	int pointtemp=pointnum;
	int ffttime=0;
	while(pointtemp!=0)
	{
		pointtemp/=2;
		ffttime++;
	}
	ffttime--;
	if(pointnum>pow(2,ffttime))
	{
		ffttime++;
	}
	for(i=pointnum+1;i<=pow(2,ffttime);i++)
	{
		x[i].real=0;
		x[i].image=0;
	}
	FFT((int)pow(2,ffttime),ffttime,x);
	DrawInPic2((int)pow(2,ffttime),x);
	

	int start=3;
	Complex *y = new Complex[1025];
	double wd=(double)phase*0.01;
	IIRFunction(X,y,pointnum,wd,m_jieshu.GetCurSel()+3);
	DrawOutPic1(pointnum,y);
	FFT((int)pow(2,ffttime),ffttime,y);
	DrawOutPic2((int)pow(2,ffttime),y);

}

void CWavefilterDlg::OnOK() 
{
	MainFunction();
	
}

void CWavefilterDlg::OnSelchangeWaveShape() 
{
	if(m_wave.GetCurSel()==1)
	{
		m_ZK1="占空比:";
		m_zhankong_msg.Format("%d%%",m_zhankong.GetScrollPos());
		CWnd *pWnd=GetDlgItem(IDC_ZHANKONG);
		pWnd->ShowWindow(SW_SHOW);
		UpdateData(false);
	}
	else
	{
		m_ZK1="";
		m_zhankong_msg="";
		CWnd *pWnd=GetDlgItem(IDC_ZHANKONG);
		pWnd->ShowWindow(SW_HIDE);
		UpdateData(false);
	}
	MainFunction();	
}

void CWavefilterDlg::OnClrpic() 
{
	DrawCoordinate();	
}

void CWavefilterDlg::OnSelchangeJieshu() 
{
	MainFunction();	
}

void CWavefilterDlg::IIRFunction(Complex A[], Complex B[], int N, double wd, int grade)
{
     for(int i = 1; i <= 256; i++)
	 {
            B[i].real=0.0;
			B[i].image=0.0;
	 }

        wd *= PI;
        double wa = (double)2 * tan(wd / (double)2);
        if(grade == 3)
        {
            double a0 = pow(wa, 3);
            double a1 = (double)3 * pow(wa, 3);
            double a2 = (double)3 * pow(wa, 3);
            double a3 = pow(wa, 3);
            double A0 = (double)8 + (double)8 * wa + (double)4 * pow(wa, 2) + pow(wa, 3);
            double A1 = ((double)-24 - (double)8 * wa) + (double)4 * pow(wa, 2) + (double)3 * pow(wa, 3);
            double A2 = ((double)24 - (double)8 * wa - (double)4 * pow(wa, 2)) + (double)3 * pow(wa, 3);
            double A3 = (((double)-8 + (double)8 * wa) - (double)4 * pow(wa, 2)) + pow(wa, 3);
            a0 /= A0;
            a1 /= A0;
            a2 /= A0;
            a3 /= A0;
            A1 /= A0;
            A2 /= A0;
            A3 /= A0;
            B[1].real =a0 * A[1].real;
            B[2].real =a0 * A[2].real + a1 * A[1].real - A1 * B[1].real;
            B[3].real = a0 * A[3].real + a1 * A[2].real + a2 * A[1].real - A1 * B[2].real - A2 * B[1].real;
            for(int n = 4; n <= N; n++)
                B[n].real =a0 * A[n].real + a1 * A[n - 1].real + a2 * A[n - 2].real + a3 * A[n - 3].real - A1 * B[n - 1].real - A2 * B[n - 2].real - A3 * B[n - 3].real;

        }
        if(grade== 4)
        {
            double a0 = pow(wa, 4);
            double a1 = (double)4 * pow(wa, 4);
            double a2 = (double)6 * pow(wa, 4);
            double a3 = (double)4 * pow(wa, 4);
            double a4 = pow(wa, 4);
            double A0 = (double)16 + 20.904 * wa + 12.564 * pow(wa, 2) + 5.226 * pow(wa, 3) + pow(wa, 4);
            double A1 = ((double)-64 - 41.808 * wa) + 10.452 * pow(wa, 3) + (double)4 * pow(wa, 4);
            double A2 = ((double)96 - 25.128 * pow(wa, 2)) + (double)6 * pow(wa, 4);
            double A3 = (((double)-64 + 41.808 * wa) - 10.452 * pow(wa, 3)) + (double)4 * pow(wa, 4);
            double A4 = ((((double)16 - 20.904 * wa) + 12.564 * pow(wa, 2)) - 5.226 * pow(wa, 3)) + pow(wa, 4);
            a0 /= A0;
            a1 /= A0;
            a2 /= A0;
            a3 /= A0;
            a4 /= A0;
            A1 /= A0;
            A2 /= A0;
            A3 /= A0;
            A4 /= A0;
            B[1].real=a0 * A[1].real;
            B[2].real=(a0 * A[2].real + a1 * A[1].real) - A1 * B[1].real;
            B[3].real=(a0 * A[3].real + a1 * A[2].real + a2 * A[1].real) - A1 * B[2].real - A2 * B[1].real;
            B[4].real=(a0 * A[4].real + a1 * A[3].real + a2 * A[2].real + a3 * A[1].real) - A1 * B[3].real - A2 * B[2].real - A3 * B[1].real;
            for(int n = 5; n <= N; n++)
                B[n].real=(a0 * A[n].real + a1 * A[n - 1].real + a2 * A[n - 2].real + a3 * A[n - 3].real + a4 * A[n - 4].real) - A1 * B[n - 1].real - A2 * B[n - 2].real - A3 * B[n - 3].real - A4 * B[n - 4].real;

        }
        if(grade == 5)
        {
            double a0 = pow(wa, 5);
            double a1 = (double)5 * pow(wa, 5);
            double a2 = (double)10 * pow(wa, 5);
            double a3 = (double)10 * pow(wa, 5);
            double a4 = (double)5 * pow(wa, 5);
            double a5 = pow(wa, 5);
            double A0 = (double)32 + 51.776000000000003 * wa + 41.887999999999998 * pow(wa, 2) + 20.943999999999999 * pow(wa, 3) + 6.4720000000000004 * pow(wa, 4) + pow(wa, 5);
            double A1 = ((double)-160 - 155.328 * wa - 41.887999999999998 * pow(wa, 2)) + 20.943999999999999 * pow(wa, 3) + 19.416 * pow(wa, 4) + (double)5 * pow(wa, 5);
            double A2 = (((double)320 + 103.55200000000001 * wa) - 83.775999999999996 * pow(wa, 2) - 41.887999999999998 * pow(wa, 3)) + 12.944000000000001 * pow(wa, 4) + (double)10 * pow(wa, 5);
            double A3 = (((double)-320 + 103.55200000000001 * wa + 83.775999999999996 * pow(wa, 2)) - 41.887999999999998 * pow(wa, 3) - 12.944000000000001 * pow(wa, 4)) + (double)10 * pow(wa, 5);
            double A4 = ((((double)160 - 155.328 * wa) + 41.887999999999998 * pow(wa, 2) + 20.943999999999999 * pow(wa, 3)) - 19.416 * pow(wa, 4)) + (double)5 * pow(wa, 5);
            double A5 = (((((double)-32 + 51.776000000000003 * wa) - 41.887999999999998 * pow(wa, 2)) + 20.943999999999999 * pow(wa, 3)) - 6.4720000000000004 * pow(wa, 4)) + pow(wa, 5);
            a0 /= A0;
            a1 /= A0;
            a2 /= A0;
            a3 /= A0;
            a4 /= A0;
            a5 /= A0;
            A1 /= A0;
            A2 /= A0;
            A3 /= A0;
            A4 /= A0;
            A5 /= A0;
            B[1].real=a0 * A[1].real;
            B[2].real=(a0 * A[2].real + a1 * A[1].real) - A1 * B[1].real;
            B[3].real=(a0 * A[3].real + a1 * A[2].real + a2 * A[1].real) - A1 * B[2].real - A2 * B[1].real;
            B[4].real=(a0 * A[4].real + a1 * A[3].real + a2 * A[2].real + a3 * A[1].real) - A1 * B[3].real - A2 * B[2].real - A3 * B[1].real;
            B[5].real=(a0 * A[5].real + a1 * A[4].real + a2 * A[3].real + a3 * A[2].real + a4 * A[1].real) - A1 * B[4].real - A2 * B[3].real - A3 * B[2].real - A4 * B[1].real;
            for(int n = 6; n <= N; n++)
                B[n].real=(a0 * A[n].real + a1 * A[n - 1].real + a2 * A[n - 2].real + a3 * A[n - 3].real + a4 * A[n - 4].real + a5 * A[n - 5].real) - A1 * B[n - 1].real - A2 * B[n - 2].real - A3 * B[n - 3].real - A4 * B[n - 4].real - A5 * B[n - 5].real;

        }
}

void CWavefilterDlg::OnCancel() 
{
	/*
	dlgtest dlg;
	if(dlg.DoModal()==IDOK)
	{
	return;
	}
	*/
	CDialog::OnCancel();
}

void CWavefilterDlg::OnMouseMove(UINT nFlags, CPoint point) 
{
	// TODO: Add your message handler code here and/or call default
	int x=point.x;
	int y=point.y;
	if(x>=30&&x<=320&&y>140&&y<=280)
	{
		m_num1.Format("%d",x);
		m_num2.Format("%d",y);
	}
	UpdateData(false);
	CDialog::OnMouseMove(nFlags, point);
}

void CWavefilterDlg::TestHighPass()
{
	int i=1;
	if(m_wave.GetCurSel()==1)
		m_zhankong_msg.Format("%d%%",zhankong);
	m_point_msg.Format("%d",pointnum);
	m_period_msg.Format("%2gT",period/10.0);
	m_phase_msg.Format("%1gPI",phase*0.01);
	m_pointperperiod_msg.Format("每周期点数：%d",(int)pointnum*10/period);
	UpdateData(false);
	DrawCoordinate();
	Complex *x=new Complex[1024];
	Complex *X=new Complex[pointnum+1];
	int ftype=m_wave.GetCurSel();
	InitData(ftype,x);

	for(i=0;i<=pointnum;i++)
		X[i]=x[i];
	x1=x;
	DrawInPic1(pointnum,x1);



	int pointtemp=pointnum;
	int ffttime=0;
	while(pointtemp!=0)
	{
		pointtemp/=2;
		ffttime++;
	}
	ffttime--;
	if(pointnum>pow(2,ffttime))
	{
		ffttime++;
	}
	for(i=pointnum+1;i<=pow(2,ffttime);i++)
	{
		x[i].real=0;
		x[i].image=0;
	}

	FFT((int)pow(2,ffttime),ffttime,x);
	x2=x;
	DrawInPic2((int)pow(2,ffttime),x2);
	

	int start=3;
	double wd=(double)phase*0.01;

        wd *= PI;
		Complex *A=new Complex[1025];
		Complex *B=new Complex[1025];
		for(i=0;i<=pointnum;i++)
			A[i]=X[i];
	  for(i = 1; i <= 256; i++)
	 {
            B[i].real=0.0;
			B[i].image=0.0;
	 }
		int N;
		N=pointnum;
        double wa = (double)2 * tan(0.1*PI);
            double B0 = pow(wa, 3);
            double B1 = (double)3 * pow(wa, 3);
            double B2 = (double)3 * pow(wa, 3);
            double B3 = pow(wa, 3);
            double A0 = (double)8 + (double)8 * wa + (double)4 * pow(wa, 2) + pow(wa, 3);
            double A1 = ((double)-24 - (double)8 * wa) + (double)4 * pow(wa, 2) + (double)3 * pow(wa, 3);
            double A2 = ((double)24 - (double)8 * wa - (double)4 * pow(wa, 2)) + (double)3 * pow(wa, 3);
            double A3 = (((double)-8 + (double)8 * wa) - (double)4 * pow(wa, 2)) + pow(wa, 3);
            B0 /= A0;
            B1 /= A0;
            B2 /= A0;
            B3 /= A0;
            A1 /= A0;
            A2 /= A0;
            A3 /= A0;
			double a0,b0,b1,b2,b3,a1,a2,a3;
			double a=-cos((0.2*PI+0.6*PI)/2.0)/cos((0.2*PI-0.6*PI)/2.0);
			a0=B0-a*B1+a*a*B2-a*a*a*B3;
			a1=3*a*B0-B1*(1+2*a*a)+B2*(a*a*a+2*a)-3*a*a*B3;
			a2=3*a*a*B0-B1*(a*a*a+2*a)+B2*(1+2*a*a)-3*a*B3;
			a3=a*a*a*B0-B1*a*a+B2*a-B3;
			b0=(double)1-a*A1+a*a*A2-a*a*a*A3;
			b1=3*a-A1*(1+2*a*a)+A2*(a*a*a+2*a)-3*a*a*A3;
			b2=(a*a*a+2*a)+A2*(1+2*a*a)-3*a*A3;
			b3=a*a*a-A1*a*a+A2*a-A3;

			b2=3*a*a-A1*(a*a*a+2*a)+A2*(1+2*a*a)-3*a*A3;
            B[1].real =(a0 * A[1].real)/b0;
            B[2].real =(a0 * A[2].real + a1 * A[1].real - b1 * B[1].real)/b0;
            B[3].real = (a0 * A[3].real + a1 * A[2].real + a2 * A[1].real - b1 * B[2].real - b2 * B[1].real)/b0;
            for(int n = 4; n <= N; n++)
                B[n].real =(a0 * A[n].real + a1 * A[n - 1].real + a2 * A[n - 2].real + a3 * A[n - 3].real - b1 * B[n - 1].real - b2 * B[n - 2].real - b3 * B[n - 3].real)/b0;
	y1=B;
	DrawOutPic1(pointnum,B);
	FFT((int)pow(2,ffttime),ffttime,y1);
	y2=B;
	DrawOutPic2((int)pow(2,ffttime),y2);
}

void CWavefilterDlg::OnHighpass() 
{
	TestHighPass();
	
}