ellipord
void ellipord(double *w, int nw, double rp, double rs, int *order, double *wn)
{
int type=0;
double* wp=(double*)malloc(sizeof(double)*nw/2);
double* ws=(double*)malloc(sizeof(double)*nw/2);
double* temp=(double*)malloc(sizeof(double)*nw/2);
double* WA=(double*)malloc(sizeof(double)*nw/2);
double* WN=(double*)malloc(sizeof(double)*nw/2);
for(int i=0;i<nw/2;i++)
{
wp[i]=w[i];
temp[i]=wp[i];
ws[i]=w[nw/2+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)
{
// low
for (int i = 0; i < nw / 2; i++)
{
WA[i] = ws[i] / wp[i];
}
*order = findelliporder(WA, nw / 2, rp, rs);
}
else if (type == 2)
{
// high
for (int i = 0; i < nw / 2; i++)
{
WA[i] = wp[i] / ws[i];
}
*order = findelliporder(WA, nw / 2, rp, rs);
}
else if (type == 3)
{
double Fpass1 = wp[0];
double Fpass2 = wp[1];
double Fstop1 = ws[0];
double Fstop2 = ws[1];
double c = sin(PI * (Fpass1 + Fpass2)) / (sin(PI * Fpass1) + sin(PI * Fpass2));
double wpa = fabs(sin(PI * Fpass2) / (cos(PI * Fpass2) - c));
double ws1 = sin(PI * Fstop1) / (cos(PI * Fstop1) - c);
double ws2 = sin(PI * Fstop2) / (cos(PI * Fstop2) - c);
double wsa = fmin(fabs(ws1), fabs(ws2));
double* w=(double*)malloc(sizeof(double)*2);
w[0]=wpa;
w[1]=wsa;
*order=ellipords(w,2,rp,rs);
free(w);
}
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]));
}
*order = findelliporder(WA, nw / 2, rp, rs);
}
for (int i = 0; i < nw / 2; i++)
{
wn[i]=temp[i];
}
free(temp);
free(wp);
free(ws);
free(WN);
free(WA);
}
其中
int ellipords(double *w, int nw, double rp, double rs)
{
int order=0;
int type=0;
double* wp=(double*)malloc(sizeof(double)*nw/2);
double* ws=(double*)malloc(sizeof(double)*nw/2);
double* WA=(double*)malloc(sizeof(double)*nw/2);
double* WN=(double*)malloc(sizeof(double)*nw/2);
double W0;
for(int i=0;i<nw/2;i++)
{
wp[i]=w[i];
ws[i]=w[nw/2+i];
}
if(nw==2)
{
type=0;
}
else if(nw==4)
{
type=2;
}
if(wp[0]<ws[0])
{
type+=1;
}
else
{
type+=2;
}
if (type == 1)
{
// low
for (int i = 0; i < nw / 2; i++)
{
WA[i] = ws[i] / wp[i];
}
order = findelliporder(WA, nw / 2, rp, rs);
}
else if (type == 2)
{
// high
for (int i = 0; i < nw / 2; i++)
{
WA[i] = wp[i] / ws[i];
}
order = findelliporder(WA, nw / 2, rp, rs);
}
else if (type == 3)
{
for (int i = 0; i < nw / 2; i++)
{
wp[i] = tan(PI * wp[i] / 2.0);
ws[i] = tan(PI * ws[i] / 2.0);
}
double Fpass1 = wp[0];
double Fpass2 = wp[1];
double Fstop1 = ws[0];
double Fstop2 = ws[1];
double c = sin(PI * (Fpass1 + Fpass2)) / (sin(PI * Fpass1) + sin(PI * Fpass2));
double wpa = fabs(sin(PI * Fpass2) / (cos(PI * Fpass2) - c));
double ws1 = sin(PI * Fstop1) / (cos(PI * Fstop1) - c);
double ws2 = sin(PI * Fstop2) / (cos(PI * Fstop2) - c);
double wsa = fmin(fabs(ws1), fabs(ws2));
double* w=(double*)malloc(sizeof(double)*2);
w[0]=wpa;
w[1]=wsa;
order=ellipords(w,2,rp,rs);
free(w);
}
else if(type ==4)
{
for (int i = 0; i < nw / 2; i++)
{
WA[i]=(ws[i]*ws[i] - wp[0]*wp[1])/(ws[0]*(wp[0]-wp[1]));
}
order = findelliporder(WA, nw / 2, rp, rs);
}
free(wp);
free(ws);
free(WN);
free(WA);
return order;
}
其中
int findelliporder(double* wa,int size,double rp,double rs)
{
int order;
double minWA = INFINITY;
double epsilon;
double k1;
double k;
// Assuming ellipke function is already defined
double *ks = (double *)malloc(sizeof(double) * 2);
double *ks1 = (double *)malloc(sizeof(double) * 2);
for (int i = 0; i < size; i++)
{
double absValue = fabs(wa[i]);
minWA = fmin(minWA, absValue);
}
epsilon = sqrt(pow(10, 0.1 * rp) - 1.0);
k1 = epsilon / sqrt(pow(10, 0.1 * rs) - 1.0);
k = 1.0 / minWA;
ks[0] = k * k;
ks[1] = 1 - k * k;
ks1[0] = k1 * k1;
ks1[1] = 1 - k1 * k1;
double *capk = ellipke(ks, 2);
double *capk1 = ellipke(ks1, 2);
order = ceil((capk[0] * capk1[1]) / (capk[1] * capk1[0]));
free(ks);
free(ks1);
free(capk);
free(capk1);
return order;
}
又出现函数
double* ellipke(double* m, int size)
{
double tol=EPS;
double* e=(double*)malloc(sizeof(double)*size);
double* k=(double*)malloc(sizeof(double)*size);
for (int i = 0; i < size; i++) {
double mi = m[i];
double classin = (mi == 0.0) ? 0.0 : (mi == 1.0) ? 1.0 : 2.0 * mi;
if (mi == 0) {
k[i] = 0;
e[i] = 0;
continue;
}
if (mi == 1) {
k[i] = INFINITY;
e[i] = INFINITY;
continue;
}
double a0 = 1;
double b0 = sqrt(1 - mi);
double c0 = NAN;
double s0 = mi;
int i1 = 0;
double mm = INFINITY;
while (mm > tol) {
double a1 = (a0 + b0) / 2;
double b1 = sqrt(a0 * b0);
double c1 = (a0 - b0) / 2;
i1++;
double w1 = pow(2, i1) * pow(c1, 2);
mm = fmax(w1, 0);
s0 += w1;
a0 = a1;
b0 = b1;
c0 = c1;
}
k[i] = PI / (2 * a0);
e[i] = k[i] * (1 - s0 / 2);
if (mi == 1) {
k[i] = INFINITY;
e[i] = 1; // or some other value based on your specific case
}
}
return k;
}