filterdesign.~cpp

脑电信号分析软件

    double Angle;
    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);

        H[0]=pow(OMEGA,2)/(1-2*OMEGA*cos(Angle)+pow(OMEGA,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];
    }
/*
	X[0] = 1;
	X_Length = 1;
	
	if(order % 2 == 0)
	{
		for (int i=1;i<=K;i++)
		{
			Angle=PI*(order-1+2*i)/(2*order);

			H[0]=pow(OMEGA,2)/(1-2*OMEGA*cos(Angle)+pow(OMEGA,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];
		}

	}
	else
	{
        H[0]=OMEGA/(OMEGA+1);
        H[1]=H[0];
        H_Length=2;

		for(int i = 1; i <= K; i++)
		{
			X_Length=Conv(H_Length,H,X_Length,X,pNum);

			for (int j=0;j<X_Length;j++)
				X[j]=pNum[j];
		}

	}*/

}

/*!----------------------------------------------------------------------------
 * \brief 生成 butterworth  heightpass 的阶数
 *
 * 做预畸和lowpass 不同 其它参考 Butter::Butter_Lowpass_N：
 * \f[
 *    \omega_{pass} = \frac{1.0}{\tan(\frac{\omega_{pass}}{2})} \qquad,\qquad
 *    \omega_{stop} = \frac{1.0}{\tan(\frac{\omega_{stop}}{2})}
 * \f]
 *
 * ----------------------------------------------------------------------------*/ 
void Butter::Butter_Highpass_N()
{
    double wpass;
    double wstop;
    double Wpass;
    double Wstop;
    double e,w;

    wpass=2*PI*heightFre/fs;
    wstop=2*PI*lowFre/fs;

    Wpass=1.0/tan(wpass/2);
    Wstop=1.0/tan(wstop/2);

    e=sqrt(pow(10,As/10)-1)/sqrt(pow(10,Ap/10)-1);
    w=Wstop/Wpass;
    order=int(log(e)/log(w))+1;

    if(order > 16)
    {
        throw std::logic_error("设计的滤波器不稳定！请重试！");
    }

    OMEGA=Wpass/pow(pow(10,Ap/10)-1,1.0/(2*order));
}

/*!----------------------------------------------------------------------------
 * \brief 求高通的分母
 *
 * 和低通的区别只在于把方程中的z 换成-z。参考 Butter::Den_Butter_Lowpass
 *
 * ----------------------------------------------------------------------------*/ 
void Butter::Den_Butter_Highpass()
{
    int K=order/2;
    double X[20];
    double H[3];
    double Angle;
    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);

        H[0]=1;
        H[1]=-2*(pow(OMEGA,2)-1)/(1-2*OMEGA*cos(Angle)+pow(OMEGA,2));
        H[2]=(1+2*OMEGA*cos(Angle)+pow(OMEGA,2))/(1-2*OMEGA*cos(Angle)+pow(OMEGA,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 求高通的分子
 *
 * 和低通的区别只在于把方程中的z 换成-z。参考 Butter::Num_Butter_Lowpass
 *
 * ----------------------------------------------------------------------------*/ 
void Butter::Num_Butter_Highpass()
{
    int K=order/2;
    double X[20];
    double H[3];
    double Angle;
    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);

        H[0]=pow(OMEGA,2)/(1-2*OMEGA*cos(Angle)+pow(OMEGA,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 最小止带衰减
 *
 * ----------------------------------------------------------------------------*/
ChebyOne::ChebyOne(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 使用默认值
 *
 * ----------------------------------------------------------------------------*/
ChebyOne::ChebyOne(PassType passType, double fs, double lowFre, double heightFre)
			: IIRFilterDesign(passType, fs)
{
	Init(lowFre, heightFre, 1, 60);

	GenerateCoe();
}

ChebyOne::~ChebyOne()
{
	ReleaseArray();
}

void ChebyOne::Init(double lowFre, double heightFre, double aPass, double aStop)
{
	this->lowFre = lowFre;
	this->heightFre = heightFre;

	Ap = aPass;
	As = aStop;
}

void ChebyOne::GenerateCoe()
{
	switch(passType)
	{
		case LOW_PASS:
			GenerateLowCoe();
			break;

		case HEIGHT_PASS:
			GenerateHeightCoe();
			break;

		default:
			break;
	}
}

void ChebyOne::GenerateLowCoe()
{
	Cheby_Lowpass_N();
                       
	AllocateArray();
    Den_Cheby_Lowpass();
	Num_Cheby_Lowpass();
}
void ChebyOne::GenerateHeightCoe()
{
	Cheby_Highpass_N();

	AllocateArray();
	Den_Cheby_Highpass();
	Num_Cheby_Highpass();
}

/*!----------------------------------------------------------------------------
 * \brief 
 *     幅频函数
 *
 * Chebyshev I 的幅频响应函数
 * \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]
 *
 * \param freq 要求的频率点
 *
 * \return 幅频 
 *
 * ----------------------------------------------------------------------------*/ 
double ChebyOne::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, freqW/Wpass), 2);

	range = 1.0/sqrt((1 + pow(Epass, 2)*CN_2));
	//range = 1.0/(1 + pow(Epass, 2)*CN_2);
	
		
	return range;
}


/*!----------------------------------------------------------------------------
 * \brief 
 *
 * 阶数N由下面公式求出：
 * \f[ 
 *   N_{exact} = \frac{\ln(e + \sqrt{e^{2} - 1})}{\ln(w + \sqrt{w^{2} - 1})}
 * \f] 
 *   
 * 定义a：
 * \f[ 
 *   a = \frac{1}{N}\ln(\frac{1}{\varepsilon_{pass}} + 
 *       \sqrt{\frac{1}{\varepsilon_{pass}^{\phantom{pass}{2}}} + 1})
 * \f] 
 * \f$ \Omega_{0} \f$ 为：
 * \f[ 
 *   \Omega_{0} = \Omega_{pass}\sinh a
 * \f] 
 *  
 * ----------------------------------------------------------------------------*/
void ChebyOne::Cheby_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(1.0/Epass+sqrt(1.0/pow(Epass,2)+1))/order;
    OMEGA=Wpass*sinh(a);
}

/*!----------------------------------------------------------------------------
 * \brief  生成低通系数分母
 *
 * 二阶基本节：
 * \f[
 * H_{i}(z) = \frac{G_{i}(1 + z^{-1})^2}{1 + a_{i1}z^{-1} +
 *            a_{i2}z^{-2}} \qquad i = 1,2,\ldots,K.
 * \f]
 *
 * 当阶数为奇的时候：
 * \f[
 * H_{0}(z) = \left\{ \begin{array}{ll}
 *            \sqrt{\frac{1}{1 + \varepsilon_{pass}^{\phantom{pass}{2}}}}& \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} \qquad, \qquad
 *   a_{01} = \frac{\Omega_{0} - 1}{\Omega_{0} + 1}
 * \f]
 *
 * 定义\f$ \theta_{i} \f$
 * \f[
 *   \theta_{i} = \frac{\pi}{2N}(N - 1 + 2i)
 * \f]
 *
 * \f$ a_{i1} a_{i2} \f$为：
 * \f[
 *   a_{i1} = \frac{2(\Omega_{0}^{\phantom{0}{2}} +
 *            \Omega_{i}^{\phantom{i}{2}} - 1)}{1 - 2\Omega_{0}\cos{\theta_{i}}
 *            + \Omega_{0}^{\phantom{0}{2}} + \Omega_{i}^{\phantom{i}{2}}}
 * \f]
 * \f[
 *   a_{i2} = \frac{1 + 2\Omega_{0}\cos\theta_{i} +
 *            \Omega_{0}^{\phantom{0}{2}} + 
 *            \Omega_{i}^{\phantom{i}{2}}}{1 -
 *            2\Omega_{0}\cos{\theta_{i}} +
 *            \Omega_{0}^{\phantom{0}{2}} + \Omega_{i}^{\phantom{i}{2}}}
 * \f]
 * 
 * ----------------------------------------------------------------------------*/
void ChebyOne::Den_Cheby_Lowpass()
{
    int K=order/2;
    double X[20];
    double H[3];
    double Angle;
    double Wi;
    double Wpass=tan(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=Wpass*sin(Angle);

        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];
    }
}

/*!----------------------------------------------------------------------------
 * \brief  求分子的系数， 参照求分母ChebyOne::Den_Cheby_Lowpass()
 *
 * 定义\f$ \Omega_{i} \f$
 * \f[
 *   \Omega_{i} = \Omega_{pass}\sin\theta_{i} \quad, \quad i= 1,2,3\ldots,K
 * \f]
 *
 * \f[
 * G_{i} = \frac{\Omega_{0}^{\phantom{0}{2}} +
 *         \Omega_{i}^{\phantom{i}{2}}}{1 - 2\Omega_{0}\cos{\theta_{i}} +
 *         \Omega_{0}^{\phantom{0}{2}} + \Omega_{i}^{\phantom{i}{2}}}\\
 * \f]
 *
 * 
 * ----------------------------------------------------------------------------*/
void ChebyOne::Num_Cheby_Lowpass()
{
    int K=order/2;
    double X[20];
    double H[3];
    double Angle;
    double Wi;
    double Wpass=tan(PI*lowFre/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];
    }
}

void ChebyOne::Cheby_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(1.0/Epass+sqrt(1.0/pow(Epass,2)+1))/order;
    OMEGA=Wpass*sinh(a);
}

void ChebyOne::Den_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]=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=Wpass*sin(Angle);