📄 remez.c
字号:
nzz = nfcns + 2;
niter = 0;
comp = 0;
luck = 0;
y1 = 0;
nut1 = 0;
loop1 = 1;
printf("Iteration");
while(niter++ < ITRMAX){ /* LOOP 1 */
if(loop1 == 0) break; /* OUT OF LOOP1 */
printf(" %d",niter);
if(!(niter%20)) putchar('\n'); /* newline every 20 iterations */
iext[nzz] = ngrid + 1;
for(j=1;j<=nz;j++) x[j] = cos(PI2 * grid[iext[j]]);
jet = (int)((nfcns - 1)/15 +1);
for(j=1;j<=nz;j++) ad[j] = d(j,nz,jet);
dnum = 0.0;
dden = 0.0;
k=1;
for(j=1;j<=nz;j++){
l=iext[j];
dnum += ad[j] * des[l];
if(wt[l] == 0.0) wt[l] = EPS;
dden += k*ad[j]/wt[l];
k = -k;
}
dev = dnum/dden;
nu = 1;
if(dev > 0.0) nu = -1;
dev *= -nu;
k = nu;
for(j=1;j<=nz;j++){
l = iext[j];
if(wt[l] == 0.0) wt[l] = EPS;
y[j] = des[l] + k*dev/wt[l];
k = -k;
}
if(dev < devl){ /* "ouch" */
printf(" ********* FAILURE TO CONVERGE **********\n");
printf("Probable cause is machine rounding error\n");
printf("The impulse response may be correct\n");
printf("Check with frequency response\n");
break; /* OUT OF LOOP 1 */
}
devl = dev;
jchnge = 0;
k1 = iext[1];
knz = iext[nz];
klow = 0;
nut = -nu;
j = 1;
/* SEARCH FOR THE EXTREMAL FREQUENCIES OF THE BEST APPROXIMATION */
while(1){ /* LOOP 2 */
if(j == nzz) ynz = comp;
if(j >= nzz){
if(j > nzz){
if(luck > 9){
kn = iext[nzz];
for(j=1;j<=nfcns;j++) iext[j] = iext[j+1];
iext[nz] = kn;
break; /* OUT OF LOOP 2 */
}
else{
if(comp > y1) y1 = comp;
k1 = iext[nzz];
lable325: l = ngrid + 1; /* label325 ? */
klow = knz;
nut = -nut1;
comp = y1 * (1.00001);
out1 = 0;
out2 = 0;
while(1){ /* LOOP 3 */
l--;
if(l <= klow){
out1 = 1;
break; /* out of LOOP 3 */
}
else{
err = (gee(l,nz) - des[l]) * wt[l];
dtemp = nut * err - comp;
if(dtemp > 0.0){ /* BREAK OUT OF LOOP 3 */
out2 = 1;
break; /* out of LOOP 3 */
}
}
} /* END OF while LOOP 3 */
if(out1){
if(luck == 6){
if(jchnge > 0) break; /* OUT OF LOOP 2 */
loop1 = 0;
break; /* OUT OF LOOP 2 and LOOP 1 */
}
for(j=1;j<=nfcns;j++) iext[nzz - j] = iext[nz - j];
iext[1] = k1;
break; /* OUT OF LOOP 2 */
} /* END OF if(out1) */
if(out2){
j = nzz;
comp = nut * err;
luck += 10;
goto lable235;
}
} /* END OF if(luck > 9) else */
} /* END OF if(j > nzz) */
else{
if(k1 > iext[1]) k1 = iext[1];
if(knz < iext[nz]) knz = iext[nz];
nut1 = nut;
nut = -nu;
l = 0;
kup = k1;
comp = ynz * (1.00001);
luck = 1;
out1 = 0;
out2 = 0;
while(1){ /* LOOP 4 */
l++;
if(l >= kup){
out1 = 1;
break; /* out of LOOP 4 */
}
err = (gee(l,nz) - des[l]) * wt[l];
dtemp = nut * err - comp;
if(dtemp > 0.0){ /* BREAK OUT OF LOOP 4 */
out2 = 1;
break; /* out of LOOP 4 */
}
} /* END OF while LOOP 4 */
if(out1){
luck = 6;
goto lable325;
}
if(out2){
j = nzz;
comp = nut * err;
goto lable210;
}
} /* END OF if(j > nzz) else */
} /* END OF if(j >= nzz) */
else{
kup = iext[j+1];
l = iext[j] + 1;
nut = - nut;
if(j == 2) y1 = comp;
comp = dev;
if(l >= kup){
lable220: l--;
out1 = 0;
out2 = 0;
out3 = 0;
while(1){ /* LOOP 5 */
l--;
if(l <= klow){
out1 = 1;
break; /* OUT OF LOOP 5 */
}
err = (gee(l,nz) - des[l]) * wt[l];
dtemp = nut * err - comp;
if(dtemp > 0.0){
out2 = 1;
break; /* OUT OF LOOP 5 */
}
if(jchnge > 0){
out3 = 1;
break; /* OUT OF LOOP5 */
}
} /* END OF while LOOP 5 */
if(out1){ /* if we exited LOOP 5 this way */
l = iext[j] + 1;
if(jchnge > 0){
iext[j] = l - 1;
j++;
klow = l - 1;
jchnge++;
continue; /* to top of LOOP 2 */
}
out1 = 0;
out2 = 0;
while(1){ /* LOOP 6 */
l++;
if(l >= kup){
out1 = 1;
break; /* OUT OF LOOP 6 */
}
err = (gee(l,nz) - des[l]) * wt[l];
dtemp = nut * err - comp;
if(dtemp > 0.0){ /* BREAK OUT OF LOOP 6 */
comp = nut * err;
out2 = 1;
break; /* OUT OF LOOP 6 */
}
} /* END OF while LOOP 6 */
if(out1){
klow = iext[j];
j++;
continue; /* to top of LOOP 2 */
}
if(out2) goto lable210;
} /* END OF if(out1) */
else if(out2){
comp = nut * err;
lable235: while(1){ /* LOOP 7 */
l--;
if(l <= klow) break; /* OUT OF LOOP 7 */
err = (gee(l,nz) - des[l]) * wt[l];
dtemp = nut * err - comp;
if(dtemp <= 0.0) break; /* OUT OF LOOP 7 */
comp = nut * err;
} /* END OF while LOOP 7 */
klow = iext[j];
iext[j] = l+1;
j++;
jchnge++;
continue; /* to top of LOOP 2 */
} /* END OF else if(out2) */
else if(out3){
klow = iext[j];
j++;
continue; /* to top of LOOP 2 */
}
} /* END OF if(l >= kup) */
else{
err = (gee(l,nz) - des[l]) * wt[l];
dtemp = nut * err - comp;
if(dtemp <= 0.0) goto lable220;
comp = nut * err;
lable210: while(1){ /* LOOP 8 */
l++;
if(l >= kup) break; /* OUT OF LOOP 8 */
err = (gee(l,nz) - des[l]) * wt[l];
dtemp = nut * err - comp;
if(dtemp <= 0.0) break; /* OUT OF LOOP 8 */
comp = nut * err;
}
iext[j] = l - 1;
j++;
klow = l - 1;
jchnge++;
continue; /* to top of LOOP 2 */
}
} /* END OF if(j >= nzz) else */
} /* END OF while LOOP 2 */
} /* END OF while LOOP 1 */
/* CALCULATION OF THE COEFFICIENTS OF THE BEST APPROXIMATION
/* USING THE INVERSE DISCRETE FOURIER TRANSFORM */
lable400:
nm1 = nfcns - 1;
fsh = 1.0e-6;
gtemp = grid[1];
x[nzz] = -2.0;
cn = 2 * nfcns - 1;
delf = 1.0 / cn;
l = 1;
kkk = 0;
if(edge[1] < 0.01 && edge[2*nbands] > 0.49) kkk = 1;
if(nfcns <= 3) kkk = 1;
if(kkk != 1) {
dtemp = cos(PI2 * grid[1]);
dnum = cos(PI2 * grid[ngrid]);
tmp = dtemp - dnum;
if(tmp == 0.0) tmp = EPS;
aa = 2.0/tmp;
bb = -(dtemp + dnum)/tmp;
}
for(j=1;j<=nfcns;j++){
ft = (j - 1) * delf;
xt = cos(PI2 * ft);
if(kkk != 1){
xt = (xt - bb)/aa;
ft = acos(xt)/PI2;
}
out1 = 0;
out2 = 0;
while(1){ /* LOOP 9 */
xe = x[l];
if(xt > xe){
out1 = 1;
break; /* OUT OF LOOP 9 */
}
if((xe - xt) < fsh){
out2 = 1;
break; /* OUT OF LOOP 9 */
}
l++;
} /* END OF while LOOP 9 */
if(out1){
if((xt - xe) < fsh) a[j] = y[l];
else{
grid[1] = ft;
a[j] = gee(1,nz);
}
}
if(out2) a[j] = y[l];
if(l > 1) l--;
} /* END for() LOOP */
grid[1] = gtemp;
dden = PI2/cn;
for(j=1;j<=nfcns;j++){
dtemp = 0.0;
dnum = (j-1) * dden;
if(nm1 >= 1) for(k=1;k<=nm1;k++)
dtemp += a[k+1] * cos(dnum * k);
alpha[j] = 2 * dtemp + a[1];
} /* END for() LOOP */
for(j=2;j<=nfcns;j++) alpha[j] *= 2/cn;
alpha[1] /= cn;
if(kkk != 1){
p[1] = 2.0 * alpha[nfcns] * bb + alpha[nm1];
p[2] = 2.0 * alpha[nfcns] * aa;
q[1] = alpha[nfcns - 2] - alpha[nfcns];
for(j=2;j<=nm1;j++){
if(j >= nm1){
aa *= 0.5;
bb *= 0.5;
}
p[j+1] = 0.0;
for(k=1;k<=j;k++){
a[k] = p[k];
p[k] = 2.0 * bb * a[k];
}
p[2] += 2.0 * aa * a[1];
jm1 = j-1;
for(k=1;k<=jm1;k++) p[k] += q[k] + aa * a[k+1];
jp1 = j + 1;
for(k=3;k<=jp1;k++) p[k] += aa * a[k-1];
if(j != nm1){
for(k=1;k<=j;k++) q[k] = - a[k];
q[1] += alpha[nfcns - 1 - j];
}
} /* END OF for() LOOP */
for(j=1;j<=nfcns;j++) alpha[j] = p[j];
} /* END OF if(kkk != 1) */
if(nfcns <= 3){
alpha[nfcns + 1] = 0.0;
alpha[nfcns + 2] = 0.0;
}
} /* END REMEZ EXCHANGE ALGORITHM */
/* FUNCTION TO CALCULATE THE LAGRANGE INTERPOLATION
/* COEFFICIENTS FOR USE IN THE FUNCTION GEE */
double d(k,n,m) int k,n,m;{
int j,el;
double de,q;
de = 1.0;
q = x[k];
for(el=1;el<=m;el++)
for(j=el;j<=n;j+=m)
if(j != k) de *= 2.0 * (q - x[j]);
if(de == 0.0) de = EPS;
return(1.0/de);
}
/* FUNCTION TO EVALUATE THE FREQUENCY RESPONSE USING THE
/* LAGRANGE INTERPOLATION FORMULA IN THE BARYCENTRIC FORM */
double gee(k,n) int k,n;{
int j;
double c,de,p,xf;
de = 0.0;
p = 0.0;
xf = cos(PI2*grid[k]);
for(j=1;j<=n;j++){
c = xf - x[j];
if(c == 0.0) c = EPS;
c = ad[j]/c;
de += c;
p += c * y[j];
}
if(de == 0.0) de = EPS;
return(p/de);
}
void rerror(error_text) char error_text[];{
printf("%s\n",error_text);
printf("...now exiting to system...\n");
exit(1);
}
/* FLOAT INPUT */
double get_float(char *title_string,double low_limit,double up_limit)
{
double i;
int error_flag;
char *endcp;
char ctemp[128];
if(low_limit > up_limit){
printf("\nLimit error lower > upper\n");
exit(1);}
do {
printf("\n%s [%G to %G] ? ", title_string, low_limit, up_limit);
gets(ctemp);
i = strtod(ctemp, &endcp);
error_flag = (ctemp == endcp) || (*endcp != '\0');
} while (i < low_limit || i > up_limit || error_flag);
return(i);
}
/* int INPUT */
int get_int(char *title_string,int low_limit,int up_limit)
{
int i;
int error_flag;
char *endcp;
char ctemp[128];
if(low_limit > up_limit){
printf("\nLimit error lower > upper\n");
exit(1);}
do {
printf("\n%s [%d to %d] ? ", title_string, low_limit, up_limit);
gets(ctemp);
i = strtol(ctemp,&endcp,10);
error_flag = (ctemp == endcp) || (*endcp != '\0');
} while (i < low_limit || i > up_limit || error_flag);
return(i);
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -