This is tools.c in view mode; [Download] [Up]
/* Standalone analysis program, reading fixed-point monaural sound files.
* Currently set for maximum of 34 poles, & max anal segment of 500 samples,
* meaning that the frame slices cannot exceed 250 samples.
* Program size expands linearly with slice size, and as square of npoles.
*/
#include "sysdep.h"
#include <stdio.h>
#include <sys/file.h>
#include <fcntl.h>
#include <math.h>
#include <pwd.h>
#define NEXT
#include "sfheader.h"
#include "lpsf.h"
#include "lp.h"
#define INCR 4
#define readin read
extern char soundFile[], lpFile[], *comment;
extern int npoles, lpframeoffset;
extern float lpinskip, lpdur;
extern void report(char *s);
static int NPOLE; /* why are these external? -dll */
static int FRAME;
static int NPP1;
static int alpol(), gauss(), dies(), die();
static int debug=0, verbose=0;
#define NDATA 4 /* number of data values stored with frame */
int anallpc ()
{
int jj, slice, counter;
float coef[POLEMAX+NDATA], inskip, dur;
char *sfname;
char *iname;
char *oname;
double errn, rms1, rms2, cc[40];
double sigi[FRAMAX];
long i, nsamps, skipbytes, outskip;
struct SFDESC sfdesc;
int anal, sfd, nbytes, nblpc, firstframe, nread;
int n, hsize, esize; /* new declarations */
char inbuf[SIZEOF_HEADER], framerpt[64];
SFHEADER *sfh;
int Hflag = 0;
LPHEADER *lph;
char *lpbuf, *tp;
char c;
int nframes;
lpbuf = (char *) calloc(LPBUFSIZ, sizeof(char));
lph = (LPHEADER *) lpbuf;
tp = lph->text;
iname = &soundFile[0]; /* '\0' ==> stdin */
oname = &lpFile[0]; /* '\0' ==> stdout */
NPOLE = npoles;
slice = lpframeoffset;
inskip = lpinskip;
dur = lpdur; /* should implement an EOF dur */
strncat(tp, comment, (LPBUFSIZ - sizeof(LPHEADER) + 4));
tp += strlen(tp);
if ((anal = creat (oname,0644)) < 0) {
dies ("cannot open output data file: %s", oname);
return 0;
}
sfname = iname;
if ((sfd = open (sfname,O_RDONLY)) < 0) {
dies("cannot open sound file %s", sfname);
close (anal);
return 0;
}
if ((n = readin (sfd, inbuf, sizeof (SNDSoundStruct)))
< sizeof (SNDSoundStruct)) {
die("soundfile header read error");
close (anal); close (sfd);
return 0;
}
sfh = (SFHEADER *) inbuf;
if (sfh->sfinfo.NeXTheader.magic != SND_MAGIC) {
/* if (sfmagic (sfh) != SF_MAGIC) { */
dies("%s is not a soundfile", iname);
close (anal); close (sfd);
return 0;
}
/* sfdesc.sf_chans = sfchans(sfh);
sfdesc.sf_class = sfclass(sfh);
sfdesc.sf_srate = sfsrate(sfh);*/
sfdesc.sf_chans = sfh->sfinfo.NeXTheader.channelCount;
switch (sfh->sfinfo.NeXTheader.dataFormat) {
case SND_FORMAT_LINEAR_16:
sfdesc.sf_class = SF_SHORT;
break;
case SND_FORMAT_FLOAT:
sfdesc.sf_class = SF_FLOAT;
break;
default:
die ("data format must be short or float");
close (anal); close (sfd);
return 0;
}
sfdesc.sf_srate = sfh->sfinfo.NeXTheader.samplingRate;
if (debug) {
fprintf(stderr, "sf_chans=%d sf_srate=%d sf_class=%d\n",
sfdesc.sf_chans, sfdesc.sf_srate, sfdesc.sf_class);
}
if(sfdesc.sf_chans != 1) {
die("takes monaural files only");
close (anal); close (sfd);
return 0;
}
if (NPOLE > POLEMAX) {
dies("poles exceeds maximum allowed");
close (anal); close (sfd);
return 0;
}
FRAME = slice * 2;
if (FRAME > FRAMAX) {
die("framesize (inter-frame-offset*2) exceeds maximum allowed");
close (anal); close (sfd);
return 0;
}
lph->lpmagic = LP_MAGIC;
lph->srate = sfdesc.sf_srate;
lph->npoles = NPOLE;
lph->nvals = NPOLE + 4;
lph->framrate = lph->srate / (float) slice;
lph->duration = dur;
/* compute header size */
hsize = tp - (char *) lph; /* sizeof text */
esize = INCR - (hsize % INCR); /* round up to 4 byte boundary */
if (esize == 4) esize = 0; /* if we need it */
lph->headersize = hsize + esize;
if (lph->headersize >= (LPBUFSIZ - sizeof(LPHEADER) + 4))
lph->headersize = LPBUFSIZ;
if (write(anal, (char *) lph, lph->headersize) < lph->headersize) {
die("anallpc: can't write header");
close (anal); close (sfd);
return 0;
}
lseek (sfd, sfh->sfinfo.NeXTheader.dataLocation, L_SET);
skipbytes = (long)(inskip * sfdesc.sf_srate) * (long)sfdesc.sf_class;
if (lseek (sfd, skipbytes, L_INCR) < 0) {
die ("bad sflseek on skip");
close (anal); close (sfd);
return 0;
}
/* fprintf(stderr,"Specify starting output frame number\t");
/* scanf("%d",&firstframe);
/* for now we'll just say to start at 0 ('C' mind set) */
firstframe = 0;
NPP1 = NPOLE+1;
nblpc = (NPOLE + NDATA)*FLOAT;
outskip = firstframe * nblpc + lph->headersize; /* headersize is new */
if ((lseek (anal,outskip,0)) < 0) {
die("Bad lseek on analysis file");
close (anal); close (sfd);
return 0;
}
nsamps = dur * sfdesc.sf_srate;
nbytes = FRAME * sfdesc.sf_class; /* bytes for whole frame */
if ((nread = readin(sfd, (char *)sigi, nbytes)) != nbytes) {
die("soundfile read error: couldn't fill first frame");
close (anal); close (sfd);
return 0;
}
if (sfdesc.sf_class == SF_SHORT) {
short *pt = (short *) &sigi[0];
for (i = FRAME - 1; i >= 0; i--) sigi[i] = (double) pt[i];
}
else { /* SF_FLOAT */
float *pt = (float *) &sigi[0];
for (i = FRAME - 1; i >= 0; i--) sigi[i] = (double) pt[i];
}
/* for(i=0;i<FRAME;i++) printf("%d ",sigi[i]); */
nbytes >>= 1; /* bytes to refill halfframe */
counter = 1;
nframes = nsamps / slice;
for (i = 0; i < nsamps; ) {
if (!alpol(sigi,&errn,&rms1,&rms2,cc)) {
close (anal); close (sfd);
return 0;
}
coef[0] = (float)rms2;
coef[1] = (float)rms1;
coef[2] = (float)errn;
coef[3] = 0.; /* for pitch of this frame */
/* reverse the coefs & change sign */
for(jj=NDATA; jj<NPOLE+NDATA; jj++)
coef[jj] = (float)-cc[NPOLE-jj+(NDATA - 1)];
/* then write the anal frame */
if ((n = write(anal, (char *)coef, nblpc)) != nblpc) {
die("write error");
close (anal); close (sfd);
return 0;
}
sprintf (framerpt, "frame %5d of %5d", counter++, nframes);
report (framerpt);
for(jj=0; jj<slice; jj++,i++) /* now move slice forward */
sigi[jj] = sigi[jj+slice]; /* & refill: */
for (jj = slice; jj < FRAME; jj++) sigi[jj] = 0.0;
if ((n = readin(sfd, (char *)(sigi+slice), nbytes)) <= 0)
break; /* ran up against eof */
if (sfdesc.sf_class == SF_SHORT) {
short *pt = (short *) &sigi[slice];
for (jj = FRAME - 1; jj >= slice; jj--)
sigi[jj] = (double) pt[jj - slice];
}
else { /* SF_FLOAT */
float *pt = (float *) &sigi[slice];
for (jj = FRAME - 1; jj >= slice; jj--)
sigi[jj] = (double) pt[jj - slice];
}
}
close (anal); close (sfd);
return 1;
}
static int alpol(sig, errn, rms1, rms2, c)
double *sig, *errn, *rms1, *rms2, *c;
{
double a[POLEMAX][POLEMAX], v[POLEMAX], b[POLEMAX];
double *vp=v, *bp=b;
double sum, sumx, sumy;
int k1, i, l, k, limit, j;
k1 = NPOLE + 1;
for (i = 0; i < NPOLE; ++i) {
sum = (double) 0.0;
for (k = NPOLE; k < FRAME; ++k) sum += sig[k-(i+1)] * sig[k];
v[i] = -sum;
if (i != NPOLE - 1) {
limit = NPOLE - (i+1);
for (l=0; l < limit ;++l) {
sum += sig[NPOLE-(i+1)-(l+1)] *
sig[NPOLE-(l+1)] -
sig[FRAME-(i+1)-(l+1)] *
sig[FRAME-(l+1)];
a[(i+1)+l][l] = a[l][(i+1)+l] = sum;
}
}
}
sum = (double) 0.0;
for (k = NPOLE; k < FRAME; ++k) sum += sig[k] * sig[k];
sumy = sumx = sum;
for (l = 0; l < NPOLE; ++l) {
sum += sig[NPOLE-(l+1)] * sig[NPOLE-(l+1)] -
sig[FRAME-(l+1)] * sig[FRAME-(l+1)];
a[l][l] = sum;
}
if (!gauss (a, v, b)) return 0;
for (i = 0; i < NPOLE; ++i) sumy = sumy - b[i] * v[i];
*rms1 = sqrt(sumx/(double) (FRAME - k1 + 1) );
*rms2 = sqrt(sumy/(double) (FRAME - k1 + 1) );
*errn = (*rms2 * *rms2) / (*rms1 * *rms1);
for (bp = b; bp - b < NPOLE; ++bp, ++c) *c = *bp;
return 1;
}
static int gauss(aold, bold, b)
double aold[POLEMAX][POLEMAX], *bold, *b;
{
double amax, dum, pivot;
double c[POLEMAX], a[POLEMAX][POLEMAX];
int i, j, k, l, istar, ii, lp;
/* aold and bold untouched by this subroutine */
for (i=0; i < NPOLE ;++i) {
c[i] = bold[i];
for (j = 0; j < NPOLE; ++j) a[i][j] = aold[i][j];
}
/* eliminate i-th unknown */
for (i = 0; i < NPOLE - 1; ++i) { /* find largest pivot */
amax = 0.0;
for (ii = i; ii < NPOLE; ++ii) {
if (fabs (a[ii][i]) >= amax) {
istar = ii;
amax = fabs (a[ii][i]);
}
}
if (amax < 1e-20) {
die("gauss: ill-conditioned");
return 0;
}
for (j = 0; j < NPOLE; ++j) { /* switch rows */
dum = a[istar][j];
a[istar][j] = a[i][j];
a[i][j] = dum;
}
dum = c[istar];
c[istar] = c[i];
c[i] = dum;
/* pivot */
for (j = i + 1; j < NPOLE; ++j) {
pivot = a[j][i] / a[i][i];
c[j] = c[j] - pivot * c[i];
for (k = 0; k < NPOLE; ++k)
a[j][k] = a[j][k] - pivot * a[i][k];
}
} /* return if last pivot is too small */
if (fabs (a[NPOLE-1][NPOLE-1]) < 1e-20) {
die("gauss: ill-conditioned");
return 0;
}
*(b+NPOLE-1) = c[NPOLE-1] / a[NPOLE-1][NPOLE-1];
/* back substitute */
for (k = 0; k < NPOLE - 1; ++k) {
l = NPOLE-1 -(k+1);
*(b+l) = c[l];
lp = l + 1;
for (j = lp; j < NPOLE; ++j) *(b+l) += -a[l][j] * *(b+j);
*(b+l) /= a[l][l];
}
return(1);
}
static char errmsg[60];
static int dies(s,t)
char *s, *t;
{
sprintf(errmsg,s,t);
return die(errmsg);
}
static int die(s)
char *s;
{
report (s);
return 0;
}
#define PFRAMSIZ 2 /* words per frame (pitch&rms */
#define PCHSIZE 16 /* record size in float words */
#define BPPFRAME FLOAT*PFRAMSIZ /* bytes per frame */
#define BPPREC BPPFRAME*PCHSIZE /* bytes per record */
#define FPPREC PCHSIZE/PFRAMSIZ /* frames per record */
static float PFPS; /* pitch frames per second */
static float sampspf;
static int pitches;
static int lpclast;
static float getfreq();
extern char ptFile[];
extern float ptframerate;
int lpconcat ()
{
int anal,j;
int npoles,lpcframe,pchframe,pchlast,nbpch;
long nskiplpc,nskippch,nblpc;
char *input, *output;
float pch[2],val,sr,time;
LPHEADER *lph;
char lpbuf[LPBUFSIZ];
char c;
input = &ptFile[0];
output = &lpFile[0];
PFPS = ptframerate; /* this is ptrack.c default SR/JSLIDE */
lpcframe = 0;
lpclast = -1; /* to EOF */
if((anal = open (output,2)) < 0) {
report ("Can't open lpc analysis file");
return 0;
}
if((pitches = open (input,0)) < 0) {
report ("Can't open pitch analysis file");
return 0;
}
if ((j = read (anal, lpbuf, LPBUFSIZ)) <= 0) {
report ("Can't read header on lpc analysis file");
return 0;
}
lph = (LPHEADER *) lpbuf;
npoles = lph->npoles;
sr = lph->srate;
sampspf = sr / lph->framrate;
nblpc = (npoles+4)*FLOAT - FLOAT;/*to beginning of next pchloc*/
nskiplpc = (long)(lpcframe)*(long)(nblpc+FLOAT) + (3*FLOAT);
/* pch is 4th data value in fr*/
if((lseek (anal,nskiplpc + lph->headersize,0)) < 0) { /* add */
report ("Bad lseek on analysis file"); /* headersize */
return 0;
}
for(j = lpcframe; ;j++) {
time = (float)j*sampspf/sr;
if ((j > lpclast && lpclast != -1) || (time > lph->duration))
break;
/*
* if((read(pitches,(char *)pch,nbpch)) != nbpch) {
* printf("Bad read on pitch analysis file\n");
* exit(1);
* }
*/
if ((val = getfreq (time)) < 0.0) return 0;
if((write(anal,(char *)&val,FLOAT)) != FLOAT) {
report ("Bad write on lpc analysis file");
return 0;
}
lseek (anal,nblpc,1);
}
return 1;
}
static float getfreq(time)
float time;
{
int i,j,frame;
float fframe,fraction,cps,error;
static int oldframe = 0;
static int endframe = 0;
static float parray[PCHSIZE*2];
fframe = time * PFPS;
frame = (int)fframe;
fraction = fframe - (float)frame;
if(!((frame >= oldframe) && (frame < endframe))) {
if(lseek(pitches,((long)frame*(long)BPPFRAME),0) < 0) {
report ("Bad lseek on pitch file");
return -1.0;
}
read (pitches,(char *)parray,BPPREC);
oldframe = frame;
endframe = oldframe + FPPREC - 1;
}
cps = (1.-fraction) * *(parray+(frame-oldframe)*PFRAMSIZ) +
fraction * *((parray+(frame-oldframe)*PFRAMSIZ)+2);
error = (1.-fraction) * *(parray+(frame-oldframe)*PFRAMSIZ+1) +
fraction * *((parray+(frame-oldframe)*PFRAMSIZ)+3);
return(cps);
}
/* PTRACK - Standalone pitch analysis program. Revised dll@ems.
* Reads a fixed-point monaural sound file (with or without header).
*
* Program accepts command-line data (see man entry),
* or with -q will querie the user for the following:
*
* -Name of soundfile.
* -Sampling rate (if no header).
* -Name of file to contain analysis data--can be an old or null file.
* -Framesize. The size of the frame to be analyzed, in samples.
* Maximum is 350, which has been good working number. Frame should
* be able to contain at least one expected pitch period of signal.
* -JSLIDE. The number of new samples to add per frame. This number
* should normally be less than half of LSLICE. 125 is good for 15k
* and 250 is good for 30k sampling rate.
* This number reflects the real time period of the frame, i.e. how
* far you creep up for each frame in the signal.
* -pchlow and pchigh-- expected lower and upper bounds of analysis.
* the better you guess at this, the better the quantization will
* be in the analysis.
* At the present time the highest fundamental analyzable is about SR/20
* cps, due to the settings of the lowpass filter. This can be fixed
* by using interpolation, as in the fortran version of this program,
* or, more sensibly, different filters.
* -inskip, amount of time to skip on input file before beginning analysis.
* The program is now written to write data sequentially beginning at the
* start of the data file, but this can be changed by simple adding an
* lseek.
* (inskip should match the -s value in soundout).
* -dur, amount of time to analyze. This should equal, but not
* exceed the -d value from soundout.
* The following is good for a 15k signal, for example:
*
* inname 15000 outname
* 350 125
* 200 300
* 0 .5
*
*/
#include "ptrack.h"
extern float ptinskip, ptdur, sr, ptlowest, pthighest;
extern int ptframesize, ptframeoffset;
extern double lowpass();
extern void init_lowpass();
int LSLICE, JSLIDE, JMAX, MM;
float SR, NYQ;
static float gtphi[50][5][MAXMM],gtpsi[50][6][MAXMM];
static float gtgamph[50][5],gtgamps[50][6];
int ptrack ()
{
double sig[PFRAMAX];
int MIDPOINT,needed,howmuch;
float freq[50];
float getpch(),getrms(), ptable();
float pchlow,pchigh,input[4];
FILE *pdata;
int i,n,jj,anal,framesize,sampsize;
long nsamps;
float data[2],inskip,dur;
char *name, *fullsfname, ptrpt[64];
char *sigp,*sigtemp;
double rmsm,freqm,rmsx,freqx;
SFHEADER *sfh;
char sfbuf[SIZEOF_HEADER];
char *sfname;
int sfd, sfn, counter;
int Hflag = 0;
char c;
int rflag = 0;
sfname = &soundFile[0]; /* '\0' => stdin */
name = &ptFile[0]; /* '\0' => stdout */
LSLICE = ptframesize; /* 350 is max! */
JSLIDE = ptframeoffset;
pchlow = ptlowest;
pchigh = pthighest;
inskip = ptinskip;
dur = ptdur; /* should implement an EOF dur! */
SR = sr;
if ((anal = creat (name,0644)) < 0) { /* create if not exists */
dies("ptrack: can't open output data file %s", name);
return 0;
}
fullsfname = sfname;
if ((sfd = open (fullsfname,O_RDONLY)) < 0) {
dies("Ptrack: can't open soundfile %s",fullsfname);
close (anal);
return 0;
}
if ((sfn = readin (sfd,sfbuf,sizeof (SNDSoundStruct)))
< sizeof (SNDSoundStruct)) {
die("Ptrack: error reading soundfile header");
close (anal);
close (sfd);
return 0;
}
sfh = (SFHEADER *) sfbuf;
/* if (!ismagic(sfh)) {*/
if (sfh->sfinfo.NeXTheader.magic != SND_MAGIC) {
dies("Ptrack: %s is not a soundfile",sfname);
close (anal);
close (sfd);
return 0;
}
switch (sfh->sfinfo.NeXTheader.dataFormat) {
case SND_FORMAT_LINEAR_16:
sampsize = SF_SHORT;
break;
case SND_FORMAT_FLOAT:
sampsize = SF_FLOAT;
break;
default:
die ("data format must be short or float");
close (anal); close (sfd);
return 0;
}
{ /* do input skipping */
long skipbytes;
int nread;
lseek (sfd, sfh->sfinfo.NeXTheader.dataLocation, L_SET);
skipbytes = (long) (inskip * SR * sampsize);
if (lseek (sfd,skipbytes,L_INCR) < 0) {
die("bad sflseek on skip");
close (anal);
close (sfd);
return 0;
}
}
NYQ = SR/2.;
JMAX = LSLICE/10;
MM = ((LSLICE/10+1)/2);
MIDPOINT = LSLICE - JSLIDE;
sigp = (char *)(sig+MIDPOINT);
ptable(pchlow,pchigh,gtphi,gtpsi,gtgamph,gtgamps,freq,n);
if ((jj = readin(sfd,(char *)sig,LSLICE*sampsize))
!= LSLICE*sampsize) {
die("Ptrack: couldn't fill first frame");
close (anal);
close (sfd);
return 0;
}
if (sampsize == SF_SHORT) {
short *pt = (short *) &sig[0];
for (i = LSLICE - 1; i >= 0; i--) sig[i] = (double) pt[i];
}
else { /* SF_FLOAT */
float *pt = (float *) &sig[0];
for (i = LSLICE - 1; i >= 0; i--) sig[i] = (double) pt[i];
}
init_lowpass ();
for (i = 0; i < LSLICE; i++) sig[i] = lowpass (sig[i]);
nsamps = SR * dur;
counter = 1;
while (nsamps) {
data[0] = getpch(sig,gtphi,gtpsi,gtgamph,gtgamps,freq,n);
data[1] = getrms(sig);
write (anal,(char *)data,8);
sprintf (ptrpt, "frame %5d, pitch = %8.2f", counter++,
data[0]);
report (ptrpt);
sigp = (char *)(sig+MIDPOINT);
for (i = 0; i < MIDPOINT; i++) sig[i] = sig[i+JSLIDE];
howmuch = 0;
needed = JSLIDE * sampsize;
if ((howmuch = readin (sfd,sigp,needed)) != needed) break;
if (sampsize == SF_SHORT) {
short *pt = (short *) &sig[MIDPOINT];
for (jj = LSLICE - 1; jj >= MIDPOINT; jj--)
sig[jj] = (double) pt[jj - MIDPOINT];
}
else { /* SF_FLOAT */
float *pt = (float *) &sig[MIDPOINT];
for (jj = LSLICE - 1; jj >= MIDPOINT; jj--)
sig[jj] = (double) pt[jj - MIDPOINT];
}
for(i = MIDPOINT; i < LSLICE; i++) sig[i] = lowpass (sig[i]);
nsamps -= JSLIDE;
}
close (anal); close (sfd);
return 1;
}
void lpdump (FILE *outfile)
{
int fd, nbytes, l, framsize, i, frmno = 0;
char lpbuf[LPBUFSIZ], dumprpt[32];
float frmbuf[MAXFRAME];
LPHEADER *lph;
if ((fd = open (lpFile, O_RDONLY)) < 0) {
dies ("lpdump: can't open %s", lpFile);
return;
}
if ((nbytes = read (fd, lpbuf, LPBUFSIZ)) <= 0) {
dies ("lpdump: read error on %s", lpFile);
close (fd);
return;
}
lph = (LPHEADER *) lpbuf;
if (lph->lpmagic != LP_MAGIC) {
dies ("lpdump: %s is not an lpfile", lpFile);
close (fd);
return;
}
if ((l = lseek (fd, lph->headersize, 0)) < 0) {
dies ("lpdump: bad lseek on %s", lpFile);
close (fd);
return;
}
framsize = lph->nvals * FLOAT;
while (++frmno) {
if ((read (fd, (char *) frmbuf, framsize)) < framsize) {
die ("End of lp file");
close (fd);
return;
}
fprintf(outfile, "\n\nFrame #%d\n",frmno);
sprintf (dumprpt, "dumping frame %d", frmno);
report (dumprpt);
fprintf(outfile, "rms1: %f rms2: %f error: %f cps: %f\n",frmbuf[0], frmbuf[1], frmbuf[2], frmbuf[3]);
fprintf(outfile, "Coeffs:");
for (i = 0; i < lph->npoles; i++) {
if ((i % 6) == 0) fprintf (outfile, "\n ");
fprintf (outfile, "%10.6f ", frmbuf[4 + i]);
}
fprintf (outfile, "\n");
}
}
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.