几个c程序，希望有用
include <math.h>
/* In the VMS compiler, M_PI is not defined in math.h */
#ifdef vms
#define M_PI 3.14159265358979323846
#endif
#ifndef KR
#define M_PI 3.14159265358979323846
#endif
#define FALSE 0
#define TRUE 1
#define BIG 1e10
#define SMALL 1e-10
#define ORDER5 1e-5
#define ORDER4 1e-4
#define ABS(x) ((x) >= 0 ? (x) : -(x))
#define MIN(a,b) ((a) <= (b) ? (a) : (b))
#define MAX(a,b) ((a) >= (b) ? (a) : (b))
/* Table of constant values */
static long neg_i1 = -1;
static long pos_i1 = 1;
static long pos_i2 = 2;
static long pos_i4 = 4;
static long pos_i5 = 5;
static float neg_f1 = -1.0;
static float pos_f2 = 2.0;
static double pos_d10 = 10.0;
typedef struct {
float r, i;
} complex;
typedef struct {
double r, i;
} doublecomplex;
/* .................End of include file.................*/
/* SPBFCT 11/14/85 */
/* GENERATES (I1)!/(I1-I2)!=I1*(I1-1)*...*(I1-I2+1). */
/* NOTE: 0!=1 AND SPBFCT(I,I)=SPBFCT(I,I-1)=I!. */
#ifndef KR
double spbfct(long *i1, long *i2)
#else
double spbfct(i1, i2)
long *i1, *i2;
#endif
{
/* Local variables */
long i;
double ret_val;
ret_val = 0.0;
if (*i1 < 0 || *i2 < 0 || *i2 > *i1)
{
return(ret_val);
}
ret_val = 1.0;
for (i = *i1 ; i >= (*i1 - *i2 + 1) ; --i)
{
ret_val *= (double) i;
}
return(ret_val);
} /* spbfct */
/* SPBILN 11/13/85 */
/* CONVERTS ANALOG H(S) TO DIGITAL H(Z) VIA BILINEAR TRANSFORM */
/* ANALOG TRANSFER FUNCTION DIGITAL TRANSFER FUNCTION */
/* D(L)*S**L+.....+D(0) B(0)+......+B(L)*Z**-L */
/* H(S)=--------------------- H(Z)=---------------------- */
/* C(L)*S**L+.....+C(0) 1+.......+A(L)*Z**-L */
/* H(S) IS ASSUMED TO BE PRE-SCALED AND PRE-WARPED */
/* LN SPECIFIES THE LENGTH OF THE COEFFICIENT ARRAYS */
/* FILTER ORDER L IS COMPUTED INTERNALLY */
/* WORK IS AN INTERNAL ARRAY (2D) SIZED TO MATCH COEF ARRAYS */
/* IERROR=0 NO ERRORS DETECTED IN TRANSFORMATION */
/* 1 ALL ZERO TRANSFER FUNCTION */
/* 2 INVALID TRANSFER FUNCTION; Y(K) COEF=0 */
#ifndef KR
void spbiln(float *d, float *c, long *ln, float *b, float *a, float *work, long *error)
#else
void spbiln(d, c, ln, b, a, work, error)
long *ln, *error;
float *d, *c, *b, *a, *work;
#endif
{
/* Local variables */
long i, j, l, zero_func, work_dim1;
float scale, tmp;
double atmp;
zero_func = TRUE;
for (i = *ln ; i >= 0 && zero_func ; --i)
{
if (c[i] != 0.0 || d[i] != 0.0)
{
zero_func = FALSE;
}
}
if ( zero_func )
{
*error = 1;
return;
}
work_dim1 = *ln + 1;
l = i + 1;
for (j = 0 ; j <= l ; ++j)
{
work[j * work_dim1] = 1.0;
}
tmp = 1.0;
for (i = 1 ; i <= l ; ++i)
{
tmp = tmp * (float) (l - i + 1) / (float) i;
work[i] = tmp;
}
for (i = 1 ; i <= l ; ++i)
{
for (j = 1 ; j <= l ; ++j)
{
work[i + j * work_dim1] = work[i + (j - 1) * work_dim1]
- work[i - 1 + j * work_dim1]
- work[i - 1 + (j - 1) * work_dim1];
}
}
for (i = l ; i >= 0 ; --i)
{
b[i] = 0.0;
atmp = 0.0;
for (j = 0 ; j <= l ; ++j)
{
b[i] += work[i + j * work_dim1] * d[j];
atmp += (double) work[i + j * work_dim1] * c[j];
}
scale = (double) atmp;
if (i != 0)
{
a[i - 1] = (double) atmp;
}
}
if (scale == 0.0)
{
*error = 2;
return;
}
b[0] /= scale;
for (i = 1 ; i <= l ; ++i)
{
b[i] /= scale;
a[i - 1] /= scale;
}
if (l < *ln)
{
for (i = l + 1 ; i <= *ln ; ++i)
{
b[i] = 0.0;
a[i - 1] = 0.0;
}
}
*error = 0;
return;
} /* spbiln */
/* SPBSSL 11/13/85 */
/* GENERATES ANALOG FILTER COEFFICIENTS FOR LTH ORDER */
/* NORMALIZED LOWPASS BESSEL FILTER */
/* COEFFICIENTS ARE RETURNED IN ARRAYS D AND C */
/* LN SPECIFIES ARRAY SIZE (LN>=L) */
/* WSCL CONTROLS FREQUENCY SCALING SUCH THAT RESPONSE AT 1 RAD/SEC */
/* IS EQUAL TO THAT OF UNSCALED H(S) AT WSCL RAD/SEC */
/* IERROR=0 NO ERRORS DETECTED */
/* 1 INVALID FILTER ORDER (L<=0 OR L>=LN) */
/* 2 INVALID SCALE PARAMETER (WSCL<=0) */
#ifndef KR
void spbssl(long *l, float *wscl, long *ln, float *d, float *c, long *error)
#else
void spbssl(l, wscl, ln, d, c, error)
long *l, *ln, *error;
float *d, *c, *wscl;
#endif
{
/* Builtin functions */
long pow_ii();
double pow_ri();
/* Local variables */
double spbfct();
long k, ll, tmp_int1, tmp_int2;
if (*l <= 0 || *l > *ln)
{
*error = 1;
return;
}
if (*wscl <= 0.0)
{
*error = 2;
return;
}
for (ll = 0 ; ll <= *ln ; ++ll)
{
d[ll] = 0.0;
c[ll] = 0.0;
}
for (k = 0 ; k <= *l ; ++k)
{
tmp_int1 = (*l * 2) - k;
tmp_int2 = *l - k;
c[k] = pow_ri(wscl, &k)
* spbfct(&tmp_int1,&tmp_int1) / (pow_ii(&pos_i2, &tmp_int2)
* spbfct(&k,&k) * spbfct(&tmp_int2,&tmp_int2));
}
d[0] = c[0];
*error = 0;
return;
} /* spbssl */
/* SPBWCF 11/13/85 */
/* GENERATES KTH SECTION COEFFICIENTS FOR LTH ORDER NORMALIZED */
/* LOWPASS BUTTERWORTH FILTER */
/* SECOND ORDER SECTIONS: K<=(L+1)/2 */
/* ODD ORDER L: FINAL SECTION WILL CONTAIN 1ST ORDER POLE */
/* LN DEFINES COEFFICIENT ARRAY SIZE */
/* ANALOG COEFFICIENTS ARE RETURNED IN D AND C */
/* IERROR=0 NO ERRORS DETECTED */
/* 1 INVALID FILTER ORDER L */
/* 2 INVALID SECTION NUMBER K */
#ifndef KR
void spbwcf(long *l, long *k, long *ln, float *d, float *c, long *error)
#else
void spbwcf(l, k, ln, d, c, error)
long *l, *k, *ln, *error;
float *d, *c;
#endif
{
/* Local variables */
long i;
float tmp;
if (*l <= 0)
{
*error = 1;
return;
}
if (*k <= 0 || *k > ((*l + 1) / 2))
{
*error = 2;
return;
}
d[0] = 1.0;
c[0] = 1.0;
for (i = 1 ; i <= *ln ; ++i)
{
d[i] = 0.0;
c[i] = 0.0;
}
tmp = (float) *k - (*l + 1.0) / 2.0;
if (tmp == 0.0)
{
c[1] = 1.0;
}
else
{
c[1] = (float) (-2.0 * cos((double) ((*k * 2) + *l - 1)
* M_PI / (double) (*l * 2)));
c[2] = 1.0;
}
*error = 0;
return;
} /* spbwcf */
还有其它的几个例子
3145
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
