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函数对比,验证无误