On Sun, 22 Jan 2023, Sergei Golovan wrote:
Hi Bastian,
Thank you for the report!
On Mon, Jan 16, 2023 at 9:33 PM Bastian Germann <b...@debian.org> wrote:
While the file part that is copied from formant.c is okay to be included in
Debian, the part from downsample.c is not.
Please repack the source package getting rid of the file.
I don't think that I can remove this code without reimplementing its
functionality somehow. So maybe it'd be better to drop the package
from Debian. There are two packages that depend on tcl-snack:
transcriber and wavesurfer. I guess they'd have to be removed from
Debian as well.
Cheers!
--
Sergei Golovan
Dear maintainers,
I would very much appreciate if this package would be kept in Debian.
That is why I've had a look into the source code and have seen that there
is yet another implementation for a downsampler in jkGetF0.c that is not
proprietary code. And so I was thinking what is good enough to be used
for the estimation of the first formant F0, cannot be bad for doing the
computation to estimate all the others. However, in contrast to the code
from downsample.c, this other downsampler is not using 16-bit integer
values but it is a floating-point implementation. So it is not a drop-in
replacement and I had to rewrite large portions of the two functions
Fdownsample() and highpass(). I also had to make the function do_ffir()
public to be used by the highpass filter implementation and to uncomment
and fix some code in it that was headlined by /* never used */ but is
needed to make the highpass filter implementation working again.
Long story short, there is a patch attached to this e-mail that I grant
under the license of either GPL2+ (like it is used for the whole project
snack) or BSD (like it is said to be in the headers of jkFormant.c and
jkGetF0.c), whatever suits you, in the hope to save this package from
being dropped.
I have attached two versions of it because I don't know whether the 170
lines of conflicting proprietary code have to be removed before or can be
replaced by the patch. So applying remove_downsample.c.diff removes the
conflicting lines from jkFormant.c, formant.patch is the actual fix to
make the implementation work again, and formant2.patch is a combination
of both.
I hope it is not too late.
Best regards,
Thomas Uhle
--- a/generic/jkFormant.c
+++ b/generic/jkFormant.c
@@ -604,175 +604,6 @@
return(TRUE);
}
-/* Copyright (c) 1987, 1988, 1989 AT&T */
-/* All Rights Reserved */
-
-/* THIS IS UNPUBLISHED PROPRIETARY SOURCE CODE OF AT&T */
-/* The copyright notice above does not evidence any */
-/* actual or intended publication of such source code. */
-
-/* downsample.c */
-/* a quick and dirty downsampler */
-
-#ifndef TRUE
-# define TRUE 1
-# define FALSE 0
-#endif
-
-#define PI 3.1415927
-
-/* ---------------------------------------------------------- */
-int lc_lin_fir(fc,nf,coef)
-/* create the coefficients for a symmetric FIR lowpass filter using the
- window technique with a Hanning window. */
-register double fc;
-double coef[];
-int *nf;
-{
- register int i, n;
- register double twopi, fn, c;
-
- if(((*nf % 2) != 1) || (*nf > 127)) {
- if(*nf <= 126) *nf = *nf + 1;
- else *nf = 127;
- }
- n = (*nf + 1)/2;
-
- /* compute part of the ideal impulse response */
- twopi = PI * 2.0;
- coef[0] = 2.0 * fc;
- c = PI;
- fn = twopi * fc;
- for(i=1;i < n; i++) coef[i] = sin(i * fn)/(c * i);
-
- /* Now apply a Hanning window to the (infinite) impulse response. */
- fn = twopi/((double)(*nf - 1));
- for(i=0;i<n;i++)
- coef[i] *= (.5 + (.5 * cos(fn * ((double)i))));
-
- return(TRUE);
-}
-
-/* ---------------------------------------------------------- */
-
-void do_fir(buf,in_samps,bufo,ncoef,ic,invert)
-/* ic contains 1/2 the coefficients of a symmetric FIR filter with unity
- passband gain. This filter is convolved with the signal in buf.
- The output is placed in buf2. If invert != 0, the filter magnitude
- response will be inverted. */
-short *buf, *bufo, ic[];
-int in_samps, ncoef, invert;
-{
- register short *buft, *bufp, *bufp2, stem;
- short co[256], mem[256];
- register int i, j, k, l, m, sum, integral;
-
- for(i=ncoef-1, bufp=ic+ncoef-1, bufp2=co, buft = co+((ncoef-1)*2),
- integral = 0; i-- > 0; )
- if(!invert) *buft-- = *bufp2++ = *bufp--;
- else {
- integral += (stem = *bufp--);
- *buft-- = *bufp2++ = -stem;
- }
- if(!invert) *buft-- = *bufp2++ = *bufp--; /* point of symmetry */
- else {
- integral *= 2;
- integral += *bufp;
- *buft-- = integral - *bufp;
- }
-/* for(i=(ncoef*2)-2; i >= 0; i--) printf("\n%4d%7d",i,co[i]); */
- for(i=ncoef-1, buft=mem; i-- > 0; ) *buft++ = 0;
- for(i=ncoef; i-- > 0; ) *buft++ = *buf++;
- l = 16384;
- m = 15;
- k = (ncoef << 1) -1;
- for(i=in_samps-ncoef; i-- > 0; ) {
- for(j=k, buft=mem, bufp=co, bufp2=mem+1, sum = 0; j-- > 0;
- *buft++ = *bufp2++ )
- sum += (((*bufp++ * *buft) + l) >> m);
-
- *--buft = *buf++; /* new data to memory */
- *bufo++ = sum;
- }
- for(i=ncoef; i-- > 0; ) { /* pad data end with zeros */
- for(j=k, buft=mem, bufp=co, bufp2=mem+1, sum = 0; j-- > 0;
- *buft++ = *bufp2++ )
- sum += (((*bufp++ * *buft) + l) >> m);
- *--buft = 0;
- *bufo++ = sum;
- }
-}
-
-/* ******************************************************************** */
-
-int get_abs_maximum(d,n)
- register short *d;
- register int n;
-{
- register int i;
- register short amax, t;
-
- if((t = *d++) >= 0) amax = t;
- else amax = -t;
-
- for(i = n-1; i-- > 0; ) {
- if((t = *d++) > amax) amax = t;
- else {
- if(-t > amax) amax = -t;
- }
- }
- return((int)amax);
-}
-
-/* ******************************************************************** */
-
-int dwnsamp(buf,in_samps,buf2,out_samps,insert,decimate,ncoef,ic,smin,smax)
- short *buf, **buf2;
- int in_samps, *out_samps, insert, decimate, ncoef, *smin, *smax;
- short ic[];
-{
- register short *bufp, *bufp2;
- short *buft;
- register int i, j, k, l, m;
- int imax, imin;
-
- if(!(*buf2 = buft = (short*)ckalloc(sizeof(short)*insert*in_samps))) {
- perror("ckalloc() in dwnsamp()");
- return(FALSE);
- }
-
- k = imax = get_abs_maximum(buf,in_samps);
- if (k == 0) k = 1;
- if(insert > 1) k = (32767 * 32767)/k; /* prepare to scale data */
- else k = (16384 * 32767)/k;
- l = 16384;
- m = 15;
-
-
- /* Insert zero samples to boost the sampling frequency and scale the
- signal to maintain maximum precision. */
- for(i=0, bufp=buft, bufp2=buf; i < in_samps; i++) {
- *bufp++ = ((k * (*bufp2++)) + l) >> m ;
- for(j=1; j < insert; j++) *bufp++ = 0;
- }
-
- do_fir(buft,in_samps*insert,buft,ncoef,ic,0);
-
- /* Finally, decimate and return the downsampled signal. */
- *out_samps = j = (in_samps * insert)/decimate;
- k = decimate;
- for(i=0, bufp=buft, imax = imin = *bufp; i < j; bufp += k,i++) {
- *buft++ = *bufp;
- if(imax < *bufp) imax = *bufp;
- else
- if(imin > *bufp) imin = *bufp;
- }
- *smin = imin;
- *smax = imax;
- *buf2 = (short*)ckrealloc((void *) *buf2,sizeof(short) * (*out_samps));
- return(TRUE);
-}
-
/* ---------------------------------------------------------- */
int ratprx(a,k,l,qlim)
Description: Replace dwnsample() and do_fir() by downsample() and do_ffir() resp.
Utilize the functions downsample() and do_ffir() from jkGetF0.c for the signal
pre-processing. The former functions are from proprietary code and, thus, can't
be used. Notice that do_fir() was a 16-bit integer implementation for an FIR
filter whereas do_ffir() is a floating-point implementation. This might cause a
higher computational effort. On the other side we gain precision. That is why
the volume of the audio signal can stay the same now and we experience much less
effects from quantization.
Author: Thomas Uhle <thomas.u...@mailbox.tu-dresden.de>
Last-Modified: Sat, 18 Feb 2023 20:33:38 +0100
Bug-Debian: https://bugs.debian.org/1029034
--- a/generic/jkFormant.c
+++ b/generic/jkFormant.c
@@ -611,13 +611,11 @@
int *l, *k, qlim;
{
double aa, af, q, em, qq = 0, pp = 0, ps, e;
- int ai, ip, i;
+ int ai, ip, result = FALSE;
aa = fabs(a);
ai = (int) aa;
-/* af = fmod(aa,1.0); */
- i = (int) aa;
- af = aa - i;
+ af = aa - ai;
q = 0;
em = 1.0;
while(++q <= qlim) {
@@ -628,97 +626,121 @@
em = e;
pp = ip;
qq = q;
+ result = TRUE;
}
};
*k = (int) ((ai * qq) + pp);
*k = (a > 0)? *k : -(*k);
*l = (int) qq;
- return(TRUE);
+ return(result);
}
/* ----------------------------------------------------------------------- */
+extern float *downsample();
+
Sound *Fdownsample(s,freq2,start,end)
double freq2;
Sound *s;
int start;
int end;
{
- short *bufin, **bufout;
- static double beta = 0.0, b[256];
- double ratio_t, maxi, ratio, beta_new, freq1;
- static int ncoeff = 127, ncoefft = 0, nbits = 15;
- static short ic[256];
- int insert, decimate, out_samps, smin, smax;
+ float *bufin, *bufout, *bufp;
+ int frame_size = 1024, act_size, first_time, last_time;
+ double ratio, freq1;
+ int ncoeff, insert, decimate, total_samps, out_samps, ndone;
Sound *so;
register int i, j;
freq1 = s->samprate;
-
- if((bufout = (short**)ckalloc(sizeof(short*)))) {
- bufin = (short *) ckalloc(sizeof(short) * (end - start + 1));
- for (i = start; i <= end; i++) {
- bufin[i-start] = (short) Snack_GetSample(s, 0, i);
- }
+ ratio = freq2/freq1;
- ratio = freq2/freq1;
- ratprx(ratio,&insert,&decimate,10);
- ratio_t = ((double)insert)/((double)decimate);
+ if (!ratprx(ratio, &insert, &decimate, 10))
+ return(NULL);
- if(ratio_t > .99) return(s);
-
- freq2 = ratio_t * freq1;
- beta_new = (.5 * freq2)/(insert * freq1);
+ if (decimate <= insert)
+ return(NULL);
+
+ freq1 *= (double)insert;
+ freq2 = freq1/((double)decimate);
- if(beta != beta_new){
- beta = beta_new;
- if( !lc_lin_fir(beta,&ncoeff,b)) {
- printf("\nProblems computing interpolation filter\n");
- return(FALSE);
- }
- maxi = (1 << nbits) - 1;
- j = (ncoeff/2) + 1;
- for(ncoefft = 0, i=0; i < j; i++){
- ic[i] = (int) (0.5 + (maxi * b[i]));
- if(ic[i]) ncoefft = i+1;
+ /* filter length used in downsample(): 5ms */
+ ncoeff = ((int)(freq1 * 0.005))/2 + 1;
+
+ total_samps = (end - start + 1) * insert;
+ if (total_samps < (ncoeff * decimate * 3)) /* signal too short */
+ return(NULL);
+
+ if ((bufin = (float *) ckalloc(sizeof(float) * total_samps))) {
+ for (bufp = bufin, i = start; i <= end; i++) {
+ *bufp++ = Snack_GetSample(s, 0, i) * ((float)insert);
+ for(j = 1; j < insert; j++)
+ *bufp++ = 0.0f; /* insert zeros to boost the sampling frequency */
+ }
+
+ if ((frame_size * 2) > total_samps)
+ frame_size = (total_samps + 1)/2;
+
+ frame_size -= frame_size % decimate;
+
+ first_time = 1; /* new filter coefficients need to be computed */
+
+ for (ndone = 0, last_time = 0; !last_time; ndone += act_size) {
+ act_size = total_samps - ncoeff - ndone;
+ if (act_size > frame_size) {
+ act_size = frame_size;
+ out_samps = act_size/decimate;
+ } else {
+ out_samps = act_size/decimate;
+ if (!first_time && ((act_size + ncoeff) <= frame_size)) {
+ act_size += ncoeff;
+ last_time = 1;
+ } else
+ act_size = out_samps * decimate;
}
- } /* endif new coefficients need to be computed */
- if(dwnsamp(bufin,end-start+1,bufout,&out_samps,insert,decimate,ncoefft,ic,
- &smin,&smax)){
- /* so->buff_size = so->file_size = out_samps;*/
- so = Snack_NewSound(0, LIN16, s->nchannels);
- Snack_ResizeSoundStorage(so, out_samps);
- for (i = 0; i < out_samps; i++) {
- Snack_SetSample(so, 0, i, (float)(*bufout)[i]);
+ if ((bufout = downsample(bufin+ndone, total_samps-ndone, act_size, freq1,
+ &out_samps, decimate, first_time, last_time))) {
+ if (first_time) {
+ first_time = 0;
+ so = Snack_NewSound((int)freq2, LIN16, s->nchannels);
+ if (!so) {
+ printf("Can't create a new Signal in downsample()\n");
+ break;
+ }
+ Snack_ResizeSoundStorage(so, total_samps/decimate);
+ so->length = 0;
+ }
+
+ Snack_PutSoundData(so, so->length, bufout, out_samps);
+ so->length += out_samps;
+ } else {
+ printf("Problems in downsample()\n");
+ break;
}
- so->length = out_samps;
- so->samprate = (int)freq2;
- ckfree((void *)*bufout);
- ckfree((void *)bufout);
- ckfree((void *)bufin);
- return(so);
- } else
- printf("Problems in dwnsamp() in downsample()\n");
+ }
+ ckfree((void *)bufin);
} else
printf("Can't create a new Signal in downsample()\n");
- return(NULL);
+ return(so);
}
/* ---------------------------------------------------------- */
-Sound
-*highpass(s)
+extern void do_ffir();
+
+Sound *highpass(s)
Sound *s;
{
- short *datain, *dataout;
- static short *lcf;
+ float *datain, *dataout;
+ static float *lcf;
static int len = 0;
double scale, fn;
register int i;
+ int frame_size = 1024, act_size, total_samps, out_samps, ndone, init;
Sound *so;
/* Header *h, *dup_header();*/
@@ -727,28 +749,62 @@
/* This assumes the sampling frequency is 10kHz and that the FIR
is a Hanning function of (LCSIZ/10)ms duration. */
- datain = (short *) ckalloc(sizeof(short) * s->length);
- dataout = (short *) ckalloc(sizeof(short) * s->length);
- for (i = 0; i < Snack_GetLength(s); i++) {
- datain[i] = (short) Snack_GetSample(s, 0, i);
- }
-
if(!len) { /* need to create a Hanning FIR? */
- lcf = (short*)ckalloc(sizeof(short) * LCSIZ);
+ lcf = (float *) ckalloc(sizeof(float) * LCSIZ);
len = 1 + (LCSIZ/2);
- fn = PI * 2.0 / (LCSIZ - 1);
- scale = 32767.0/(.5 * LCSIZ);
- for(i=0; i < len; i++)
- lcf[i] = (short) (scale * (.5 + (.4 * cos(fn * ((double)i)))));
+ fn = M_PI * 2.0 / (LCSIZ - 1);
+ scale = 1.0/(.5 * LCSIZ);
+ for (i=0; i < len; i++)
+ lcf[i] = (float) (scale * (.5 + (.4 * cos(fn * ((double)i)))));
}
- do_fir(datain,s->length,dataout,len,lcf,1); /* in downsample.c */
+
+ total_samps = s->length;
+ if (total_samps < (len * 3))
+ total_samps = len * 3;
+
+ datain = (float *) ckalloc(sizeof(float) * total_samps);
+ dataout = (float *) ckalloc(sizeof(float) * total_samps);
+ if (!datain || !dataout) {
+ printf("Can't create a new Signal in highpass()\n");
+ return(NULL);
+ }
+
+ Snack_GetSoundData(s, 0, datain, s->length);
+
+ for (i = s->length; i < total_samps; i++)
+ datain[i] = 0.0;
+
+ if (frame_size > total_samps)
+ frame_size = total_samps;
+
+ for (ndone = 0, init = 1; !(init & 2); ndone += act_size) {
+ act_size = total_samps - len - ndone;
+ if (act_size > frame_size) {
+ out_samps = act_size = frame_size;
+ } else {
+ out_samps = act_size;
+ if (!(init & 1) && ((act_size + len) <= frame_size)) {
+ act_size += len;
+ init = 2;
+ }
+ }
+
+ do_ffir(datain+ndone, total_samps-ndone, dataout+ndone,
+ &out_samps, act_size, len, lcf, 1, 1, init);
+
+ if (init & 1)
+ init = 0;
+ }
+
so = Snack_NewSound(s->samprate, LIN16, s->nchannels);
- if (so == NULL) return(NULL);
- Snack_ResizeSoundStorage(so, s->length);
- for (i = 0; i < s->length; i++) {
- Snack_SetSample(so, 0, i, (float)dataout[i]);
+ if (!so) {
+ printf("Can't create a new Signal in highpass()\n");
+ } else {
+ Snack_ResizeSoundStorage(so, s->length);
+ Snack_PutSoundData(so, 0, dataout, s->length);
+ so->length = s->length;
}
- so->length = s->length;
+
ckfree((void *)dataout);
ckfree((void *)datain);
return(so);
--- a/generic/jkGetF0.c
+++ b/generic/jkGetF0.c
@@ -305,8 +305,10 @@
return(error);
}
-static void get_cand(), peak(), do_ffir();
+static void get_cand(), peak();
static int lc_lin_fir(), downsamp();
+extern void do_ffir();
+extern float *downsample();
/* ----------------------------------------------------------------------- */
void get_fast_cands(fdata, fdsdata, ind, step, size, dec, start, nlags, engref, maxloc, maxval, cp, peaks, locs, ncand, par)
@@ -503,7 +505,7 @@
}
/* ---------------------------------------------------------- */
-static void do_ffir(buf,in_samps,bufo,out_samps,idx, ncoef,fc,invert,skip,init)
+void do_ffir(buf, in_samps, bufo, out_samps, idx, ncoef, fc, invert, skip, init)
/* fc contains 1/2 the coefficients of a symmetric FIR filter with unity
passband gain. This filter is convolved with the signal in buf.
The output is placed in buf2. If(invert), the filter magnitude
@@ -524,10 +526,11 @@
register float *buf1;
buf1 = buf;
- if(ncoef > fsize) {/*allocate memory for full coeff. array and filter memory */ fsize = 0;
+ if(ncoef > fsize) { /* allocate memory for full coeff. array and filter memory */
i = (ncoef+1)*2;
if(!((co = (float *)ckrealloc((void *)co, sizeof(float)*i)) &&
(mem = (float *)ckrealloc((void *)mem, sizeof(float)*i)))) {
+ fsize = 0;
fprintf(stderr,"allocation problems in do_fir()\n");
return;
}
@@ -553,63 +556,63 @@
*dp1 = integral - *dp3;
}
- for(i=ncoef-1, dp1=mem; i-- > 0; ) *dp1++ = 0;
+ /* initialize 1st half with zeros */
+ for(i=ncoef-1, dp1=mem; i-- > 0; ) *dp1++ = 0.0;
}
- else
+ else {
+ /* initialize 1st half with last data from previous run */
for(i=ncoef-1, dp1=mem, sp=state; i-- > 0; ) *dp1++ = *sp++;
-
- i = in_samps;
- resid = 0;
+ }
k = (ncoef << 1) -1; /* inner-product loop limit */
- if(skip <= 1) { /* never used */
-/* *out_samps = i;
- for( ; i-- > 0; ) {
+ if(skip <= 1) {
+ for(i = *out_samps; i-- > 0; ) {
for(j=k, dp1=mem, dp2=co, dp3=mem+1, sum = 0.0; j-- > 0;
*dp1++ = *dp3++ )
sum += *dp2++ * *dp1;
- *--dp1 = *buf++;
- *bufo++ = (sum < 0.0)? sum -0.5 : sum +0.5;
+ *--dp1 = *buf++; /* new data to memory */
+ *bufo++ = (sum < 0.0)? sum -0.5f : sum +0.5f;
}
- if(init & 2) {
+ if(init & 2) {
for(i=ncoef; i-- > 0; ) {
for(j=k, dp1=mem, dp2=co, dp3=mem+1, sum = 0.0; j-- > 0;
*dp1++ = *dp3++ )
sum += *dp2++ * *dp1;
- *--dp1 = 0.0;
- *bufo++ = (sum < 0)? sum -0.5 : sum +0.5;
+ *--dp1 = 0.0; /* pad end with zeros */
+ *bufo++ = (sum < 0.0)? sum -0.5f : sum +0.5f;
}
*out_samps += ncoef;
}
- return;
-*/
- }
+ }
else { /* skip points (e.g. for downsampling) */
/* the buffer end is padded with (ncoef-1) data points */
for( l=0 ; l < *out_samps; l++ ) {
- for(j=k-skip, dp1=mem, dp2=co, dp3=mem+skip, sum=0.0; j-- >0;
+ for(j=k-skip, dp1=mem, dp2=co, dp3=mem+skip, sum=0.0; j-- > 0;
*dp1++ = *dp3++)
sum += *dp2++ * *dp1;
- for(j=skip; j-- >0; *dp1++ = *buf++) /* new data to memory */
+ for(j=skip; j-- > 0; *dp1++ = *buf++) /* new data to memory */
sum += *dp2++ * *dp1;
- *bufo++ = (sum<0.0) ? sum -0.5f : sum +0.5f;
+ *bufo++ = (sum < 0.0) ? sum -0.5f : sum +0.5f;
}
if(init & 2){
resid = in_samps - *out_samps * skip;
for(l=resid/skip; l-- >0; ){
- for(j=k-skip, dp1=mem, dp2=co, dp3=mem+skip, sum=0.0; j-- >0;
+ for(j=k-skip, dp1=mem, dp2=co, dp3=mem+skip, sum=0.0; j-- > 0;
*dp1++ = *dp3++)
sum += *dp2++ * *dp1;
- for(j=skip; j-- >0; *dp1++ = 0.0)
+ for(j=skip; j-- > 0; *dp1++ = 0.0) /* pad end with zeros */
sum += *dp2++ * *dp1;
- *bufo++ = (sum<0.0) ? sum -0.5f : sum +0.5f;
+ *bufo++ = (sum < 0.0) ? sum -0.5f : sum +0.5f;
(*out_samps)++;
}
}
- else
- for(dp3=buf1+idx-ncoef+1, l=ncoef-1, sp=state; l-- >0; ) *sp++ = *dp3++;
+ }
+
+ if(!(init & 2)) { /* unless this is already the end of the signal */
+ /* keep last (ncoef-1) data points for the next initialization */
+ for(dp3=buf1+idx-ncoef+1, l=ncoef-1, sp=state; l-- > 0; ) *sp++ = *dp3++;
}
}
@@ -936,7 +939,7 @@
float **f0p_pt, **vuvp_pt, **rms_speech_pt, **acpkp_pt;
int *vecsize, last_time;
{
- float maxval, engref, *sta, *rms_ratio, *dsdata, *downsample();
+ float maxval, engref, *sta, *rms_ratio, *dsdata;
register float ttemp, ftemp, ft1, ferr, err, errmin;
register int i, j, k, loc1, loc2;
int nframes, maxloc, ncand, ncandp, minloc,
Description: Replace dwnsample() and do_fir() by downsample() and do_ffir() resp.
Utilize the functions downsample() and do_ffir() from jkGetF0.c for the signal
pre-processing. The former functions are from proprietary code and, thus, can't
be used. Notice that do_fir() was a 16-bit integer implementation for an FIR
filter whereas do_ffir() is a floating-point implementation. This might cause a
higher computational effort. On the other side we gain precision. That is why
the volume of the audio signal can stay the same now and we experience much less
effects from quantization.
Author: Thomas Uhle <thomas.u...@mailbox.tu-dresden.de>
Last-Modified: Sat, 18 Feb 2023 20:33:38 +0100
Bug-Debian: https://bugs.debian.org/1029034
--- a/generic/jkFormant.c
+++ b/generic/jkFormant.c
@@ -604,175 +604,6 @@
return(TRUE);
}
-/* Copyright (c) 1987, 1988, 1989 AT&T */
-/* All Rights Reserved */
-
-/* THIS IS UNPUBLISHED PROPRIETARY SOURCE CODE OF AT&T */
-/* The copyright notice above does not evidence any */
-/* actual or intended publication of such source code. */
-
-/* downsample.c */
-/* a quick and dirty downsampler */
-
-#ifndef TRUE
-# define TRUE 1
-# define FALSE 0
-#endif
-
-#define PI 3.1415927
-
-/* ---------------------------------------------------------- */
-int lc_lin_fir(fc,nf,coef)
-/* create the coefficients for a symmetric FIR lowpass filter using the
- window technique with a Hanning window. */
-register double fc;
-double coef[];
-int *nf;
-{
- register int i, n;
- register double twopi, fn, c;
-
- if(((*nf % 2) != 1) || (*nf > 127)) {
- if(*nf <= 126) *nf = *nf + 1;
- else *nf = 127;
- }
- n = (*nf + 1)/2;
-
- /* compute part of the ideal impulse response */
- twopi = PI * 2.0;
- coef[0] = 2.0 * fc;
- c = PI;
- fn = twopi * fc;
- for(i=1;i < n; i++) coef[i] = sin(i * fn)/(c * i);
-
- /* Now apply a Hanning window to the (infinite) impulse response. */
- fn = twopi/((double)(*nf - 1));
- for(i=0;i<n;i++)
- coef[i] *= (.5 + (.5 * cos(fn * ((double)i))));
-
- return(TRUE);
-}
-
-/* ---------------------------------------------------------- */
-
-void do_fir(buf,in_samps,bufo,ncoef,ic,invert)
-/* ic contains 1/2 the coefficients of a symmetric FIR filter with unity
- passband gain. This filter is convolved with the signal in buf.
- The output is placed in buf2. If invert != 0, the filter magnitude
- response will be inverted. */
-short *buf, *bufo, ic[];
-int in_samps, ncoef, invert;
-{
- register short *buft, *bufp, *bufp2, stem;
- short co[256], mem[256];
- register int i, j, k, l, m, sum, integral;
-
- for(i=ncoef-1, bufp=ic+ncoef-1, bufp2=co, buft = co+((ncoef-1)*2),
- integral = 0; i-- > 0; )
- if(!invert) *buft-- = *bufp2++ = *bufp--;
- else {
- integral += (stem = *bufp--);
- *buft-- = *bufp2++ = -stem;
- }
- if(!invert) *buft-- = *bufp2++ = *bufp--; /* point of symmetry */
- else {
- integral *= 2;
- integral += *bufp;
- *buft-- = integral - *bufp;
- }
-/* for(i=(ncoef*2)-2; i >= 0; i--) printf("\n%4d%7d",i,co[i]); */
- for(i=ncoef-1, buft=mem; i-- > 0; ) *buft++ = 0;
- for(i=ncoef; i-- > 0; ) *buft++ = *buf++;
- l = 16384;
- m = 15;
- k = (ncoef << 1) -1;
- for(i=in_samps-ncoef; i-- > 0; ) {
- for(j=k, buft=mem, bufp=co, bufp2=mem+1, sum = 0; j-- > 0;
- *buft++ = *bufp2++ )
- sum += (((*bufp++ * *buft) + l) >> m);
-
- *--buft = *buf++; /* new data to memory */
- *bufo++ = sum;
- }
- for(i=ncoef; i-- > 0; ) { /* pad data end with zeros */
- for(j=k, buft=mem, bufp=co, bufp2=mem+1, sum = 0; j-- > 0;
- *buft++ = *bufp2++ )
- sum += (((*bufp++ * *buft) + l) >> m);
- *--buft = 0;
- *bufo++ = sum;
- }
-}
-
-/* ******************************************************************** */
-
-int get_abs_maximum(d,n)
- register short *d;
- register int n;
-{
- register int i;
- register short amax, t;
-
- if((t = *d++) >= 0) amax = t;
- else amax = -t;
-
- for(i = n-1; i-- > 0; ) {
- if((t = *d++) > amax) amax = t;
- else {
- if(-t > amax) amax = -t;
- }
- }
- return((int)amax);
-}
-
-/* ******************************************************************** */
-
-int dwnsamp(buf,in_samps,buf2,out_samps,insert,decimate,ncoef,ic,smin,smax)
- short *buf, **buf2;
- int in_samps, *out_samps, insert, decimate, ncoef, *smin, *smax;
- short ic[];
-{
- register short *bufp, *bufp2;
- short *buft;
- register int i, j, k, l, m;
- int imax, imin;
-
- if(!(*buf2 = buft = (short*)ckalloc(sizeof(short)*insert*in_samps))) {
- perror("ckalloc() in dwnsamp()");
- return(FALSE);
- }
-
- k = imax = get_abs_maximum(buf,in_samps);
- if (k == 0) k = 1;
- if(insert > 1) k = (32767 * 32767)/k; /* prepare to scale data */
- else k = (16384 * 32767)/k;
- l = 16384;
- m = 15;
-
-
- /* Insert zero samples to boost the sampling frequency and scale the
- signal to maintain maximum precision. */
- for(i=0, bufp=buft, bufp2=buf; i < in_samps; i++) {
- *bufp++ = ((k * (*bufp2++)) + l) >> m ;
- for(j=1; j < insert; j++) *bufp++ = 0;
- }
-
- do_fir(buft,in_samps*insert,buft,ncoef,ic,0);
-
- /* Finally, decimate and return the downsampled signal. */
- *out_samps = j = (in_samps * insert)/decimate;
- k = decimate;
- for(i=0, bufp=buft, imax = imin = *bufp; i < j; bufp += k,i++) {
- *buft++ = *bufp;
- if(imax < *bufp) imax = *bufp;
- else
- if(imin > *bufp) imin = *bufp;
- }
- *smin = imin;
- *smax = imax;
- *buf2 = (short*)ckrealloc((void *) *buf2,sizeof(short) * (*out_samps));
- return(TRUE);
-}
-
/* ---------------------------------------------------------- */
int ratprx(a,k,l,qlim)
@@ -611,13 +611,11 @@
int *l, *k, qlim;
{
double aa, af, q, em, qq = 0, pp = 0, ps, e;
- int ai, ip, i;
+ int ai, ip, result = FALSE;
aa = fabs(a);
ai = (int) aa;
-/* af = fmod(aa,1.0); */
- i = (int) aa;
- af = aa - i;
+ af = aa - ai;
q = 0;
em = 1.0;
while(++q <= qlim) {
@@ -628,97 +626,121 @@
em = e;
pp = ip;
qq = q;
+ result = TRUE;
}
};
*k = (int) ((ai * qq) + pp);
*k = (a > 0)? *k : -(*k);
*l = (int) qq;
- return(TRUE);
+ return(result);
}
/* ----------------------------------------------------------------------- */
+extern float *downsample();
+
Sound *Fdownsample(s,freq2,start,end)
double freq2;
Sound *s;
int start;
int end;
{
- short *bufin, **bufout;
- static double beta = 0.0, b[256];
- double ratio_t, maxi, ratio, beta_new, freq1;
- static int ncoeff = 127, ncoefft = 0, nbits = 15;
- static short ic[256];
- int insert, decimate, out_samps, smin, smax;
+ float *bufin, *bufout, *bufp;
+ int frame_size = 1024, act_size, first_time, last_time;
+ double ratio, freq1;
+ int ncoeff, insert, decimate, total_samps, out_samps, ndone;
Sound *so;
register int i, j;
freq1 = s->samprate;
-
- if((bufout = (short**)ckalloc(sizeof(short*)))) {
- bufin = (short *) ckalloc(sizeof(short) * (end - start + 1));
- for (i = start; i <= end; i++) {
- bufin[i-start] = (short) Snack_GetSample(s, 0, i);
- }
+ ratio = freq2/freq1;
- ratio = freq2/freq1;
- ratprx(ratio,&insert,&decimate,10);
- ratio_t = ((double)insert)/((double)decimate);
+ if (!ratprx(ratio, &insert, &decimate, 10))
+ return(NULL);
- if(ratio_t > .99) return(s);
-
- freq2 = ratio_t * freq1;
- beta_new = (.5 * freq2)/(insert * freq1);
+ if (decimate <= insert)
+ return(NULL);
+
+ freq1 *= (double)insert;
+ freq2 = freq1/((double)decimate);
- if(beta != beta_new){
- beta = beta_new;
- if( !lc_lin_fir(beta,&ncoeff,b)) {
- printf("\nProblems computing interpolation filter\n");
- return(FALSE);
- }
- maxi = (1 << nbits) - 1;
- j = (ncoeff/2) + 1;
- for(ncoefft = 0, i=0; i < j; i++){
- ic[i] = (int) (0.5 + (maxi * b[i]));
- if(ic[i]) ncoefft = i+1;
+ /* filter length used in downsample(): 5ms */
+ ncoeff = ((int)(freq1 * 0.005))/2 + 1;
+
+ total_samps = (end - start + 1) * insert;
+ if (total_samps < (ncoeff * decimate * 3)) /* signal too short */
+ return(NULL);
+
+ if ((bufin = (float *) ckalloc(sizeof(float) * total_samps))) {
+ for (bufp = bufin, i = start; i <= end; i++) {
+ *bufp++ = Snack_GetSample(s, 0, i) * ((float)insert);
+ for(j = 1; j < insert; j++)
+ *bufp++ = 0.0f; /* insert zeros to boost the sampling frequency */
+ }
+
+ if ((frame_size * 2) > total_samps)
+ frame_size = (total_samps + 1)/2;
+
+ frame_size -= frame_size % decimate;
+
+ first_time = 1; /* new filter coefficients need to be computed */
+
+ for (ndone = 0, last_time = 0; !last_time; ndone += act_size) {
+ act_size = total_samps - ncoeff - ndone;
+ if (act_size > frame_size) {
+ act_size = frame_size;
+ out_samps = act_size/decimate;
+ } else {
+ out_samps = act_size/decimate;
+ if (!first_time && ((act_size + ncoeff) <= frame_size)) {
+ act_size += ncoeff;
+ last_time = 1;
+ } else
+ act_size = out_samps * decimate;
}
- } /* endif new coefficients need to be computed */
- if(dwnsamp(bufin,end-start+1,bufout,&out_samps,insert,decimate,ncoefft,ic,
- &smin,&smax)){
- /* so->buff_size = so->file_size = out_samps;*/
- so = Snack_NewSound(0, LIN16, s->nchannels);
- Snack_ResizeSoundStorage(so, out_samps);
- for (i = 0; i < out_samps; i++) {
- Snack_SetSample(so, 0, i, (float)(*bufout)[i]);
+ if ((bufout = downsample(bufin+ndone, total_samps-ndone, act_size, freq1,
+ &out_samps, decimate, first_time, last_time))) {
+ if (first_time) {
+ first_time = 0;
+ so = Snack_NewSound((int)freq2, LIN16, s->nchannels);
+ if (!so) {
+ printf("Can't create a new Signal in downsample()\n");
+ break;
+ }
+ Snack_ResizeSoundStorage(so, total_samps/decimate);
+ so->length = 0;
+ }
+
+ Snack_PutSoundData(so, so->length, bufout, out_samps);
+ so->length += out_samps;
+ } else {
+ printf("Problems in downsample()\n");
+ break;
}
- so->length = out_samps;
- so->samprate = (int)freq2;
- ckfree((void *)*bufout);
- ckfree((void *)bufout);
- ckfree((void *)bufin);
- return(so);
- } else
- printf("Problems in dwnsamp() in downsample()\n");
+ }
+ ckfree((void *)bufin);
} else
printf("Can't create a new Signal in downsample()\n");
- return(NULL);
+ return(so);
}
/* ---------------------------------------------------------- */
-Sound
-*highpass(s)
+extern void do_ffir();
+
+Sound *highpass(s)
Sound *s;
{
- short *datain, *dataout;
- static short *lcf;
+ float *datain, *dataout;
+ static float *lcf;
static int len = 0;
double scale, fn;
register int i;
+ int frame_size = 1024, act_size, total_samps, out_samps, ndone, init;
Sound *so;
/* Header *h, *dup_header();*/
@@ -727,28 +749,62 @@
/* This assumes the sampling frequency is 10kHz and that the FIR
is a Hanning function of (LCSIZ/10)ms duration. */
- datain = (short *) ckalloc(sizeof(short) * s->length);
- dataout = (short *) ckalloc(sizeof(short) * s->length);
- for (i = 0; i < Snack_GetLength(s); i++) {
- datain[i] = (short) Snack_GetSample(s, 0, i);
- }
-
if(!len) { /* need to create a Hanning FIR? */
- lcf = (short*)ckalloc(sizeof(short) * LCSIZ);
+ lcf = (float *) ckalloc(sizeof(float) * LCSIZ);
len = 1 + (LCSIZ/2);
- fn = PI * 2.0 / (LCSIZ - 1);
- scale = 32767.0/(.5 * LCSIZ);
- for(i=0; i < len; i++)
- lcf[i] = (short) (scale * (.5 + (.4 * cos(fn * ((double)i)))));
+ fn = M_PI * 2.0 / (LCSIZ - 1);
+ scale = 1.0/(.5 * LCSIZ);
+ for (i=0; i < len; i++)
+ lcf[i] = (float) (scale * (.5 + (.4 * cos(fn * ((double)i)))));
}
- do_fir(datain,s->length,dataout,len,lcf,1); /* in downsample.c */
+
+ total_samps = s->length;
+ if (total_samps < (len * 3))
+ total_samps = len * 3;
+
+ datain = (float *) ckalloc(sizeof(float) * total_samps);
+ dataout = (float *) ckalloc(sizeof(float) * total_samps);
+ if (!datain || !dataout) {
+ printf("Can't create a new Signal in highpass()\n");
+ return(NULL);
+ }
+
+ Snack_GetSoundData(s, 0, datain, s->length);
+
+ for (i = s->length; i < total_samps; i++)
+ datain[i] = 0.0;
+
+ if (frame_size > total_samps)
+ frame_size = total_samps;
+
+ for (ndone = 0, init = 1; !(init & 2); ndone += act_size) {
+ act_size = total_samps - len - ndone;
+ if (act_size > frame_size) {
+ out_samps = act_size = frame_size;
+ } else {
+ out_samps = act_size;
+ if (!(init & 1) && ((act_size + len) <= frame_size)) {
+ act_size += len;
+ init = 2;
+ }
+ }
+
+ do_ffir(datain+ndone, total_samps-ndone, dataout+ndone,
+ &out_samps, act_size, len, lcf, 1, 1, init);
+
+ if (init & 1)
+ init = 0;
+ }
+
so = Snack_NewSound(s->samprate, LIN16, s->nchannels);
- if (so == NULL) return(NULL);
- Snack_ResizeSoundStorage(so, s->length);
- for (i = 0; i < s->length; i++) {
- Snack_SetSample(so, 0, i, (float)dataout[i]);
+ if (!so) {
+ printf("Can't create a new Signal in highpass()\n");
+ } else {
+ Snack_ResizeSoundStorage(so, s->length);
+ Snack_PutSoundData(so, 0, dataout, s->length);
+ so->length = s->length;
}
- so->length = s->length;
+
ckfree((void *)dataout);
ckfree((void *)datain);
return(so);
--- a/generic/jkGetF0.c
+++ b/generic/jkGetF0.c
@@ -305,8 +305,10 @@
return(error);
}
-static void get_cand(), peak(), do_ffir();
+static void get_cand(), peak();
static int lc_lin_fir(), downsamp();
+extern void do_ffir();
+extern float *downsample();
/* ----------------------------------------------------------------------- */
void get_fast_cands(fdata, fdsdata, ind, step, size, dec, start, nlags, engref, maxloc, maxval, cp, peaks, locs, ncand, par)
@@ -503,7 +505,7 @@
}
/* ---------------------------------------------------------- */
-static void do_ffir(buf,in_samps,bufo,out_samps,idx, ncoef,fc,invert,skip,init)
+void do_ffir(buf, in_samps, bufo, out_samps, idx, ncoef, fc, invert, skip, init)
/* fc contains 1/2 the coefficients of a symmetric FIR filter with unity
passband gain. This filter is convolved with the signal in buf.
The output is placed in buf2. If(invert), the filter magnitude
@@ -524,10 +526,11 @@
register float *buf1;
buf1 = buf;
- if(ncoef > fsize) {/*allocate memory for full coeff. array and filter memory */ fsize = 0;
+ if(ncoef > fsize) { /* allocate memory for full coeff. array and filter memory */
i = (ncoef+1)*2;
if(!((co = (float *)ckrealloc((void *)co, sizeof(float)*i)) &&
(mem = (float *)ckrealloc((void *)mem, sizeof(float)*i)))) {
+ fsize = 0;
fprintf(stderr,"allocation problems in do_fir()\n");
return;
}
@@ -553,63 +556,63 @@
*dp1 = integral - *dp3;
}
- for(i=ncoef-1, dp1=mem; i-- > 0; ) *dp1++ = 0;
+ /* initialize 1st half with zeros */
+ for(i=ncoef-1, dp1=mem; i-- > 0; ) *dp1++ = 0.0;
}
- else
+ else {
+ /* initialize 1st half with last data from previous run */
for(i=ncoef-1, dp1=mem, sp=state; i-- > 0; ) *dp1++ = *sp++;
-
- i = in_samps;
- resid = 0;
+ }
k = (ncoef << 1) -1; /* inner-product loop limit */
- if(skip <= 1) { /* never used */
-/* *out_samps = i;
- for( ; i-- > 0; ) {
+ if(skip <= 1) {
+ for(i = *out_samps; i-- > 0; ) {
for(j=k, dp1=mem, dp2=co, dp3=mem+1, sum = 0.0; j-- > 0;
*dp1++ = *dp3++ )
sum += *dp2++ * *dp1;
- *--dp1 = *buf++;
- *bufo++ = (sum < 0.0)? sum -0.5 : sum +0.5;
+ *--dp1 = *buf++; /* new data to memory */
+ *bufo++ = (sum < 0.0)? sum -0.5f : sum +0.5f;
}
- if(init & 2) {
+ if(init & 2) {
for(i=ncoef; i-- > 0; ) {
for(j=k, dp1=mem, dp2=co, dp3=mem+1, sum = 0.0; j-- > 0;
*dp1++ = *dp3++ )
sum += *dp2++ * *dp1;
- *--dp1 = 0.0;
- *bufo++ = (sum < 0)? sum -0.5 : sum +0.5;
+ *--dp1 = 0.0; /* pad end with zeros */
+ *bufo++ = (sum < 0.0)? sum -0.5f : sum +0.5f;
}
*out_samps += ncoef;
}
- return;
-*/
- }
+ }
else { /* skip points (e.g. for downsampling) */
/* the buffer end is padded with (ncoef-1) data points */
for( l=0 ; l < *out_samps; l++ ) {
- for(j=k-skip, dp1=mem, dp2=co, dp3=mem+skip, sum=0.0; j-- >0;
+ for(j=k-skip, dp1=mem, dp2=co, dp3=mem+skip, sum=0.0; j-- > 0;
*dp1++ = *dp3++)
sum += *dp2++ * *dp1;
- for(j=skip; j-- >0; *dp1++ = *buf++) /* new data to memory */
+ for(j=skip; j-- > 0; *dp1++ = *buf++) /* new data to memory */
sum += *dp2++ * *dp1;
- *bufo++ = (sum<0.0) ? sum -0.5f : sum +0.5f;
+ *bufo++ = (sum < 0.0) ? sum -0.5f : sum +0.5f;
}
if(init & 2){
resid = in_samps - *out_samps * skip;
for(l=resid/skip; l-- >0; ){
- for(j=k-skip, dp1=mem, dp2=co, dp3=mem+skip, sum=0.0; j-- >0;
+ for(j=k-skip, dp1=mem, dp2=co, dp3=mem+skip, sum=0.0; j-- > 0;
*dp1++ = *dp3++)
sum += *dp2++ * *dp1;
- for(j=skip; j-- >0; *dp1++ = 0.0)
+ for(j=skip; j-- > 0; *dp1++ = 0.0) /* pad end with zeros */
sum += *dp2++ * *dp1;
- *bufo++ = (sum<0.0) ? sum -0.5f : sum +0.5f;
+ *bufo++ = (sum < 0.0) ? sum -0.5f : sum +0.5f;
(*out_samps)++;
}
}
- else
- for(dp3=buf1+idx-ncoef+1, l=ncoef-1, sp=state; l-- >0; ) *sp++ = *dp3++;
+ }
+
+ if(!(init & 2)) { /* unless this is already the end of the signal */
+ /* keep last (ncoef-1) data points for the next initialization */
+ for(dp3=buf1+idx-ncoef+1, l=ncoef-1, sp=state; l-- > 0; ) *sp++ = *dp3++;
}
}
@@ -936,7 +939,7 @@
float **f0p_pt, **vuvp_pt, **rms_speech_pt, **acpkp_pt;
int *vecsize, last_time;
{
- float maxval, engref, *sta, *rms_ratio, *dsdata, *downsample();
+ float maxval, engref, *sta, *rms_ratio, *dsdata;
register float ttemp, ftemp, ft1, ferr, err, errmin;
register int i, j, k, loc1, loc2;
int nframes, maxloc, ncand, ncandp, minloc,
