matlab cheb2ord函数 C实现

cheb2ord

void cheb2ord(double* w,int nw,double rp,double rs,int* order,double* wn)
{
    int type=0;
	double* wp=(double*)malloc(sizeof(double)*nw/2);
	double* temp=(double*)malloc(sizeof(double)*nw/2);
	double* ws=(double*)malloc(sizeof(double)*nw/2);
	double* WA=(double*)malloc(sizeof(double)*nw/2);
	double minWA;

	for(int i=0;i<nw/2;i++)
	{
		wp[i]=w[i];
		ws[i]=w[nw/2+i];
        temp[i]=ws[i];
	}

	if(nw==2)
	{
		type=0;
	}
	else if(nw==4)
	{
		type=2;
	}
	if(wp[0]<ws[0])
	{
		type+=1;
	}
	else
	{
		type+=2;
	}


	for(int i=0;i<nw/2;i++)
	{
		wp[i]=tan(PI*wp[i]/2.0);
		ws[i]=tan(PI*ws[i]/2.0);
	}

	if (type == 1)
	{
		for (int i = 0; i < nw / 2; i++)
		{
			WA[i] = ws[i] / wp[i];
		}
	}
	else if(type==2)
	{
		for (int i = 0; i < nw / 2; i++)
		{
			WA[i] = wp[i] / ws[i];
		}
	}
	else if (type == 3)
	{
		wp[0] = lclfminbnd(wp[0], ws[0] - 1e-12,0, wp, ws, rs, rp, 1);
		wp[1] = lclfminbnd(ws[1] + 1e-12, wp[1],1, wp, ws, rs, rp, 1);
		for (int i = 0; i < nw / 2; i++)
		{
			WA[i] = (ws[i] * (wp[0] - wp[1])) / (ws[i] * ws[i] - wp[0] * wp[1]);
		}
	}
	else if(type==4)
	{
		for (int i = 0; i < nw / 2; i++)
		{
			WA[i] = (ws[i] * ws[i] - wp[0] * wp[1]) / (ws[i] * (wp[0] - wp[1]));
		}
	}

	minWA=INFINITY;
	for (int i = 0; i < nw / 2; i++)
	{
		double absValue = fabs(WA[i]);
		minWA = fmin(minWA, absValue);
	}

    double argument = sqrt((pow(10, 0.1 * fabs(rs)) - 1) / (pow(10, 0.1 * fabs(rp)) - 1));
	*order =  ceil(acosh(argument) / acosh(minWA));

    for (int i = 0; i < nw / 2; i++)
    {
        wn[i]=temp[i];
    }
    

	free(wp);
	free(ws);
	free(WA);
}

其中lclfminbnd在matlab buttord函数 C实现中已经实现。

和matlab函数对比，验证无误