📄 filterdesign.~cpp
字号:
H[0]=1;
H[1]=-2*(pow(OMEGA,2)+pow(Wi,2)-1)/(1-2*OMEGA*cos(Angle)+pow(OMEGA,2)+pow(Wi,2));
H[2]=(1+2*OMEGA*cos(Angle)+pow(OMEGA,2)+pow(Wi,2))/
(1-2*OMEGA*cos(Angle)+pow(OMEGA,2)+pow(Wi,2));
H_Length=3;
X_Length=Conv(H_Length,H,X_Length,X,pDen);
for (int j=0;j<X_Length;j++)
X[j]=pDen[j];
}
}
void ChebyOne::Num_Cheby_Highpass()
{
int K=order/2;
double X[20];
double H[3];
double Angle;
double Wi;
double Wpass=cos(PI*heightFre/fs)/sin(PI*heightFre/fs);
int X_Length,H_Length;
if (order%2==0)
{
X[0]=sqrt(1.0/pow(10,Ap/10));
X_Length=1;
}
else
{
X[0]=OMEGA/(OMEGA+1);
X[1]=-X[0];
X_Length=2;
pNum[0] = X[0];
pNum[1] = X[1];
}
for (int i=1;i<=K;i++)
{
Angle=PI*(order-1+2*i)/(2*order);
Wi=Wpass*sin(Angle);
H[0]=(pow(OMEGA,2)+pow(Wi,2))/(1-2*OMEGA*cos(Angle)+pow(OMEGA,2)+pow(Wi,2));
H[1]=-2*H[0];
H[2]=H[0];
H_Length=3;
X_Length=Conv(H_Length,H,X_Length,X,pNum);
for (int j=0;j<X_Length;j++)
X[j]=pNum[j];
}
}
/*!----------------------------------------------------------------------------
* \param passType 通带类型， 高通 低通
* \param fs 采样率
* \param lowFre 过渡带的低频
* \param heightFre 过渡带的高频
* \param aPass 最大通带衰减
* \param aStop 最小止带衰减
*
* ----------------------------------------------------------------------------*/
ChebyTwo::ChebyTwo(PassType passType, double fs, double lowFre, double heightFre,
double aPass, double aStop)
: IIRFilterDesign(passType, fs)
{
Init(lowFre, heightFre, aPass, aStop);
GenerateCoe();
}
/*!----------------------------------------------------------------------------
* \param passtype 通带类型， 高通 低通
* \param fs 采样率
* \param lowfre 过渡带的低频
* \param heightfre 过渡带的高频
* \param apass 使用默认 1 db
* \param astop 使用默认 60 db
*
* \warning apass astop 使用默认值
*
* ----------------------------------------------------------------------------*/
ChebyTwo::ChebyTwo(PassType passType, double fs, double lowFre, double heightFre)
: IIRFilterDesign(passType, fs)
{
Init(lowFre, heightFre, 1, 60);
GenerateCoe();
}
ChebyTwo::~ChebyTwo()
{
ReleaseArray();
}
void ChebyTwo::Init(double lowFre, double heightFre, double aPass, double aStop)
{
this->lowFre = lowFre;
this->heightFre = heightFre;
Ap = aPass;
As = aStop;
}
void ChebyTwo::GenerateCoe()
{
switch(passType)
{
case LOW_PASS:
GenerateLowCoe();
break;
case HEIGHT_PASS:
GenerateHeightCoe();
break;
default:
break;
}
}
/*!----------------------------------------------------------------------------
* \brief 生成低通系数
*
* ----------------------------------------------------------------------------*/
void ChebyTwo::GenerateLowCoe()
{
Cheby2_Lowpass_N();
AllocateArray();
Den_Cheby2_Lowpass();
Num_Cheby2_Lowpass();
}
/*!----------------------------------------------------------------------------
* \brief 生成高通系数
*
* ----------------------------------------------------------------------------*/
void ChebyTwo::GenerateHeightCoe()
{
Cheby2_Highpass_N();
AllocateArray();
Den_Cheby2_Highpass();
Num_Cheby2_Highpass();
}
/*!----------------------------------------------------------------------------
* \brief
* 幅频函数
*
* Chebyshev II 的幅频响应函数
* \f[
* |H(\Omega)|^2 = \frac{C_{N}^{\phantom{N}{2}}(\frac{\Omega_{stop}}{\Omega})}
* {C_{N}^{\phantom{N}{2}}(\frac{\Omega_{stop}}{\Omega}) + \varepsilon_{stop}
* ^{\phantom{stop}{2}}}
*
* \f]
*
* \f$C_{N} \f$ 为切比雪夫多项式， 见IIRFilterDesign::CN(int N,double x);
*
* \param freq 要求的频率点
*
* \return 幅频
*
* ----------------------------------------------------------------------------*/
double ChebyTwo::MagResponse(double freq)
{
double range = 0.0;
double freqW = 0.0;
// double wpass = 0.0;
double wstop = 0.0;
// double Wpass = 0.0;
double Wstop = 0.0;
// double Epass = 0.0;
double Estop = 0.0;
//TODO: 这里的计算应该和计算阶数的前面部分合起来提到上一层。 做个 xxx(fstop, fpass)的函数;
if(passType == LOW_PASS)
{
// wpass=2*PI*lowFre/fs;
wstop=2*PI*heightFre/fs;
// Wpass=tan(wpass/2);
Wstop=tan(wstop/2);
freqW = tan((PI*freq)/fs);
}
else
{
// wpass=2*PI*heightFre/fs;
wstop=2*PI*lowFre/fs;
// Wpass=cos(wpass/2)/sin(wpass/2);
Wstop=cos(wstop/2)/sin(wstop/2);
freqW = cos((PI*freq)/fs)/sin((PI*freq)/fs);
}
// Epass= sqrt(pow(10,Ap/10)-1);
Estop= sqrt(pow(10,As/10)-1);
//切比雪夫多项式的平方
double CN_2 = pow(CN(order, Wstop/freqW), 2);
range = CN_2/(CN_2 + pow(Estop, 2));
//range = sqrt(CN_2/(CN_2 + pow(Estop, 2)));
return range;
}
/*!----------------------------------------------------------------------------
* \brief 计算阶数
*
* a 和 \f$ \Omega_{0} \f$ 由下面公式算出
* \f[
* a = \frac{1}{N}\ln(\varepsilon_{stop} +
* \sqrt{\varepsilon_{stop}^{\phantom{stop}{2}} + 1})
* \f]
* \f[
* \Omega_{0} = \frac{\Omega_{stop}}{\sinh a}
* \f]
*
* 其它参考 ChebyOne::Cheby_Lowpass_N()
*
* ----------------------------------------------------------------------------*/
void ChebyTwo::Cheby2_Lowpass_N()
{
double wpass;
double wstop;
double Wpass;
double Wstop;
double Epass;
double Estop;
double e,w;
double a;
wpass=2*PI*lowFre/fs;
wstop=2*PI*heightFre/fs;
Wpass=tan(wpass/2);
Wstop=tan(wstop/2);
Epass= sqrt(pow(10,Ap/10)-1);
Estop= sqrt(pow(10,As/10)-1);
e=Estop/Epass;
w=Wstop/Wpass;
order=int(log(e+sqrt(pow(e,2)-1))/log(w+sqrt(pow(w,2)-1)))+1;
if(order > 16)
{
throw std::logic_error("设计的滤波器不稳定！请重试！");
}
a=log(Estop+sqrt(pow(Estop,2)+1))/order;
OMEGA=Wstop/sinh(a);
}
/*!----------------------------------------------------------------------------
* \brief Chebyshev II 低通分母
*
* 先根据阶数的奇偶进行初始化:
* 如果阶数为奇：
* \f[
* H_{0}(z) = \left\{ \begin{array}{ll} $1$& \textrm{当~N = 2K}\\
* \frac{G_{0}(1 + z^{-1})}{1 + a_{01}z^{-1}} & \textrm{当~N = 2K + 1}\\
* \end{array} \right.
* \f]
* 其中：
* \f[
* G_{0} = \frac{\Omega_{0}}{\Omega{0} + 1} \quad,\quad
* a_{01} = \frac{\Omega_{0} - 1}{\Omega{0} + 1}
* \f]
*
* 如果为偶，二阶节公式（求分母部分）：
* \f[
* H_{i}(z) = \frac{G_{i}(1+b_{i1}z^{-1}+z^{-2})}{1+a_{i1}z^{-1}+a_{i2}z^{-2}}
* \quad,\quad i = 1,2,\ldots,K
* \f]
*
* 变量H[1]:
* \f[
* a_{i1} = \frac{2(1 - \Omega_{0}^{\phantom{0}{-2}} -
* \Omega_{i}^{\phantom{i}{-2}})}{1 - 2\Omega_{0}^{\phantom{0}{-1}}\cos\theta_{i} +
* \Omega_{0}^{\phantom{0}{-2}} + \Omega_{i}^{\phantom{i}{-2}}}
* \f]
* 变量H[2]:
* \f[
* a_{i2} = \frac{1 + 2\Omega_{0}^{\phantom{0}{-1}}\cos\theta_{i} +
* \Omega_{0}^{\phantom{0}{-2}} + \Omega_{i}^{\phantom{i}{-2}}}{1 -
* 2\Omega_{0}^{\phantom{0}{-1}}\cos\theta_{i} +
* \Omega_{0}^{\phantom{0}{-2}} + \Omega_{i}^{\phantom{i}{-2}}}
* \f]
* 变量Wi:
* \f[
* \Omega_{i} = \frac{\Omega_{stop}}{\sin \theta_{i}} \quad,\quad
* i = 1,2,\ldots,K
* \f]
* 最后卷积求出分母系数
*
* ----------------------------------------------------------------------------*/
void ChebyTwo::Den_Cheby2_Lowpass()
{
int K=order/2;
double X[20];
double H[3];
double Angle;
double Wi;
double Wstop=tan(PI*heightFre/fs);
int X_Length,H_Length;
if (order%2==0)
{
X[0]=1;
X_Length=1;
}
else
{
X[0]=1;
X[1]=(OMEGA-1)/(OMEGA+1);
X_Length=2;
pDen[0] = X[0];
pDen[1] = X[1];
}
for (int i=1;i<=K;i++)
{
Angle=PI*(order-1+2*i)/(2*order);
Wi=Wstop/sin(Angle);
H[0]=1;
H[1]=2*(1-pow(OMEGA,-2)-pow(Wi,-2))/(1-2*cos(Angle)/OMEGA+pow(OMEGA,-2)+pow(Wi,-2));
H[2]=(1+2*cos(Angle)/OMEGA+pow(OMEGA,-2)+pow(Wi,-2))/
(1-2*cos(Angle)/OMEGA+pow(OMEGA,-2)+pow(Wi,-2));
H_Length=3;
X_Length=Conv(H_Length,H,X_Length,X,pDen);
for (int j=0;j<X_Length;j++)
X[j]=pDen[j];
}
}
/*!----------------------------------------------------------------------------
* \brief 生成低通分子
*
* 参考 ChebyTwo::Den_Cheby2_Lowpass()
*
* 分子中：
* \f[
* G_{i} = \frac{1 + \Omega_{i}^{\phantom{i}{-2}}}{1 -
* 2\Omega_{0}^{\phantom{0}{-1}}\cos\theta_{i} +
* \Omega_{0}^{\phantom{0}{-2}} + \Omega_{i}^{\phantom{i}{-2}}}
* \f]
* \f[
* b_{i1} = 2\frac{1 - \Omega_{i}^{\phantom{i}{-2}}}{1 + \Omega_{i}^{\phantom{i}{-2}}}
* \f]
*
* ----------------------------------------------------------------------------*/
void ChebyTwo::Num_Cheby2_Lowpass()
{
int K=order/2;
double X[20];
double H[3];
double Angle;
double Wi;
double Wstop=tan(PI*heightFre/fs);
int X_Length,H_Length;
if (order%2==0)
{
X[0]=1;
X_Length=1;
}
else
{
X[0]=OMEGA/(OMEGA+1);
X[1]=X[0];
X_Length=2;
pNum[0] = X[0];
pNum[1] = X[1];
}
for (int i=1;i<=K;i++)
{
Angle=PI*(order-1+2*i)/(2*order);
Wi=Wstop/sin(Angle);
H[0]=(1+pow(Wi,-2))/(1-2*cos(Angle)/OMEGA+pow(OMEGA,-2)+pow(Wi,-2));
H[1]=H[0]*2*(1-pow(Wi,-2))/(1+pow(Wi,-2));
H[2]=H[0];
H_Length=3;
X_Length=Conv(H_Length,H,X_Length,X,pNum);
for (int j=0;j<=X_Length;j++)
X[j]=pNum[j];
}
}
void ChebyTwo::Cheby2_Highpass_N()
{
double wpass;
double wstop;
double Wpass;
double Wstop;
double Epass;
double Estop;
double e,w;
double a;
wpass=2*PI*heightFre/fs;
wstop=2*PI*lowFre/fs;
Wpass=cos(wpass/2)/sin(wpass/2);
Wstop=cos(wstop/2)/sin(wstop/2);
Epass= sqrt(pow(10,Ap/10)-1);
Estop= sqrt(pow(10,As/10)-1);
e=Estop/Epass;
w=Wstop/Wpass;
order=int(log(e+sqrt(pow(e,2)-1))/log(w+sqrt(pow(w,2)-1)))+1;
if(order > 16)
{
throw std::logic_error("设计的滤波器不稳定！请重试！");
}
a=log(Estop+sqrt(pow(Estop,2)+1))/order;
OMEGA=Wstop/sinh(a);
}
void ChebyTwo::Den_Cheby2_Highpass()
{
int K=order/2;
double X[20];
double H[3];
double Angle;
double Wi;
double Wstop=cos(PI*lowFre/fs)/sin(PI*lowFre/fs);
int X_Length,H_Length;
if (order%2==0)
{
X[0]=1;
X_Length=1;
}
else
{
X[0]=1;
X[1]=-(OMEGA-1)/(OMEGA+1);
X_Length=2;
pDen[0] = X[0];
pDen[1] = X[1];
}
for (int i=1;i<=K;i++)
{
Angle=PI*(order-1+2*i)/(2*order);
Wi=Wstop/sin(Angle);
H[0]=1;
H[1]=-2*(1-pow(OMEGA,-2)-pow(Wi,-2))/(1-2*cos(Angle)/OMEGA+pow(OMEGA,-2)+pow(Wi,-2));
H[2]=(1+2*cos(Angle)/OMEGA+pow(OMEGA,-2)+pow(Wi,-2))/
(1-2*cos(Angle)/OMEGA+pow(OMEGA,-2)+pow(Wi,-2));
H_Length=3;
X_Length=Conv(H_Length,H,X_Length,X,pDen);
for (int j=0;j<X_Length;j++)
X[j]=pDen[j];
}
}
void ChebyTwo::Num_Cheby2_Highpass()
{
int K=order/2;
double X[20];
double H[3];
double Angle;
double Wi;
double Wstop=cos(PI*lowFre/fs)/sin(PI*lowFre/fs);
int X_Length,H_Length;
if (order%2==0)
{
X[0]=1;
X_Length=1;
}
else
{
X[0]=OMEGA/(OMEGA+1);
X[1]=-X[0];
X_Length=2;
pNum[0] = X[0];
pNum[1] = X[1];
}
for (int i=1;i<=K;i++)
{
Angle=PI*(order-1+2*i)/(2*order);
Wi=Wstop/sin(Angle);
H[0]=(1+pow(Wi,-2))/(1-2*cos(Angle)/OMEGA+pow(OMEGA,-2)+pow(Wi,-2));
H[1]=-2*H[0]*(1-pow(Wi,-2))/(1+pow(Wi,-2));
H[2]=H[0];
H_Length=3;
X_Length=Conv(H_Length,H,X_Length,X,pNum);
for (int j=0;j<X_Length;j++)
X[j]=pNum[j];
}
}
Notch::Notch(PassType passType, double fs, double fre)
:IIRFilterDesign(passType, fs)
{
this->fre = fre;
GenerateCoe();
}
Notch::~Notch()
{
ReleaseArray();
}
void Notch::GenerateCoe()
{
order = 2;
AllocateArray();
GenerateNotch();
}
/*!----------------------------------------------------------------------------
* \brief
*
* 陷波器的中心频率是：
* \f[
* \omega_{0} = \frac{2\pi f_{0}}{f_{s}}
* \f]
*
* 陷波器只有两阶：
* \f[
* b_{1} = -2r\cos\omega_{0} \quad,\quad b_{2} = r^{2}
* \f]
* \f[
* a_{1} = -2R\cos\omega_{0} \quad,\quad a_{2} = R^{2}
* \f]
*
* ----------------------------------------------------------------------------*/
void Notch::GenerateNotch()
{
int f=fre;
double R=SETR; //越接近1，陷波器宽度越窄
double r=SETr; //越接近1，陷波器衰减越厉害
OMEGA=2*PI*f/fs;
pDen[1]=-2*R*cos(OMEGA);
pDen[2]=pow(R,2);
pNum[0]=1;
pNum[1]=-2*r*cos(OMEGA);
pNum[2]=pow(r,2);
}
/*!----------------------------------------------------------------------------
* \brief
* 幅频函数
*
* Notch 的幅频响应函数
* \f[
* \frac{(1 - 2r\cos(\omega - \Omega_{0}) + r^2)
* (1 - 2r\cos(\omega + \Omega_{0}) + r^2)}
* {(1 - 2R\cos(\omega - \Omega_{0}) + R^2)
* (1 - 2R\cos(\omega + \Omega_{0}) + R^2)}
* \f]
*
* \param freq 要求的频率点
*
* \return 幅频
*
* ----------------------------------------------------------------------------*/
double Notch::MagResponse(double freq)
{
double range = 0.0;
double r = SETr;
double R = SETR;
double freqw = (2*PI*freq)/fs;
range = (1-2*r*cos(freqw-OMEGA)+pow(r,2))*(1-2*r*cos(freqw+OMEGA)+pow(r,2))
/((1-2*R*cos(freqw-OMEGA)+pow(R,2))*(1-2*R*cos(freqw+OMEGA)+pow(R,2)));
return range;
}
#endif
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -