filterdesigndlg.cpp

IIR 滤波器设计 输入设计参数 直接生成数字滤波器

// filterDesignDlg.cpp : implementation file
//

#include "stdafx.h"
#include "filterDesign.h"
#include "filterDesignDlg.h"

#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

//yyyyyyyyyyyyyyyyyyyy 复数结构及复数运算 yyyyyyyyyyyyyyyyyyyyy
#include "math.h"
double PI = atan(1)*4;

struct complexNumber//定义复数结构z=a+bi
{
	double a;
	double b;
};

complexNumber complexAdd(complexNumber z1, complexNumber z2)//复数加法
{
	complexNumber z;
	z.a=z1.a + z2.a;
	z.b=z1.b + z2.b;
	return z;
}

complexNumber complexSub(complexNumber z1, complexNumber z2)//复数减法
{
	complexNumber z;
	z.a=z1.a - z2.a;
	z.b=z1.b - z2.b;
	return z;
}

complexNumber complexMultiply(complexNumber z1, complexNumber z2)//复数乘法
{
	complexNumber z;
	z.a=z1.a*z2.a - z1.b*z2.b;
	z.b=z1.a*z2.b + z1.b*z2.a;
	return z;
}

complexNumber complexGonge(complexNumber z1)//共轭复数
{
	complexNumber z;
	z.a=  z1.a;
	z.b= -z1.b;
	return z;
}

double complexModule2(complexNumber z)//求复数模的平方
{
	return (z.a*z.a + z.b*z.b);
}

complexNumber complexDivide(complexNumber z1, complexNumber z2)//复数除法
{
	complexNumber z=complexMultiply(z1,complexGonge(z2));
	double module=complexModule2(z2);

	z.a/=module;
	z.b/=module;
	return z;
}

complexNumber complexInvert(complexNumber z1)//复数倒法
{
	complexNumber z=complexGonge(z1);
	double module=complexModule2(z1);

	z.a/=module;
	z.b/=module;
	return z;
}
//yyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyy

//yyyyyyyyyyyyyyyy 多项式运算 yyyyyyyyyyyyyyyyyyyyyyyyyy
void polyAdd(double* a, int m, double* b, int n, double* p, int* np)//多项式加法
{	//a、b代表两个多项式序列的系数, m、n分别是它们的自变量最高/最低次数, 它们比多项式的项数少1
	//p代表乘积多项式的系数, N是其自变量最高/最低次数, 它们比多项式的项数少1

	int k;

	*np=max(m,n);//和的长度等于长序列的长度

	if (m>n)//如果a序列更长
	{
		for (k=0; k<=n; k++)
			*(p+k)=*(a+k) + *(b+k);//前半截取a、b序列之和
		for (k=n+1; k<=*np; k++)
			*(p+k)=*(a+k);//后半截取a序列
	}
	else
	{
		for (k=0; k<=m; k++)
			*(p+k)=*(a+k) + *(b+k);//前半截取a、b序列之和
		for (k=m+1; k<=*np; k++)
			*(p+k)=*(b+k);//后半截取b序列
	}
}

void polyMultiplyNumber(double* a, int m, double c, double* p, int* np)//多项式乘以系数c
{
	for (int k=0; k<=m; k++)
		*(p+k)=*(a+k) * c;//各项都乘以系数c
	*np=m;//项数不变
}

void polyMultiplyPoly(double* a, int m, double* b, int n, double* p, int* np)//多项式乘法
{
	int k;
	double aa[100], bb[100];

	*np=m+n;//积的长度等于两序列长度之和

	for (k=0; k<=*np; k++)//分别延长序列a、b成aa、bb, 使其长度都达到m+n, 延长部份补0
		{aa[k]=0;  bb[k]=0;}
	for (k=0; k<=m; k++)//将序列a复制进aa的前半截
		aa[k]=*(a+k);
	for (k=0; k<=n; k++)//将序列b复制进bb的前半截
		bb[k]=*(b+k);

	for (k=0;k<=*np;k++)//计算各次项的系数
	{
		double sum=0;
		for (int i=0;i<=k;i++)//第k次项的系数由k+1个a[i]*b[k-i]相加而成
			sum=sum + aa[i] * bb[k-i];
		*(p+k)=sum;
	}
}

void polyPower(double* a, int m, int k, double* p, int* np)//多项式的幂(a的k次幂)
{
	int i=1, N1;
	double p1[100];
	
	if (k<=0)
	{ AfxMessageBox("多项式的幂不可求, 因为幂指数k<=0 !"); exit(0);}

	polyMultiplyNumber(a, m, 1, p1, &N1);//p1=a

	do//循环累乘a(共乘k-1个a)
	{
		polyMultiplyNumber(p1, N1, 1, p, np);//p1→p

		i++;
		if (i>k) break;

		polyMultiplyPoly(p, *np, a, m, p1, &N1);//p1=p*a
	}while(1);
}
//yyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyyy

/////////////////////////////////////////////////////////////////////////////
// CFilterDesignDlg dialog

CFilterDesignDlg::CFilterDesignDlg(CWnd* pParent /*=NULL*/)
	: CDialog(CFilterDesignDlg::IDD, pParent)
{
	//{{AFX_DATA_INIT(CFilterDesignDlg)
		// NOTE: the ClassWizard will add member initialization here
	//}}AFX_DATA_INIT
	// Note that LoadIcon does not require a subsequent DestroyIcon in Win32
	m_hIcon = AfxGetApp()->LoadIcon(IDR_MAINFRAME);
}

void CFilterDesignDlg::DoDataExchange(CDataExchange* pDX)
{
	CDialog::DoDataExchange(pDX);
	//{{AFX_DATA_MAP(CFilterDesignDlg)
		// NOTE: the ClassWizard will add DDX and DDV calls here
	//}}AFX_DATA_MAP
}

BEGIN_MESSAGE_MAP(CFilterDesignDlg, CDialog)
	//{{AFX_MSG_MAP(CFilterDesignDlg)
	ON_WM_SYSCOMMAND()
	ON_WM_PAINT()
	ON_WM_QUERYDRAGICON()
	ON_BN_CLICKED(IDC_Button_Calculate, OnButtonCalculate)
	ON_WM_LBUTTONDBLCLK()
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CFilterDesignDlg message handlers

BOOL CFilterDesignDlg::OnInitDialog()
{
	CDialog::OnInitDialog();
	
	// TODO: Add extra initialization here
	ShowWindow(SW_MAXIMIZE);//yyyyyyyyyyyyyyyy
	
	return TRUE;  // return TRUE  unless you set the focus to a control
}

void CFilterDesignDlg::OnSysCommand(UINT nID, LPARAM lParam)
{
	CDialog::OnSysCommand(nID, lParam);
}

void CFilterDesignDlg::OnPaint() 
{
	CDialog::OnPaint();
}

// The system calls this to obtain the cursor to display while the user drags
//  the minimized window.
HCURSOR CFilterDesignDlg::OnQueryDragIcon()
{
	return (HCURSOR) m_hIcon;
}

void CFilterDesignDlg::OnLButtonDblClk(UINT nFlags, CPoint point) 
{
	// TODO: Add your message handler code here and/or call default
	exit(0);//yyyyyyyyyyyyyyyyyyyyy

	CDialog::OnLButtonDblClk(nFlags, point);
}

void CFilterDesignDlg::OnButtonCalculate() 
{
	CDC* pDC=GetDC();
	pDC->SetBkMode(TRANSPARENT);
	CString astr;
	int positionY=10;
	int x0=500, y0=300;
	pDC->MoveTo(x0-50, y0);
	pDC->LineTo(x0+500, y0);
	pDC->MoveTo(x0, y0-200);
	pDC->LineTo(x0, y0+200);

	//输入参数. 主要是截止频率和衰减分贝(根据所需要的滤波器进行人为设定)
	double fs, fc, fst;
	fs=50;//10*1000;//1;//采样频率(Hz)
	fc=3;//1*1000;//0.1;//低通段终止频率(Hz)
	fst=4;//1.5*1000;//0.15;//衰减段终止频率(Hz)
	
	double DBc, DBst;
	DBc=1;//低通段内允许的最大衰减(DB)
	DBst=10;//15;//衰减段必须达到的衰减程度(DB)

	//计算周期和角频率
	double T, wc, wst;
	T=1/fs;//采样周期(即采样时间间隔)
	wc=2*PI*fc*T;//低通段终止点的角频率
	wst=2*PI*fst*T;//衰减段终止点的角频率
	
	//计算原型滤波器的的阶数N和低通截止频率Ωc
	double n=(1/2.0)*log((pow(10,DBc/10)-1)/(pow(10,DBst/10)-1)) 
		           / log(wc/wst);//求滤波器的阶数n
	int N=(int)ceil(n);//对n向上取整
	if (N%2 != 0)//如果N是奇数, 则取比之略大的偶数
		N=N+1;

				astr.Format("         N = %d", N);
				pDC->TextOut(20, positionY, astr);
				positionY+=20;

	if (N>=20)
		{MessageBox("滤波器的阶数太高! 请调整设计要求,重新设计."); exit(0);}

	double aumigac=(wc/T)/pow(pow(10,DBc/10)-1, 1.0/(2*N));//Ωc
    
	//计算模拟滤波器的极点s[k](=a+bi)
	complexNumber s[20];
	int k;
    for (k=1; k<=N; k++)
	{
        s[k].a = aumigac * cos((1.0/2 + (2*k-1)/(2.0*N)) * PI);
        s[k].b = aumigac * sin((1.0/2 + (2*k-1)/(2.0*N)) * PI);
        
				astr.Format("%.9f + %.9fj", s[k].a, s[k].b);
				pDC->TextOut(20, positionY, astr);
				positionY+=20;
    }
	positionY+=30;

	//模拟传递函数Ha(s)分式的分子
	double 	K0=pow(aumigac,N);
    astr.Format("%f", K0);
	pDC->TextOut(20, positionY, astr);
	positionY+=20;
	
	//计算两共轭极点所对应部份分式的分母, 即(s-sk)×(s-sk*)
	double c[20][6];//二次分式的系数. 从分子到分母、从左到右依次为c[1]、c[2]、…、c[5]
	for (k=1; k<=N/2; k++)
	{
		c[k][4]=-2*s[k].a;//[4]代表分母的第二项
		c[k][5]=s[k].a*s[k].a + s[k].b*s[k].b;//[5]代表分母的最后一项

		astr.Format("s^2 + %.4fs + %.4f", c[k][4], c[k][5]);	
		pDC->TextOut(20, positionY, astr);
		positionY+=20;
	}
	positionY+=30;

	//计算Ha(s)展开成部份分式后各分式的分子A[k](=a+bi)
	complexNumber A[20];
	for (k=1; k<=N; k++)
	{
		complexNumber z;
		z.a=1;	z.b=0;
		for (int j=1; j<=N; j++)
			if (j!=k)
				z=complexMultiply(z, complexSub(s[k],s[j]));
		
		A[k]=complexInvert(z);
		A[k].a*=K0*T;
		A[k].b*=K0*T;
	}

	//计算数字传递函数H(z). 两个共轭极点组成一个二次多项式分式
	for (k=1; k<=N/2; k++)
	{
		c[k][1]=2*A[k].a;
		c[k][2]=-(2*A[k].a*exp(s[k].a*T)*cos(s[k].b*T) + 2*A[k].b*exp(s[k].a*T)*sin(s[k].b*T));
		c[k][3]=1;
		c[k][4]=-2*exp(s[k].a*T)*cos(s[k].b*T);
		c[k][5]=exp(2*s[k].a*T);

		astr.Format("%7.4f %12.4f %12.0f %16.4f %12.4f"
			      , c[k][1], c[k][2], c[k][3], c[k][4], c[k][5]);
		pDC->TextOut(20, positionY, astr);
		positionY+=20;
	} 

	//将H(z)中的各个多项式分式合并
	double fenzi[20][2], fenmu[20][3]; //分式多项式的分子、分母, 分子肯定是1次, 分母肯定是2次
	                                   //假设总共有不超过20个分式多项式
	double fenziSum[41], fenmuSum[41];//各分式通分并合并同类项后所得的分子、分母多项式. 
	                                   //20个二次多项式最多合并成40阶多项式
	int Nfenzi, Nfenmu;//分子分母多项式的项数

	double polytemp[20][41], polytemp2[41];
	int ntemp[20], ntemp2, i;

	for (k=1; k<=N/2; k++)//构造分子分母多项式
	{
		fenzi[k][0]=c[k][1];	fenzi[k][1]=c[k][2];//c[1]~c[5]中的前两个是分子系数, 后三个是分母系数
		fenmu[k][0]=c[k][3];	fenmu[k][1]=c[k][4];	fenmu[k][2]=c[k][5];
	}

	for (k=1; k<=N/2; k++)//遍历各个分式, 取其分子. 然后乘以其余分式的分母, 构成分子的一项. 
	{
		polyMultiplyNumber(fenzi[k], 1, 1, polytemp[k], &ntemp[k]);//取一个分式的分子

		for (i=1; i<=N/2; i++)//遍历各个分式, 取其分母, 以与上面的分子相乘
			if (i!=k)
			{
				polyMultiplyPoly(fenmu[i], 2, polytemp[k], ntemp[k], polytemp2, &ntemp2);//累乘一个分母
				polyMultiplyNumber(polytemp2, ntemp2, 1, polytemp[k], &ntemp[k]);//结果转回给polytemp[k]

			}
	}

				for (k=1; k<=N/2; k++)//遍历各个分式, 取其分子. 然后乘以其余分式的分母, 构成分子的一项. 
					if (ntemp[k] != N-1)
						{AfxMessageBox("多项式计算有误(阶次不符), 请检查!");  exit(0);}

	polyMultiplyNumber(polytemp[1], ntemp[1], 1, fenziSum, &Nfenzi);//把各个分子项累加起来, 得总的分子
	for (k=2; k<=N/2; k++)
		polyAdd(fenziSum, Nfenzi, polytemp[k], ntemp[k], fenziSum, &Nfenzi);

	polyMultiplyNumber(fenmu[1], 2, 1, fenmuSum, &Nfenmu);//累乘各个分母累, 得总的分母
	for (k=2; k<=N/2; k++)
	{
		polyMultiplyPoly(fenmuSum, Nfenmu, fenmu[k], 2, polytemp2, &ntemp2);//累乘一个分母
		polyMultiplyNumber(polytemp2, ntemp2, 1, fenmuSum, &Nfenmu);//结果转回给fenmuSum
	}

				if (Nfenzi!=N-1 || Nfenmu!=N)
					{AfxMessageBox("多项式计算有误(阶次不符), 请检查!");  exit(0);}

	//将H(z)用z=exp(jw)化成H(jw),以便画频率响应图
	double w, moduleHw;
	complexNumber z1, z2, Hw;
	double ratioX=200, ratioY=100;

	for (w=0; w<=PI; w+=0.001)
	{
		z1.a=0;		z1.b=0;//H(jw)的分子
		for (k=0; k<=N-1; k++)//将分子按exp(-jw)=cos(w)-jsin(w)展开后, 进行实部和虚部的合并
		{
			z1.a+= fenziSum[k]*cos(k*w);
			z1.b+=-fenziSum[k]*sin(k*w);
		}
		
		z2.a=0;		z2.b=0;//H(jw)的分母
		for (k=0; k<=N; k++)//将分母按exp(-jw)=cos(w)-jsin(w)展开后, 进行实部和虚部的合并
		{
			z2.a+= fenmuSum[k]*cos(k*w);
			z2.b+=-fenmuSum[k]*sin(k*w);
		}

		Hw=complexDivide(z1,z2);//分子除分母,得H(jw)

		moduleHw=sqrt(complexModule2(Hw));//求|H(jw)|

		if (moduleHw>1.1)
			{MessageBox("幅度大于1!"); exit(0);}

		//画幅度曲线|H(jw)|=f(w)
		pDC->SetPixel(x0+w*ratioX, y0-moduleHw*ratioY, RGB(255,0,0));

		//将|H(jw)|对w画图, 得幅度曲线
		pDC->SetPixel(x0+w*ratioX, y0-moduleHw*ratioY, RGB(255,0,0));

		//画截止频率竖线
		pDC->SetPixel(x0+wc*ratioX, y0-moduleHw*ratioY, RGB(0,0,255));

		//标截止频率
		astr.Format("%.1fHz", fc);
		pDC->TextOut(x0+wc*ratioX-20, y0+5, astr);
	}

	//将设计所得滤波器的参数输出成数据文件
	CStdioFile toFile;
	CString fileName="E:\\VC程序集\\工作\\滤波器设计\\滤波器设计.txt";

	if(!toFile.Open(fileName, CFile::modeCreate|CFile::modeNoTruncate|CFile::modeWrite))//打开输出文件. 
		    //如该文件不存在, 则modeCreate使之被创建; 如果已经存在, 则modeNoTruncate使其内容不被清空
		{MessageBox("文件打开失败，请检查!"); exit(0);}

	toFile.SeekToEnd();//定位到文件尾继续写
	toFile.WriteString("\n\n---------------------------------------------------------------------------------------------------------\n");
	toFile.WriteString("低通滤波器: \n");

	CString astrPart;
	astr="    分子:";
	for (k=0; k<=N-1; k++)//将H(z)的分子多项式系数写成串
		{astrPart.Format("%15.10f", fenziSum[k]);	astr+=astrPart;}
	astr+="  //\n";
	toFile.WriteString(astr);
	
	astr="    分母:";
	for (k=0; k<=N; k++)//将H(z)的分母多项式系数写成串
		{astrPart.Format("%15.10f", fenmuSum[k]);	astr+=astrPart;}
	astr+="  //\n";
	toFile.WriteString(astr);

	//将传递函数的分子分母组合成差分方程的形式
	toFile.WriteString("    差分方程: \n");
	astr="        y[n]=";
	for (k=0; k<=N-1; k++)//组合分子,即Σbr*x[n-r]项
	{
		astrPart.Format("%+15.10f*x[n-%d]", fenziSum[k], k);
		astr+=astrPart;
	}
	
	for (k=1; k<=N; k++)//组合分母,即Σak*y[n-k]项
	{
		astrPart.Format("%+15.10f*y[n-%d]", -fenmuSum[k],k);
		astr+=astrPart;
	}
	
	toFile.WriteString(astr);
}