📄 filterdesigndlg.cpp
字号:
// 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_BN_CLICKED(IDC_Button_Test, OnButtonTest)
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);//yyyyyyyyyyyyyyy
CDialog::OnLButtonDblClk(nFlags, point);
}
void CFilterDesignDlg::OnButtonCalculate()
{
CDC* pDC=GetDC();
pDC->SetBkMode(TRANSPARENT);
CString astr;
int positionY=10;
int x0=200, y0=500;
pDC->MoveTo(x0-50, y0);
pDC->LineTo(x0+700, y0);
pDC->MoveTo(x0, y0-200);
pDC->LineTo(x0, y0+200);
//巴特沃思多项式的系数
double a[11][11];
a[1][0]=1; a[1][1]=1;
a[2][0]=1; a[2][1]=1.4142136; a[2][2]=1;
a[3][0]=1; a[3][1]=2.0000000; a[3][2]=2.0000000; a[3][3]=1;
a[4][0]=1; a[4][1]=2.6131259; a[4][2]=3.4142136; a[4][3]=2.6131259; a[4][4]=1;
a[5][0]=1; a[5][1]=3.2360680; a[5][2]=5.2360680; a[5][3]=5.2360680; a[5][4]=3.2360680; a[5][5]=1;
a[6][0]=1; a[6][1]=3.8637033; a[6][2]=7.4641016; a[6][3]=9.1416202; a[6][4]=7.4641016; a[6][5]=3.8637033; a[6][6]=1;
a[7][0]=1; a[7][1]=4.4939592; a[7][2]=10.0978347; a[7][3]=14.5917939; a[7][4]=14.5917939; a[7][5]=10.0978347; a[7][6]=4.4939592; a[7][7]=1;
a[8][0]=1; a[8][1]=5.1258909; a[8][2]=13.1370712; a[8][3]=21.8461510; a[8][4]=25.6883559; a[8][5]=21.8461510; a[8][6]=13.1370712; a[8][7]=5.1258309; a[8][8]=1;
a[9][0]=1; a[9][1]=5.7587705; a[9][2]=16.5817187; a[9][3]=31.1634375; a[9][4]=41.9863857; a[9][5]=41.9863857; a[9][6]=31.1634375; a[9][7]=16.5817187; a[9][8]=5.7587705; a[9][9]=1;
a[10][0]=1; a[10][1]=6.3924532; a[10][2]=20.4317291; a[10][3]=42.8020611; a[10][4]=64.8823963; a[10][5]=74.2334292; a[10][6]=64.8823963; a[10][7]=42.8020611; a[10][8]=20.4317291; a[10][9]=6.3924532; a[10][10]=1;
//输入参数. 主要是截止频率和衰减分贝(根据所需要的滤波器进行人为设定)
double fs, f0, f1, f3, f4, DB1, DB2;
fs=50;//1000;//2*1000;//采样频率(Hz)
f0=2;//100;//200;//带通滤波器的低阻段终止频率(Hz)
f1=3;//200;//300;//带通滤波器的带通段起始频率(Hz)
f3=10;//250;//400;//带通滤波器的带通段终止频率(Hz)
f4=12;//400;//500;//带通滤波器的高阻段起始频率(Hz)
DB1=1;//3;//带通段内允许的最大衰减(DB)
DB2=10;//18;//20;//高、低阻段内必须达到的衰减程度(DB)
//计算周期和模拟角频率
double T, D, aumigapc, aumiga3, aumiga1, E;
T=1/fs;//采样周期(即采样时间间隔)
aumiga3=2*PI*f3;//Ω3
aumiga1=2*PI*f1;//Ω1
aumigapc=1;//Ωpc. 它一般取为常数1(rad/s)
D=aumigapc/tan((aumiga3-aumiga1)/2*T);//求D
E=2*cos((aumiga3+aumiga1)/2*T) / cos((aumiga3-aumiga1)/2*T);//求E
//计算并定Ωps
double aumigaps1, aumigaps2, aumigaps;
aumigaps1=-D*(E/2-cos(2*PI*f0*T)) / sin(2*PI*f0*T);//根据低阻段终止频率求Ωps
aumigaps2=D*(E/2-cos(2*PI*f4*T)) / sin(2*PI*f4*T);//根据高阻段起始频率求Ωps
aumigaps=min(aumigaps1, aumigaps2);//取两个Ωps中的小者作为终定的Ωps, 这样可以使滤波器更陡
//确定原型滤波器阶数N
double aumigap=D*(E/2-cos(2*PI*f3*T)) / sin(2*PI*f3*T);//根据带通段终止频率求Ωp
double n=log10((pow(10, DB2/10)-1)/(pow(10, DB1/10)-1))/(2*log10(aumigaps/aumigap));
int N=(int)ceil(n);
astr.Format(" N = %9.6f", n);
pDC->TextOut(20,20, astr);
if (N>10)
{AfxMessageBox("带通滤波器的阶数要求太高(>10), 请调整设计要示, 重新设计!"); exit(0);}
//求数字带通滤波器的传递函数H(z)
double fenzi[3], fenmu[3];//p=f(z)表达式中的分子多项式D(1 - EZ^-1 + Z^-2)
//分母多项式(1 - Z^-2)
double p[12][21], pSum[21];//H(z)分子多项式(1个)、分母中的多项式(1+N个)及它们的和.
//由于只有N<=10的滤波器有设计数据可查, 所以此处最多有12
//个多项式, 每个多项式最多有1+2N<=21项
int Np[12], Nsum;//各式项式p[k]所含的项数及它们和所含的项数
double fenziPower[21], fenmuPower[21];//分子幂和分母幂. 最多1+2*10=21项
int Nfenzi, Nfenmu;//分子幂的项数和分母幂的项数
int k;
fenzi[0]=D; fenzi[1]=-D*E; fenzi[2]=D;
fenmu[0]=1; fenmu[1]=0; fenmu[2]=-1;
polyPower(fenmu, 2, N, p[0], &Np[0]);//计算H(z)分子多项式(=p(z)分母的N次幂)
for (k=1; k<=N+1; k++)//计算H(z)分母中的N+1个多项式.(按H(p)中按p的降序排列)
{
if (k==1)//第一个多项式p[1](=p(z)分子的N 次幂)
polyPower(fenzi, 2, N, p[k], &Np[k]);
else if(k==N+1)//最后一个多项式p[N+1](=p(z)分母的N 次幂)
polyPower(fenmu, 2, N, p[k], &Np[k]);
else//中间各个多项式p[k](=p(z)分子的N-k+1次幂 * 分母的k-1次幂)
{
polyPower(fenzi, 2, N-k+1, fenziPower, &Nfenzi);//分子的N-k+1次幂
polyPower(fenmu, 2, k-1, fenmuPower, &Nfenmu);//分母的k-1次幂
polyMultiplyPoly(fenziPower, Nfenzi, fenmuPower, Nfenmu, p[k], &Np[k]);//两者相乘
if (Np[k] != 2*N)
{AfxMessageBox("多项式计算有误(阶次不符), 请检查!"); exit(0);}
}
polyMultiplyNumber(p[k], Np[k], a[N][k-1], p[k], &Np[k]);//各多项式乘以对应的巴特沃思多项式系数
}
polyMultiplyNumber(p[1], Np[1], 1, pSum, &Nsum);//把各多项式加起来(同类项合并)
for (k=2; k<=N+1; k++)
polyAdd(pSum, Nsum, p[k], Np[k], pSum, &Nsum);
if (Np[0] != 2*N || Nsum != 2*N)
{AfxMessageBox("多项式计算有误(阶次不符), 请检查!"); exit(0);}
astr.Format("%15.8f", p[0][0]/pSum[0]);
pDC->TextOut(20,70, astr);
for (k=0; k<=2*N; k++)//将H(z)的分母多项式系数写成串
{
astr.Format("%15.8f", pSum[k]/pSum[0]);
pDC->TextOut(20,100+k*30, astr);
}
//将H(z)用z=exp(jw)化成H(jw),以便画频率响应图
double w, moduleHw;
complexNumber z1, z2, Hw;
double ratioX=150, ratioY=100;
for (w=0; w<=PI; w+=0.001)
{
z1.a=0; z1.b=0;//H(jw)的分子
for (k=0; k<=2*N; k++)//将分子按exp(-jw)=cos(w)-jsin(w)展开后, 进行实部和虚部的合并
{
z1.a+= p[0][k]*cos(k*w);
z1.b+=-p[0][k]*sin(k*w);
}
z2.a=0; z2.b=0;//H(jw)的分母
for (k=0; k<=2*N; k++)//将分母按exp(-jw)=cos(w)-jsin(w)展开后, 进行实部和虚部的合并
{
z2.a+= pSum[k]*cos(k*w);
z2.b+=-pSum[k]*sin(k*w);
}
Hw=complexDivide(z1,z2);//分子除分母,得H(jw)
moduleHw=sqrt(complexModule2(Hw));//求|H(jw)|
//画幅度曲线|H(jw)|=f(w)
pDC->SetPixel(x0+w*ratioX, y0-moduleHw*ratioY, RGB(255,0,0));
//画通带起止频率竖线
double wc1=(2*PI*f1*T);//标起始频率
pDC->SetPixel(x0+wc1*ratioX, y0-moduleHw*ratioY, RGB(0,0,255));
astr.Format("%.1fHz", f1);
pDC->TextOut(x0+wc1*ratioX-20, y0+5, astr);
double wc3=(2*PI*f3*T);//标终止频率
pDC->SetPixel(x0+wc3*ratioX, y0-moduleHw*ratioY, RGB(0,0,255));
astr.Format("%.1fHz", f3);
pDC->TextOut(x0+wc3*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<=2*N; k++)//将H(z)的分子多项式系数写成串
{astrPart.Format("%15.10f", p[0][k]/pSum[0]); astr+=astrPart;}
astr+=" //\n";
toFile.WriteString(astr);
astr=" 分母:";
for (k=0; k<=2*N; k++)//将H(z)的分母多项式系数写成串
{astrPart.Format("%15.10f", pSum[k]/pSum[0]); astr+=astrPart;}
astr+=" //\n";
toFile.WriteString(astr);
//将传递函数的分子分母组合成差分方程的形式
toFile.WriteString(" 差分方程: \n");
astr=" y[n]=";
for (k=0; k<=2*N; k++)//组合分子,即Σbr*x[n-r]项
{
astrPart.Format("%+15.10f*x[n-%d]", p[0][k]/pSum[0], k);
astr+=astrPart;
}
for (k=1; k<=2*N; k++)//组合分母,即Σak*y[n-k]项
{
astrPart.Format("%+15.10f*y[n-%d]", -pSum[k]/pSum[0],k);
astr+=astrPart;
}
toFile.WriteString(astr);
}
void CFilterDesignDlg::OnButtonTest()
{
CDC* pDC=GetDC();
pDC->SetBkMode(TRANSPARENT);
CString astr;
double a[100]={1,1}, b[100]={1,2}, p[100];
int np;
polyAdd(a, 1, b, 1, p, &np);
astr.Format("两多项式之和为:%10.0f %10.0f", *(p+0), *(p+1));
pDC->TextOut(20,20,astr);
polyMultiplyNumber(b, 1, 10, p, &np);
astr.Format("两多项式之和为:%10.0f %10.0f", *(p+0), *(p+1));
pDC->TextOut(20,50,astr);
polyMultiplyPoly(a, 1, b, 1, p, &np);
astr.Format("两多项式之和为:%10.0f %10.0f %10.0f", *(p+0), *(p+1), *(p+2));
pDC->TextOut(20,80,astr);
polyPower(a, 1, 3, p, &np);
astr.Format("两多项式之和为:%10.0f %10.0f %10.0f %10.0f", *(p+0), *(p+1), *(p+2), *(p+3), *(p+4));
pDC->TextOut(20,110,astr);
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -