📄 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_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);//yyyyyyyyyyyy
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;
//输入参数. 主要是截止频率和衰减分贝(根据所需要的滤波器进行人为设定)
int N=10;//巴特沃思滤波器的阶数
double fs=50;//10*1000;//采样频率(Hz)
double fc=9;//3*1000;//2*1000;//高通3DB的截止频率(Hz)
//计算周期和截止处数字角频率
double T, wc;
T=1/fs;//采样周期(即采样时间间隔)
wc=2*PI*fc*T;
//计算c1
double aumigapc=1;//原型滤波器的截止频率. 此值总是取为常数1(rad/s)
double c1=aumigapc*tan(wc/2);
if (N>10)
{AfxMessageBox("滤波器的阶数要求太高(>10), 请调整设计要示, 重新设计!"); exit(0);}
astr.Format(" N = %d", N);
pDC->TextOut(20,10, astr);
//求数字高通滤波器的传递函数H(z)
double fenzi[2], fenmu[2];//p=f(z)表达式中的分子多项式(1 + Z^-1)
//分母多项式(1 - Z^-1)
double p[12][11], pSum[11];//H(z)分子多项式(1个)、分母中的多项式(1+N个)及它们的和.
//由于只有N<=10的滤波器有设计数据可查, 所以此处最多有12
//个多项式, 每个多项式最多有1+N<=11项
int Np[12], Nsum;//各式项式p[k]所含的项数及它们和所含的项数
double fenziPower[11], fenmuPower[11];//分子幂和分母幂. 最多1+10=11项
int Nfenzi, Nfenmu;//分子幂的项数和分母幂的项数
int k;
fenzi[0]=c1; fenzi[1]=c1;
fenmu[0]=1; fenmu[1]=-1;
polyPower(fenmu, 1, 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, 1, N, p[k], &Np[k]);
else if(k==N+1)//最后一个多项式p[N+1](=p(z)分母的N 次幂)
polyPower(fenmu, 1, N, p[k], &Np[k]);
else//中间各个多项式p[k](=p(z)分子的N-k+1次幂 * 分母的k-1次幂)
{
polyPower(fenzi, 1, N-k+1, fenziPower, &Nfenzi);//分子的N-k+1次幂
polyPower(fenmu, 1, k-1, fenmuPower, &Nfenmu);//分母的k-1次幂
polyMultiplyPoly(fenziPower, Nfenzi, fenmuPower, Nfenmu, p[k], &Np[k]);//两者相乘
if (Np[k] != 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] != N || Nsum != N)
{AfxMessageBox("多项式计算有误(阶次不符), 请检查!"); exit(0);}
astr.Format("%15.8f", p[0][0]/pSum[0]);
pDC->TextOut(20,50, astr);
for (k=0; k<=N; k++)//将H(z)的分母多项式系数写成串
{
astr.Format("%15.8f", pSum[k]/pSum[0]);
pDC->TextOut(20,80+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.002)
{
z1.a=0; z1.b=0;//H(jw)的分子
for (k=0; k<=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<=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));
//画高通段起始频率竖线
pDC->SetPixel(x0+wc*ratioX, y0-moduleHw*ratioY, RGB(0,0,255));
astr.Format("%.1fHz", fc+1);
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; k++)//将H(z)的分子多项式系数写成串
{astrPart.Format("%15.10f", p[0][k]/pSum[0]); astr+=astrPart;}
astr+=" //\n";
toFile.WriteString(astr);
astr=" 分母:";
for (k=0; k<=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<=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<=N; k++)//组合分母,即Σak*y[n-k]项
{
astrPart.Format("%+15.10f*y[n-%d]", -pSum[k]/pSum[0],k);
astr+=astrPart;
}
toFile.WriteString(astr);
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -