99 lines
2.0 KiB
C
99 lines
2.0 KiB
C
// C port of http://ldesoras.free.fr/prod.html#src_hiir
|
|
#include <stdlib.h>
|
|
#include <stdio.h>
|
|
#include <math.h>
|
|
#include <errno.h>
|
|
|
|
const double PI = 3.1415926535897932384626433832795;
|
|
|
|
void
|
|
compute_transition_param(double *kp, double *qp, double transition)
|
|
{
|
|
double k = *kp;
|
|
double q = *qp;
|
|
k = tan((1 - transition*2)*PI/4);
|
|
k *= k;
|
|
double kksqrt = pow(1 - k*k, 0.25);
|
|
double e = 0.5*(1 - kksqrt)/(1 + kksqrt);
|
|
double e4 = e*e*e*e;
|
|
q = e*(1 + e4*(2 + e4*(15 + 150*e4)));
|
|
*kp = k;
|
|
*qp = q;
|
|
}
|
|
|
|
double
|
|
compute_acc_numden(double q, int order, int c, int den)
|
|
{
|
|
den = den == 1;
|
|
int i = den;
|
|
int sign = den ? -1 : 1;
|
|
double sum = 0;
|
|
double q_ii;
|
|
do {
|
|
int i2 = i + 1 - den;
|
|
q_ii = pow(q, i*i2);
|
|
q_ii *= sin((i + i2)*c*PI/order + den*PI/2);
|
|
q_ii *= sign;
|
|
sum += q_ii;
|
|
|
|
sign = -sign;
|
|
++i;
|
|
} while (fabs(q_ii) > 1e-100);
|
|
return sum;
|
|
}
|
|
|
|
double
|
|
compute_coef(int c, double k, double q, int order)
|
|
{
|
|
double num = compute_acc_numden(q, order, c, 0)*pow(q, 0.25);
|
|
double den = compute_acc_numden(q, order, c, 1) + 0.5;
|
|
double ww = num/den;
|
|
ww *= ww;
|
|
|
|
double x = sqrt((1 - ww*k)*(1 - ww/k))/(1 + ww);
|
|
double coef = (1 - x)/(1 + x);
|
|
|
|
return coef;
|
|
}
|
|
|
|
void
|
|
compute_coefs_spec_order_tbw(double coef_arr[], int n, double transition)
|
|
{
|
|
double k;
|
|
double q;
|
|
compute_transition_param(&k, &q, transition);
|
|
const int order = n*2 + 1;
|
|
|
|
/*printf("k: %.18f\nq: %.18f\n", k, q);*/
|
|
|
|
for (int i = 0; i < n; i++)
|
|
coef_arr[i] = compute_coef(i + 1, k, q, order);
|
|
}
|
|
|
|
int
|
|
main(int argc, char **argv)
|
|
{
|
|
if (argc != 3) {
|
|
fputs("usage: halfband COUNT TRANSITION\n", stderr);
|
|
return 1;
|
|
}
|
|
|
|
double tv[2];
|
|
for (int i = 1; i <= 2; i++) {
|
|
tv[i - 1]=strtold(argv[i], NULL);
|
|
if (errno) {
|
|
fprintf(stderr, "arg #%i failed to convert to double\n", i);
|
|
return 1;
|
|
}
|
|
}
|
|
int count = (int) tv[0];
|
|
|
|
double *coefs = malloc(count*sizeof(double));
|
|
compute_coefs_spec_order_tbw(coefs, count, tv[1]);
|
|
for (int i = 0; i < count; i++)
|
|
printf("%.18f\n", coefs[i]);
|
|
free(coefs);
|
|
|
|
return 0;
|
|
}
|