This is ptrackfuns.c in view mode; [Download] [Up]
#include <stdio.h>
#include "ptrack.h"
#include "lpsf.h"
extern void report (char *s);
float ptable(fmin,fmax,tphi,tpsi,tgamph,tgamps,freq,n)
float fmin,fmax,tphi[50][5][MAXMM],tpsi[50][6][MAXMM],tgamph[50][5];
float tgamps[50][6],freq[50];
int n;
{
int i,j,k,np,nt;
static float phi[5][MAXMM],psi[6][MAXMM],gamphi[5],gampsi[6];
float omega,t,trigpo();
for (i=0 ; i<50 ; ++i)
{
freq[i] = fmin + ((float) (i)) * (fmax - fmin )/50.;
t = ( (float) (NYQ/10.) )/ freq[i];
nt = (int) t;
n = (nt < 5 ? nt : 5);
np=n+1;
omega = (freq[i]*(2.*PI))/ (SR/10.);
trigpo(omega,phi,psi,gampsi,gamphi,n);
for ( j = 0 ; j < n ; ++j)
{
for ( k= 0 ; k < MM ; ++k)
tphi[i][j][k] = phi[j][k];
tgamph[i][j] = gamphi[j];
}
for (j=0 ; j<np ; ++j)
{
for (k=0 ; k < MM ; ++k)
tpsi[i][j][k]=psi[j][k];
tgamps[i][j]=gampsi[j];
}
}
return;
}
float
trigpo(omega,phi,psi,gampsi,gamphi,n)
float omega,psi[6][MAXMM],phi[5][MAXMM];
float gamphi[5],gampsi[6];
int n;
{
int j=0;
int k,np;
double alpha,beta,gamma;
double wcos[MAXMM],wsin[MAXMM];
double p,z,a,b,yy;
double sin(),cos();
np=n+1;
for (k=0 ; k<MM ; ++k)
{
yy = omega* (float) (k);
wcos[k]=cos(yy);
wsin[k]=sin(yy);
}
beta=gamma= ZERO;
for (k=0 ; k<MM ; ++k)
{
p=wsin[k];
z=p*p;
beta+=z*wcos[k];
gamma+=z;
phi[0][k]=p;
}
gamphi[0]=gamma;
a=TWO * beta/gamma;
alpha=beta=gamma=ZERO;
for (k=0 ; k<MM ; ++k)
{
p= (TWO * wcos[k]-a) * phi[0][k];
alpha += wcos[k] * p * phi[0][k];
beta += wcos[k]* ( p * p ) ;
gamma += ( p * p );
phi[1][k] = p ;
}
gamphi[1]=gamma;
a = TWO * beta/gamma;
b = TWO *alpha/gamphi[0];
for (j=2 ; j<n ; ++j)
{
alpha=beta=gamma = ZERO;
for (k=0 ; k< MM ; ++k)
{
p = (TWO * wcos[k] - a ) * phi[j-1][k] - b * phi[j-2][k];
alpha += wcos[k] * p * phi[j-1][k];
beta += wcos[k] * (p * p);
gamma += (p * p) ;
phi[j][k] = p ;
}
gamphi[j] = gamma;
a = TWO * beta/gamma;
b = TWO *alpha/gamphi[j-1];
}
beta = ZERO ;
gamma = (double) MM;
for ( k=0; k < MM ; ++k)
{
beta += wcos[k];
psi[0][k]= ONE;
}
gampsi[0]=gamma;
a=beta/gamma;
alpha=beta=gamma= ZERO;
for ( k=0 ; k < MM ; ++k)
{
p = wcos[k]-a;
alpha += wcos[k] * p*psi[0][k];
beta += wcos[k] * ( p * p );
gamma += (p * p);
psi[1][k] = p ;
}
gampsi[1] = gamma ;
a = TWO * beta / gamma ;
b = TWO * alpha / gampsi[0];
for (j=2 ; j<np ;++j)
{
alpha=beta=gamma= ZERO;
for ( k=0; k < MM ; ++k)
{
p=(TWO * wcos[k]-a)* psi[j-1][k]-b*psi[j-2][k];
alpha += wcos[k]*p*psi[j-1][k];
beta += wcos[k]* (p*p);
gamma += (p*p);
psi[j][k] = p;
}
gampsi[j]=gamma;
a=TWO*beta/gamma;
b=TWO*alpha/gampsi[j-1];
}
return;
}
float
getrms(sig)
double *sig; /* from float */
{
int step;
float rms = ZERO ;
double sqrt();
for (step = 0; step < LSLICE; ++step) rms += sig[step]*sig[step];
return (sqrt ((double) (rms / (double) LSLICE)));
}
float
getpch(sigf,tphi,tpsi,tgamph,tgamps,freq,n)
float tphi[50][5][MAXMM],tpsi[50][6][MAXMM],tgamph[50][5],tgamps[50][5];
double *sigf;
float *freq;
int n;
{
int i, j,tenj;
float s[PFRAMAX/10],g[MAXMM],h[MAXMM];
float fm,xx,qsum,search();
for (j=0 ;j < JMAX ; ++j)
{
tenj = 10 * (j+1) - 1 ;
s[j]=sigf[tenj];
}
for (i=0 ; i<MM ; ++i)
{
g[i] = .5 * ( s[MM+i-1] - s[MM-i-1] );
h[i] = .5 * ( s[MM+i-1] + s[MM-i-1] );
}
qsum = 0.;
for (i=0 ; i<MM ; ++i)
qsum += (g[i]*g[i]) + (h[i]*h[i]);
xx=search(&fm,qsum,tphi,tpsi,tgamph,tgamps,g,h,freq,n);
return(xx);
}
float search(fm,qsum,tphi,tpsi,tgamph,tgamps,g,h,freq,n)
float *fm,qsum;
float tphi[50][5][MAXMM],tpsi[50][6][MAXMM],tgamph[50][5];
float tgamps[50][6],g[],h[],freq[];
int n;
{
float funmin, filfun();
float fun[50];
float f1,f2,f3,x0,x1,x2,x3,a,b,c,ftemp;
int istar,i;
filfun (tphi,tgamph,tpsi,tgamps,fun,freq,g,h,qsum,n);
funmin = fun[istar = 0];
for (i = 1; i < 50; ++i) {
ftemp = fun[i];
if(ftemp < funmin) {
funmin = ftemp;
istar = i;
}
}
if( istar == 0 || istar == 49 ) {
*fm = fun[istar];
return (freq[istar]);
}
else {
x1 = freq[istar-1];
f1 = fun[istar-1];
x2 = freq[istar];
f2 = fun[istar];
x3 = freq[istar+1];
f3 = fun[istar+1];
a = f3/((x3-x1)*(x3-x2));
b = f2/((x2-x1)*(x2-x3));
c = f1/((x1-x2)*(x1-x3));
x0 = .5*(a*(x1+x2)+b*(x1+x3)+c*(x2+x3))/(a+b+c);
*fm = a*(x0-x1)*(x0-x2)+b*(x0-x1)*(x0-x3)+c*(x0-x2)*(x0-x3);
return(x0);
}
}
float
filfun(tphi,tgamph,tpsi,tgamps,fun,freq,g,h,qsum,n)
/* Too many magic numbers */
float tphi[50][5][MAXMM],tpsi[50][6][MAXMM],tgamph[50][5];
float tgamps[50][6],*g,*h,*freq,*fun,qsum;
int n;
{
float c,sum;
int i,np,j,k;
for(i=0 ; i < 50 ;++i)
{
n = (NYQ/10.)/freq[i];
n = (n < 5) ? n : 5;
/* n = (int)((((float)(NYQ/10.))/freq[i]) < 5) ? n : 5; */
np = n+1;
sum = ZERO;
for(j=0 ; j < n ; ++j)
{
c = ZERO;
for(k=0 ; k< MM ; ++k)
c += g[k]*tphi[i][j][k];
sum += (c*c)/tgamph[i][j];
}
for (j=0 ; j<np ; ++j)
{
c = ZERO;
for (k=0 ; k < MM ; ++k)
c += h[k] * tpsi[i][j][k];
sum += (c*c)/tgamps[i][j];
}
fun[i] = qsum-sum;
}
return;
}
static double w1,w11,w12,w2,w21,w22,w3,w31,w32,w4,w41,w42;
void init_lowpass ()
{
w1 = w11 = w12 = w2 = w21 = w22 = w3 = w31 = w32 = w4 = w41
= w42 = ZERO;
return;
}
double lowpass(x)
double x; /* from float */
{
static double c= .00048175311, a1= -1.89919924, c1= -1.923248041,
d1= .985720370, a2= -1.86607670, c2= -1.90075003,
d2= .948444690, a3= -1.66423461, c3= -1.87516686,
d3= .896241023, c4= -.930449120;
double temp,y;
w1=c*x-c1*w11-d1*w12;
temp=w1+a1*w11+w12;
w12=w11;
w11=w1;
w2=temp-c2*w21-d2*w22;
temp=w2+a2*w21+w22;
w22=w21;
w21=w2;
w3=temp-c3*w31-d3*w32;
temp=w3+a3*w31+w32;
w32=w31;
w31=w3;
w4=temp-c4*w41;
y=w4+w41;
/* w42 set but not used in lowpass */
w42=w41;
w41=w4;
return(y);
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.