This is window.c in view mode; [Download] [Up]
#include <math.h>
#include <stdio.h>
/*
#include <malloc.h>
*/
/*-----------------------------------------------------------------------
WINDOW Subroutines
These are essentially C versions of the classical window subroutines in
the IEEE Programs for Digital Signal Processing. In order that
the code be more easily portable to PC's, complicated arithmetic
expressions have been broken into parts, dynamic array allocation
has replaced static declaration, and error checking has been
improved.
-----------------------------------------------------------------------*/
/*-----------------------------------------------------------------------
subroutine: triang
triangular window
-----------------------------------------------------------------------*/
triang(nf,w,n,ieo)
int nf, n, ieo;
float *w;
{
/*
nf = filter length in samples
w = window coefficients for half the window
n = half window length=(nf+1)/2
ieo = even - odd indication--ieo=0 for nf even
*/
float fn, xi;
int i;
fn = n;
for (i = 0; i < n; i++){
xi = i;
if (ieo == 0)
xi += 0.5;
w[i] = 1. - xi / fn;
}
return;
}
/*-----------------------------------------------------------------------
subroutine: hammin
generalized hamming window routine
window is w(n) = alpha + beta * cos( twopi*(n-1) / (nf-1) )
-----------------------------------------------------------------------*/
hammin(nf,w,n,ieo,alpha,beta)
int nf, n, ieo;
float *w, alpha, beta;
{
/*
nf = filter length in samples
w = window array of size n
n = half length of filter=(nf+1)/2
ieo = even odd indicator--ieo=0 if nf even
alpha = constant of window
beta = constant of window--generally beta=1-alpha
*/
float pi2, fn, fi;
int i;
pi2 = 8.0 * atan(1.0);
fn = nf - 1;
for (i = 0; i < n; i++) {
fi = i;
if (ieo == 0)
fi += 0.5;
w[i] = cos((pi2 * fi) / fn);
w[i] = alpha + beta * w[i];
}
return;
}
/*-----------------------------------------------------------------------
subroutine: kaiser
kaiser window
-----------------------------------------------------------------------*/
kaiser(nf,w,n,ieo,beta)
int nf, n, ieo;
float *w, beta;
{
/*
nf = filter length in samples
w = window array of size n
n = filter half length=(nf+1)/2
ieo = even odd indicator--ieo=0 if nf even
beta = parameter of kaiser window
*/
extern float ino();
float bes, xind, xi;
int i;
bes = ino(beta);
xind = (float)(nf-1)*(nf-1);
for (i = 0; i < n; i++) {
xi = i;
if (ieo == 0)
xi += 0.5;
xi = 4. * xi * xi;
xi = sqrt(1. - xi / xind);
w[i] = ino(beta * xi);
w[i] /= bes;
}
return;
}
/*-----------------------------------------------------------------------
function: ino
bessel function for kaiser window
-----------------------------------------------------------------------*/
float ino(x)
float x;
{
float y, t, e, de, sde, xi;
int i;
y = x / 2.;
t = 1.e-08;
e = 1.;
de = 1.;
for (i = 1; i <= 25; i++) {
xi = i;
de = de * y / xi;
sde = de * de;
e += sde;
if (e * t > sde)
break;
}
return(e);
}
/*-----------------------------------------------------------------------
subroutine: chebc
subroutine to generate chebyshev window parameters when
one of the three parameters nf,dp and df is unspecified
note: parameters must be passed as addresses so they can be changed
-----------------------------------------------------------------------*/
chebc(nf,dp,df,n,x0,xn)
int *nf, *n;
float *dp, *df, *x0, *xn;
{
/*
nf = filter length (in samples)
dp = filter ripple (absolute scale)
df = normalized transition width of filter
n = (nf+1)/2 = filter half length
x0 = (3-c0)/(1+c0) with c0=cos(pi*df) = chebyshev window constant
xn = nf-1
*/
extern float arccos(), coshin();
float pi, c0, c1, c2, x;
pi = 4. * atan(1.0);
if (*nf == 0) { /* dp,df specified, determine nf */
c1 = coshin((1. + *dp) / *dp);
c0 = cos(pi * *df);
x = coshin(1. / c0);
x = 1. + c1 / x;
*nf = x + 1.0;
*n = (*nf + 1) / 2;
*xn = *nf - 1.;
}
else if (*df != 0.0) { /* nf,df specified, determine dp */
*xn = *nf - 1.;
c0 = cos(pi * *df);
c1 = *xn * coshin(1. / c0);
*dp = cosh(c1) - 1.;
*dp = 1. / *dp;
}
else { /* nf,dp specified, determine df */
*xn = *nf - 1.;
c1 = coshin((1. + *dp) / *dp);
c2 = cosh(c1 / *xn);
*df = arccos(1. / c2) / pi;
}
*x0 = 3. - cos(2. * pi * *df);
*x0 /= (1. + cos(2. * pi * *df));
return;
}
/*-----------------------------------------------------------------------
subroutine: cheby
dolph chebyshev window design
-----------------------------------------------------------------------*/
cheby(nf,w,n,ieo,dp,df,x0,xn)
int nf, n, ieo;
float *w, dp, df, x0, xn;
{
/*
nf = filter length in samples
w = window array of size n
n = half length of filter = (nf+1)/2
ieo = even-odd indicator--ieo=0 for nf even
dp = window ripple on an absolute scale
df = normalized transition width of window
x0 = window parameter related to transition width
xn = nf-1
*/
float *pr, *pi;
float pie, fnf, alpha, beta, twopi, c1, c2, f, x, p, twn, sum,
xi, xj, t1, t2;
int i, j;
if ((pr = (float *) calloc(1024,sizeof(float))) == NULL){
fprintf(stderr,"cheby: can't calloc memory\n");
return;
}
if ((pi = (float *) calloc(1024,sizeof(float))) == NULL){
fprintf(stderr,"cheby: can't calloc memory\n");
return;
}
pie = 4. * atan(1.0);
xn = nf - 1;
fnf = nf;
alpha = (x0 + 1.) / 2.;
beta = (x0 - 1.) / 2.;
twopi = 2. * pie;
c2 = xn / 2.;
for (i = 0; i < nf; i++) {
xi = i;
f = xi / fnf;
x = cos(twopi * f);
x = alpha * x + beta;
if (x < -1.) /* in case of arithmetic round-off error */
x = -1.;
if (x > 1.){
p = coshin(x);
p = dp * cosh(c2 * p);
}
else {
p = arccos(x);
p = dp * cos(c2 * p);
}
pi[i] = 0.;
pr[i] = p;
/* for even length filters use a one-half sample delay,
also the frequency response is antisymmetric in frequency */
if (ieo != 1) {
pr[i] = p * cos(pie * f);
pi[i] = -p * sin(pie * f);
if (i > nf / 2 + 1){
pr[i] = -pr[i];
pi[i] = -pi[i];
}
}
}
/* use dft to give window */
twn = twopi / fnf;
for (i = 0; i < n; i++) {
xi = i;
sum = 0.;
for (j = 0; j < nf; j++) {
xj = j;
t1 = cos(twn * xj * xi);
t2 = sin(twn * xj * xi);
sum += pr[j] * t1 + pi[j] * t2;
}
w[i] = sum;
}
c1 = w[1];
for (i = 0; i < n; i++)
w[i] /= c1;
free(pi);
free(pr);
return;
}
/*-----------------------------------------------------------------------
function: coshin
function for hyperbolic inverse cosine of x
-----------------------------------------------------------------------*/
float coshin(x)
float x;
{
float c;
c = sqrt(x * x - 1.);
c = log(x + c);
return(c);
}
/*-----------------------------------------------------------------------
function: arccos
function for inverse cosine of x
-----------------------------------------------------------------------*/
float arccos(x)
float x;
{
float c, a;
if (x < 0.) {
a = sqrt(1. - x * x);
a /= x;
c = atan(a);
c += 4. * atan(1.0);
}
else if (x == 0.)
c = 2. * atan(1.0);
else {
a = sqrt(1. - x * x);
a /= x;
c = atan(a);
}
return(c);
}
/*-----------------------------------------------------------------------
function: cosh
function for hyperbolic cosine of x
-----------------------------------------------------------------------*/
/*
float cosh(x)
float x;
{
float c;
c = exp(x);
c += exp(-x);
c *= .5;
return(c);
}
*/
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.