ref: 829553d297fa56cb0f89ae959f87be048b54c9ad
parent: fc4337819992eafd2cc9a9fead11c67d30c5595e
author: cbagwell <cbagwell>
date: Tue Nov 16 17:38:38 EST 1999
Lots more good resample updates.
--- a/Changelog
+++ b/Changelog
@@ -14,6 +14,8 @@
was with 16-bit integer arithmetic. The older version
in sox-12.16 shifted frequencies slightly and was less
accurate.
+ o Extensive rewrite of polyphas.c, should be faster and use
+ less memory now. The sox-12.16 polyphase code had some bugs.
o New effect 'filter' which is a high-quality DSP lowpass/
highpass/bandpass filter using windowed sinc function
methods, like polyphase and resample.
--- a/resdefs.h
+++ b/resdefs.h
@@ -1,39 +1,0 @@
-
-/*
- * FILE: stdefs.h
- * BY: Christopher Lee Fraley
- * DESC: Defines standard stuff for inclusion in C programs.
- * DATE: 6-JUN-88
- * VERS: 1.0 (6-JUN-88, 2:00pm)
- */
-
-
-#define TRUE 1
-#define FALSE 0
-
-/* Some files include both this file and math.h which will most
- * likely already have PI defined.
- */
-#ifndef PI
-#define PI (3.14159265358979232846)
-#endif
-#ifndef PI2
-#define PI2 (6.28318530717958465692)
-#endif
-#define D2R (0.01745329348) /* (2*pi)/360 */
-#define R2D (57.29577951) /* 360/(2*pi) */
-
-#define MAX(x,y) ((x)>(y) ?(x):(y))
-#define MIN(x,y) ((x)<(y) ?(x):(y))
-#define ABS(x) ((x)<0 ?(-(x)):(x))
-#define SGN(x) ((x)<0 ?(-1):((x)==0?(0):(1)))
-
-typedef char BOOL;
-typedef short HWORD;
-typedef unsigned short UHWORD;
-typedef int IWORD;
-#ifndef WORD
-typedef int WORD;
-#endif
-typedef unsigned int UWORD;
-
--- a/sox.1
+++ b/sox.1
@@ -85,7 +85,7 @@
.br
pick
.br
- polyphase [ \fI-w \fR< \fInum\fR / \fIham\fR > ]
+ polyphase [ \fI-w \fR< \fInut\fR / \fIham\fR > ]
[ \fI -width \fR< \fI long \fR / \fIshort \fR / \fI# \fR> ]
[ \fI-cutoff # \fR ]
.br
@@ -139,7 +139,7 @@
.TP 10
\fB-r \fIrate\fR
Give sample rate in Hertz of file. To cause the output file to have
-a different sample rate then the input file, include this option
+a different sample rate than the input file, include this option
with the appropriate rate value along with the output options.
If the input and output files have
different rates then a sample rate change effect must be ran. If a
@@ -176,7 +176,7 @@
\fB-c \fIchannels\fR
The number of sound channels in the data file.
This may be 1, 2, or 4; for mono, stereo, or quad sound data. To cause
-the output file to have a different number of channels then the input
+the output file to have a different number of channels than the input
file, include this option with the approraite value with the output
file options.
If the input and output file have a different number of channels then the
@@ -198,7 +198,7 @@
\fB-p\fR
Run in preview mode and run fast. This will somewhat speed up
sox when the output format has a different number of channels and
-a different rate then the input file. The order that the effects
+a different rate than the input file. The order that the effects
are run in will be arranged for maximum speed and not quality.
.TP 10
\fB-v \fIvolume\fR
@@ -354,7 +354,7 @@
.B ossdsp
OSS /dev/dsp device driver
.br
-This is a psuedo-file type and can be optionally compiled into Sox. Run
+This is a pseudo-file type and can be optionally compiled into Sox. Run
.B sox -h
to see if you have support for this file type. When this driver is used
it allows you to open up the OSS /dev/dsp file and configure it to
@@ -383,7 +383,7 @@
.B sunau
Sun /dev/audio device driver
.br
-This is a psuedo-file type and can be optionally compiled into Sox. Run
+This is a pseudo-file type and can be optionally compiled into Sox. Run
.B sox -h
to see if you have support for this file type. When this driver is used
it allows you to open up a Sun /dev/audio file and configure it to
@@ -488,7 +488,7 @@
Reduce the number of channels by averaging the samples,
or duplicate channels to increase the number of channels.
This effect is automatically used when the number of input
-samples differ then the number of output channels. When reducing
+samples differ from the number of output channels. When reducing
the number of channels it is possible to manually specify the
avg effect and use the \fI-l\fR and \fI-r\fR options to select only
the left or right channel for the output instead of averaging the
@@ -526,7 +526,7 @@
frequency and settling around it.
See \fBfilter\fR for a bandpass effect with steeper shoulders.
.TP
-chorus \fIgain-in gain-out delay decay speed deptch
+chorus \fIgain-in gain-out delay decay speed depth
.TP 10
-s \fR| \fI-t [ \fIdelay decay speed depth -s \fR| \fI-t ... \fR]
Add a chorus to a sound sample. Each quadtuple
@@ -597,7 +597,7 @@
A lowpass filter is obtained by leaving \fIlow\fR unspecified, or 0.
A highpass filter is obtained by leaving \fIhigh\fR unspecified, or 0,
-or greater than or equal to the Nyquist freq.
+or greater than or equal to the Nyquist frequency.
The \fIwindow-len\fR, if unspecified, defaults to 128.
Longer windows give a sharper cutoff, smaller windows a more gradual cutoff.
@@ -657,7 +657,7 @@
Select the left or right channel of a stereo sample,
or one of four channels in a quadrophonic sample.
.TP
-polyphase [ \fI-w \fR< \fInum\fR / \fIham\fR > ]
+polyphase [ \fI-w \fR< \fInut\fR / \fIham\fR > ]
.TP
[ \fI -width \fR< \fI long \fR / \fIshort \fR / \fI# \fR> ]
.TP 10
@@ -664,24 +664,27 @@
[ \fI-cutoff # \fR ]
Translate input sampling rate to output sampling rate via polyphase
interpolation, a DSP algorithm. This method is slow and uses lots
-of RAM, but gives much better results then
+of RAM, but gives much better results than
.B rate.
.br
-w < nut / ham > : select either a Nuttal (~90 dB stopband) or Hamming
-(~43 dB stopband) window.
-.B Warning:
-Nuttall windows require 2x length than Hamming windows. Default is
+(~43 dB stopband) window. Default is
.I nut.
.br
--width long / short / # : specify the width of the filter.
+-width long / short / # : specify the (approximate) width of the filter.
.I long
is 1024 samples;
.I short
is 128 samples. Alternatively, an exact number can be used. Default is
.I long.
+The
+.I short
+option is
+.B not
+recommended, as it produces poor quality results.
.br
-cutoff # : specify the filter cutoff frequency in terms of fraction of
-bandwidth. If upsampling, then this is the fraction of the orignal signal
+bandwidth. If upsampling, then this is the fraction of the original signal
that should go through. If downsampling, this is the fraction of the
signal left after downsampling. Default is 0.95. Remember that
this is a float.
@@ -704,9 +707,9 @@
.B polyphase.
If you are wondering which of
.B Sox's
-rate changing effects to ues, you will want to read a
+rate changing effects to use, you will want to read a
detailed analysis of all of them at http://eakaw2.et.tu-dresden.de/~andreas/resample/resample.html
-[Nov,1999: These tests need to be updated for sox-12.18, which has bugfixes to the
+[Nov,1999: These tests need to be updated for sox-12.17, which has bugfixes to the
resample and polyphase code.]
.TP 10
resample [ \fI-qs\fB | \fI-q\fB | \fI-ql\fB ] [ \fIrolloff\fB [ \fIbeta\fB ] ]\fR
@@ -718,13 +721,31 @@
The \fI-qs\fR, \fI-q\fR, or \fI-ql\fR options specify increased accuracy
at the cost of lower execution speed. By default, linear interpolation
-is used, with a window width about 37 samples at the lower rate.
+is used, with a window width about 45 samples at the lower rate.
This gives an accuracy of about 16 bits, but insufficient stopband rejection
-in the case that you want to have rolloff greater than about 0.85 of
+in the case that you want to have rolloff greater than about 0.80 of
the Nyquist frequency.
The \fI-q*\fR options use quadratic interpolation of filter
-coefficients, resulting in about 22 bits precision.
-\fI-qs\fR, \fI-q\fR, or \fI-ql\fR use window lengths of 37, 75, or 150
+coefficients, resulting in about 24 bits precision.
+.br
+Following is a table of the reasonable defaults which are built-in to sox:
+.br
+ \fBOption Window rolloff beta interpolation\fR
+.br
+ \fB------ ------ ------- ---- -------------\fR
+.br
+ (none) 45 0.80 16 linear
+.br
+ -qs 45 0.80 16 quadratic
+.br
+ -q 75 0.875 16 quadratic
+.br
+ -ql 149 0.94 16 quadratic
+.br
+ \fB------ ------ ------- ---- -------------\fR
+.br
+.br
+\fI-qs\fR, \fI-q\fR, or \fI-ql\fR use window lengths of 45, 75, or 149
samples, respectively, at the lower sample-rate of the two files.
This means progressively sharper stop-band rejection, at proportionally
slower execution times.
@@ -732,7 +753,8 @@
\fIrolloff\fR refers to the cut-off frequency of the
low pass filter and is given in terms of the
Nyquist frequency for the lower sample rate. rolloff therefore should
-be something between 0. and 1., in practice 0.8-0.95. The default is 0.8.
+be something between 0. and 1., in practice 0.8-0.95. The defaults are
+indicated above.
The \fIbeta\fR parameter
determines the type of filter window used. Any value greater than 2.0 is
@@ -739,7 +761,7 @@
the beta for a Kaiser window. Beta <= 2.0 selects a Nuttall window.
If unspecified, the default is a Kaiser window with beta 16.
-In the case of Kaiser window beta > 2.0, lower betas produce a somewhat
+In the case of Kaiser window (beta > 2.0), lower betas produce a somewhat
faster transition from passband to stopband, at the cost of noticeable artifacts.
A beta of 16 is the default, beta less than 10 is not recommended. If you want
a sharper cutoff, don't use low beta's, use a longer sample window.
@@ -749,17 +771,30 @@
just curious about comparing the effects of Nuttall vs. Kaiser windows.
This is the default effect if the two files have different sampling rates.
-Default parameters are Kaiser window of length 37, rolloff 0.80, beta 16,
-linear interpolation. \fI-qs\fR is only slightly slower, but more accurate for
+Default parameters are, as indicated above, Kaiser window of length 45,
+rolloff 0.80, beta 16, linear interpolation.
+
+\fBNOTE: \fI-qs\fR is only slightly slower, but more accurate for
16-bit or higher precision.
+\fBNOTE:\fR In many cases of up-sampling, no interpolation is needed,
+as exact filter coefficients can be computed in a reasonable amount of space.
+To be precise, this is done when
+
+.br
+ input_rate < output_rate
+.br
+ &&
+.br
+ output_rate/gcd(input_rate,output_rate) <= 511
+.br
.TP 10
reverb \fIgain-out delay \fR[ \fIdelay ... \fR]
-Add reverbation to a sound sample. Each delay is given
+Add reverberation to a sound sample. Each delay is given
in milliseconds and its feedback is depending on the
reverb-time in milliseconds. Each delay should be in
the range of half to quarter of reverb-time to get
-a realistic reverbation. Gain-out is the volume of the
+a realistic reverberation. Gain-out is the volume of the
output.
.TP 10
reverse
--- a/src/filter.c
+++ b/src/filter.c
@@ -24,6 +24,7 @@
#define memmove(dest,src,len) bcopy((src),(dest),(len))
#endif
+/* this Float MUST match that in resample.h */
#define Float double/*float*/
#define ISCALE 0x10000
--- a/src/resample.c
+++ b/src/resample.c
@@ -46,12 +46,12 @@
#include "st.h"
/* resample includes */
-#include "resdefs.h"
#include "resampl.h"
#define Float double/*float*/
#define ISCALE 0x10000
+#define NQMAX 511
#define BUFFSIZE 8192 /*16384*/ /* Total buffer size */
#define L64 long long
@@ -65,13 +65,18 @@
LONG Nwing;
LONG Nq;
Float *Imp; /* impulse [Nwing+1] Filter coefficients */
+
double Time; /* Current time/pos in input sample */
LONG dhb;
+
+ LONG a,b; /* gcd-reduced input,output rates */
+ LONG t; /* Current time/pos for exact-coeff's method */
+
LONG Xh; /* number of past/future samples needed by filter */
- LONG Xoff; /* Xh plus some room for creep */
+ LONG Xoff; /* Xh plus some room for creep */
LONG Xread; /* X[Xread] is start-position to enter new samples */
LONG Xp; /* X[Xp] is position to start filter application */
- LONG Xsize,Ysize; /* size (Floats) of X[],Y[] */
+ LONG Xsize,Ysize; /* size (Floats) of X[],Y[] */
Float *X, *Y; /* I/O buffers */
} *resample_t;
@@ -88,13 +93,9 @@
LONG Num));
static LONG SrcUD(P2(resample_t r, LONG Nx));
+static LONG SrcEX(P2(resample_t r, LONG Nx));
-#if 0
-static u_int32_t iprodC;
-static u_int32_t iprodM;
-#endif
-
/*
* Process options
*/
@@ -106,10 +107,10 @@
resample_t r = (resample_t) effp->priv;
/* These defaults are conservative with respect to aliasing. */
- r->rolloff = 0.8;
+ r->rolloff = 0.80;
r->beta = 16; /* anything <=2 means Nutall window */
r->quadr = 0;
- r->Nmult = 39;
+ r->Nmult = 45;
/* This used to fail, but with sox-12.15 it works. AW */
if ((n >= 1)) {@@ -118,13 +119,13 @@
n--; argv++;
}
else if (!strcmp(argv[0], "-q")) {- r->rolloff = 0.9;
+ r->rolloff = 0.875;
r->quadr = 1;
r->Nmult = 75;
n--; argv++;
}
else if (!strcmp(argv[0], "-ql")) {- r->rolloff = 0.9;
+ r->rolloff = 0.94;
r->quadr = 1;
r->Nmult = 149;
n--; argv++;
@@ -154,13 +155,23 @@
eff_t effp;
{resample_t r = (resample_t) effp->priv;
- LONG Xoff;
+ LONG Xoff, gcdrate;
int i;
+ extern long st_gcd(P2(long a,long b));
r->Factor = (double)effp->outinfo.rate / (double)effp->ininfo.rate;
- r->Nq = Nc; /* for now */
+ gcdrate = st_gcd((long)effp->ininfo.rate, (long)effp->outinfo.rate);
+ r->a = effp->ininfo.rate / gcdrate;
+ r->b = effp->outinfo.rate / gcdrate;
+ if (r->a <= r->b && r->b <= NQMAX) {+ r->quadr = -1; /* exact coeff's */
+ r->Nq = r->b; /* MAX(r->a,r->b); */
+ } else {+ r->Nq = Nc; /* for now */
+ }
+
/* Check for illegal constants */
# if 0
if (Lp >= 16) fail("Error: Lp>=16");@@ -178,12 +189,17 @@
if (i <= 0)
fail("resample: Unable to make filter\n");- report("Nmult: %ld, Nwing: %ld, Nq: %ld\n",r->Nmult,r->Nwing,r->Nq);+ /*report("Nmult: %ld, Nwing: %ld, Nq: %ld\n",r->Nmult,r->Nwing,r->Nq);*/- r->dhb = Np; /* Fixed-point Filter sampling-time-increment */
- if (r->Factor<1.0) r->dhb = r->Factor*Np + 0.5;
- r->Xh = (r->Nwing<<La)/r->dhb;
- /* (Xh * dhb)>>La is max index into Imp[] */
+ if (r->quadr < 0) { /* exact coeff's method */+ r->Xh = r->Nwing/r->b;
+ report("resample: rate ratio %ld:%ld, coeff interpolation not needed\n", r->a, r->b);+ } else {+ r->dhb = Np; /* Fixed-point Filter sampling-time-increment */
+ if (r->Factor<1.0) r->dhb = r->Factor*Np + 0.5;
+ r->Xh = (r->Nwing<<La)/r->dhb;
+ /* (Xh * dhb)>>La is max index into Imp[] */
+ }
/* reach of LP filter wings + some creeping room */
Xoff = r->Xh + 10;
@@ -195,7 +211,9 @@
r->Xread = Xoff;
/* Current-time pointer for converter */
r->Time = Xoff;
-
+ if (r->quadr < 0) { /* exact coeff's method */+ r->t = Xoff*r->Nq;
+ }
i = BUFFSIZE - 2*Xoff;
if (i < r->Factor + 1.0/r->Factor) /* Check input buffer size */
fail("Factor is too small or large for BUFFSIZE");@@ -202,7 +220,7 @@
r->Xsize = 2*Xoff + i/(1.0+r->Factor);
r->Ysize = BUFFSIZE - r->Xsize;
- report("Xsize %d, Ysize %d, Xoff %d",r->Xsize,r->Ysize,r->Xoff);+ /* report("Xsize %d, Ysize %d, Xoff %d",r->Xsize,r->Ysize,r->Xoff); */r->X = (Float *) malloc(sizeof(Float) * (BUFFSIZE));
r->Y = r->X + r->Xsize;
@@ -226,11 +244,11 @@
LONG i, last, Nout, Nx, Nproc;
/* constrain amount we actually process */
- /*fprintf(stderr,"Xp %d, Xread %d, isamp %d, ",r->Xp, r->Xread,*isamp); */
+ /*fprintf(stderr,"Xp %d, Xread %d, isamp %d, ",r->Xp, r->Xread,*isamp);*/
Nproc = r->Xsize - r->Xp;
- i = MIN(r->Ysize, *osamp);
+ i = (r->Ysize < *osamp)? r->Ysize : *osamp;
if (Nproc * r->Factor >= i)
Nproc = i / r->Factor;
@@ -258,21 +276,38 @@
*osamp = 0;
return;
}
- Nout = SrcUD(r, Nproc);
- /*fprintf(stderr,"Nproc %d --> %d\n",Nproc,Nout);*/
- /* Move converter Nproc samples back in time */
- r->Time -= Nproc;
- /* Advance by number of samples processed */
- r->Xp += Nproc;
- /* Calc time accumulation in Time */
- {- LONG creep = r->Time - r->Xoff;
- if (creep)
- {- r->Time -= creep; /* Remove time accumulation */
- r->Xp += creep; /* and add it to read pointer */
- /* fprintf(stderr,"Nproc %ld, creep %ld\n",Nproc,creep); */
- }
+ if (r->quadr < 0) { /* exact coeff's method */+ LONG creep;
+ Nout = SrcEX(r, Nproc);
+ /*fprintf(stderr,"Nproc %d --> %d\n",Nproc,Nout);*/
+ /* Move converter Nproc samples back in time */
+ r->t -= Nproc * r->b;
+ /* Advance by number of samples processed */
+ r->Xp += Nproc;
+ /* Calc time accumulation in Time */
+ creep = r->t/r->b - r->Xoff;
+ if (creep)
+ {+ r->t -= creep * r->b; /* Remove time accumulation */
+ r->Xp += creep; /* and add it to read pointer */
+ /*fprintf(stderr,"Nproc %ld, creep %ld\n",Nproc,creep);*/
+ }
+ } else { /* approx coeff's method */+ LONG creep;
+ Nout = SrcUD(r, Nproc);
+ /*fprintf(stderr,"Nproc %d --> %d\n",Nproc,Nout);*/
+ /* Move converter Nproc samples back in time */
+ r->Time -= Nproc;
+ /* Advance by number of samples processed */
+ r->Xp += Nproc;
+ /* Calc time accumulation in Time */
+ creep = r->Time - r->Xoff;
+ if (creep)
+ {+ r->Time -= creep; /* Remove time accumulation */
+ r->Xp += creep; /* and add it to read pointer */
+ /* fprintf(stderr,"Nproc %ld, creep %ld\n",Nproc,creep); */
+ }
}
{@@ -305,29 +340,29 @@
LONG *osamp;
{resample_t r = (resample_t) effp->priv;
- LONG i, Nout, Nx;
-
- /*fprintf(stderr,"Xp %d, Xread %d <--- DRAIN\n",r->Xp, r->Xread);*/
- if (r->Xsize - r->Xread < r->Xoff)
- fail("resample_drain: Problem!\n");+ LONG isamp_res, *Obuf, osamp_res;
- /* fill out end with Xoff zeros */
- for(i = 0; i < r->Xoff; i++)
- r->X[i + r->Xread] = 0;
+ /* fprintf(stderr,"Xoff %d, Xt %d <--- DRAIN\n",r->Xoff, r->Xt); */
- Nx = r->Xread - r->Xp;
-
- if (Nx * r->Factor >= *osamp)
- fail("resample_drain: Overran output buffer!\n");-
- /* Resample stuff in input buffer */
- Nout = SrcUD(r, Nx);
- /*fprintf(stderr,"Nproc %d --> %d\n",Nx,Nout);*/
-
- for(i = 0; i < Nout; i++)
- *obuf++ = r->Y[i] * ISCALE;
-
- *osamp = Nout;
+ /* stuff end with Xoff zeros */
+ isamp_res = r->Xoff;
+ osamp_res = *osamp;
+ Obuf = obuf;
+ while (isamp_res>0 && osamp_res>0) {+ LONG Isamp, Osamp;
+ Isamp = isamp_res;
+ Osamp = osamp_res;
+ resample_flow(effp, NULL, Obuf, &Isamp, &Osamp);
+ /* fprintf(stderr,"DRAIN isamp,osamp (%d,%d) -> (%d,%d)\n",
+ isamp_res,osamp_res,Isamp,Osamp); */
+ Obuf += Osamp;
+ osamp_res -= Osamp;
+ isamp_res -= Isamp;
+ }
+ *osamp -= osamp_res;
+ /* fprintf(stderr,"DRAIN osamp %d\n", *osamp); */
+ if (isamp_res)
+ warn("drain overran obuf by %d\n", isamp_res); fflush(stderr);}
/*
@@ -342,7 +377,6 @@
free(r->Imp - 1);
free(r->X);
/* free(r->Y); Y is in same block starting at X */
- /*report("iC %d, iM %d, ratio %d", iprodC, iprodM, iprodM/iprodC);*/}
/* over 90% of CPU time spent in this iprodUD() function */
@@ -356,15 +390,15 @@
double v;
LONG Ho;
- Ho = T0 * dhb;
- Ho += (ct-1)*dhb; /* so Float sum starts with smallest coef's */
- Xp += (ct-1)*Inc;
- v = 0;
+ Ho = T0 * dhb;
+ Ho += (ct-1)*dhb; /* so Float sum starts with smallest coef's */
+ Xp += (ct-1)*Inc;
+ v = 0;
do {Float coef;
LONG Hoh;
Hoh = Ho>>La;
- coef = Imp[Hoh];
+ coef = Imp[Hoh];
{Float dm,dp,t;
dm = coef - Imp[Hoh-1];
@@ -391,9 +425,9 @@
LONG Ho;
Ho = T0 * dhb;
- Ho += (ct-1)*dhb; /* so Float sum starts with smallest coef's */
+ Ho += (ct-1)*dhb; /* so Float sum starts with smallest coef's */
Xp += (ct-1)*Inc;
- v = 0;
+ v = 0;
do {Float coef;
LONG Hoh;
@@ -428,9 +462,11 @@
dt = 1.0/Factor; /* Output sampling period */
/*fprintf(stderr,"Factor %f, dt %f, ",Factor,dt); */
/*fprintf(stderr,"Time %f, ",r->Time);*/
- /* (Xh * dhb)>>La is max index into Imp[] */
- /*fprintf(stderr,"ct=%d\n",ct);*/
+ /* (Xh * dhb)>>La is max index into Imp[] */
+ /*fprintf(stderr,"ct=%d\n",ct);*/
/*fprintf(stderr,"ct=%.2f %d\n",(double)r->Nwing*Na/r->dhb, r->Xh);*/
+ /*fprintf(stderr,"ct=%ld, T=%.6f, dhb=%6f, dt=%.6f\n",
+ r->Xh, time-floor(time),(double)r->dhb/Na,dt);*/
Ystart = Y = r->Y;
n = (int)ceil((double)Nx/dt);
while(n--)
@@ -439,7 +475,7 @@
double v;
double T;
T = time-floor(time); /* fractional part of Time */
- Xp = r->X + (LONG)time; /* Ptr to current input sample */
+ Xp = r->X + (LONG)time; /* Ptr to current input sample */
/* Past inner product: */
v = (*prodUD)(r->Imp, Xp, -1, T, r->dhb, r->Xh); /* needs Np*Nmult in 31 bits */
@@ -455,6 +491,62 @@
return (Y - Ystart); /* Return the number of output samples */
}
+/* exact coeff's */
+static double prodEX(Imp, Xp, Inc, T0, dhb, ct)
+const Float Imp[], *Xp;
+LONG Inc, T0, dhb, ct;
+{+ double v;
+ const Float *Cp;
+
+ Cp = Imp + (ct-1)*dhb + T0; /* so Float sum starts with smallest coef's */
+ Xp += (ct-1)*Inc;
+ v = 0;
+ do {+ v += *Cp * *Xp; /* sum coeff * input sample */
+ Cp -= dhb; /* IR step */
+ Xp -= Inc; /* Input signal step. */
+ } while(--ct);
+ return v;
+}
+
+static LONG SrcEX(r, Nx)
+resample_t r;
+LONG Nx;
+{+ Float *Ystart, *Y;
+ double Factor;
+ LONG a,b;
+ LONG time;
+ int n;
+
+ Factor = r->Factor;
+ time = r->t;
+ a = r->a;
+ b = r->b;
+ Ystart = Y = r->Y;
+ n = (Nx*b + (a-1))/a;
+ while(n--)
+ {+ Float *Xp;
+ double v;
+ LONG T;
+ T = time % b; /* fractional part of Time */
+ Xp = r->X + (time/b); /* Ptr to current input sample */
+
+ /* Past inner product: */
+ v = prodEX(r->Imp, Xp, -1, T, b, r->Xh);
+ /* Future inner product: */
+ v += prodEX(r->Imp, Xp+1, 1, b-T, b, r->Xh);
+
+ if (Factor < 1) v *= Factor;
+ *Y++ = v; /* Deposit output */
+ time += a; /* Move to next sample by time increment */
+ }
+ r->t = time;
+ return (Y - Ystart); /* Return the number of output samples */
+}
+
int makeFilter(Imp, Nwing, Froll, Beta, Num)
Float Imp[];
LONG Nwing, Num;
@@ -486,7 +578,7 @@
for (i=Dh; i<Mwing; i+=Dh)
DCgain += ImpR[i];
DCgain = 2*DCgain + ImpR[0]; /* DC gain of real coefficients */
- report("DCgain err=%.12f",DCgain-1.0);+ /*report("DCgain err=%.12f",DCgain-1.0);*/ for (i=0; i<Mwing; i++) {Imp[i] = ImpR[i]/DCgain;
@@ -496,21 +588,6 @@
/* Imp[Nwing] and Imp[-1] needed for quadratic interpolation */
Imp[-1] = Imp[1];
-# if 0
- {- double v = 0;
- int s=-1;
- for (i=Num; i<Mwing; i+=Num, s=-s)
- v += s*Imp[i];
- report("NYQgain %.12f vs %.12f",Imp[0], -2*v);- v = -2*v - Imp[0];
- v *= (double)Nc/(double)(2*Mwing+1);
- for (i=0; i<Mwing; i++)
- Imp[i] += v*cos(M_PI*i/Nc);
-
- }
-# endif
-
return(Mwing);
}
--- a/test/ding.c
+++ b/test/ding.c
@@ -42,6 +42,15 @@
# define MINSAMPL -0x7fffffff
#endif
+struct _Env {+ int r; /* rise */
+ int m; /* middle */
+ int f; /* fall */
+ int e; /* end */
+ FLOAT v; /* volume */
+ FLOAT d; /* decay */
+} Env = {0,0,0,0,0.5,1.0};+
static void Usage(void)__attribute__((noreturn));
static void Usage(void)
@@ -101,22 +110,13 @@
int poflo,moflo;
FLOAT Vol0=1;
FLOAT Freq=0; /* of nyquist */
- int Offset=0;
int Pad=0;
- struct _Env {- int r; /* rise */
- int c; /* center */
- int f; /* fall */
- FLOAT v; /* volume */
- FLOAT d; /* decay */
- } Env = {0,0,0,0.5,1.0};-
double x=1,y=0;
double thx=1,thy=0;
static inline void a_init(double frq)
{- if (frq != 0) {+ if (0<frq && frq<1) {x = 1; y = 0;
} else {x = 0; y = 1;
@@ -125,15 +125,6 @@
thy = sin(frq*M_PI);
}
- static inline const double a(int k, int L)
- {- double a, u;
- u = (double)k/L; /* in [0,1] */
- a = 0.5*(1 + cos(M_PI * u));
- //printf("a(%d,%d) = %.5f, l=%.1f, u=%.6f\n",k,L,a,l,u); - return a;
- }
-
static inline void a_post(int k)
{double x1;
@@ -148,6 +139,27 @@
}
}
+ /* rise/fall function, monotonic [0,1] -> [1-0] */
+ static inline const double a(double s)
+ {+ return 0.5*(1 + cos(M_PI * s));
+ }
+
+ static double env(struct _Env *Env, int k)
+ {+ double u = 0;
+ //if (k >= Env->r && k < Env->e) {+ if (k<Env->m) {+ u = Env->v * a((double)(Env->m-k)/(Env->m-Env->r));
+ }else if (k<Env->f) {+ u = Env->v;
+ }else{+ u = Env->v * a((double)(k-Env->f)/(Env->e-Env->f));
+ }
+ //}
+ return u;
+ }
+
/* Parse the options */
while ((optc = getopt(argct, argv, "d:o:e:t:p:f:v:h")) != -1) {char *p;
@@ -168,22 +180,22 @@
break;
case 'e':
p = optarg;
- Env.c = strtol(p,&p,10);
+ Env.f = strtol(p,&p,10);
if (*p++ == ':') {- Env.r = Env.f = Env.c;
- Env.c = strtol(p,&p,10);
+ Env.m = Env.e = Env.f;
+ Env.f = strtol(p,&p,10);
}
if (*p++ == ':') {- Env.f = strtol(p,&p,10);
+ Env.e = strtol(p,&p,10);
}
- if (*p || Env.c<=0 || Env.r<0 || Env.f<0) {+ if (*p || Env.f<=0 || Env.m<0 || Env.e<0) {fprintf(stderr,"option -%c not valid (%s)\n", optc, optarg);
Usage();
}
break;
case 'o':
- Offset = strtol(optarg,&p,10);
- if (p==optarg || *p || Offset<0) {+ Env.r = strtol(optarg,&p,10);
+ if (p==optarg || *p || Env.r<0) {fprintf(stderr,"option -%c expects int value (%s)\n", optc, optarg);
Usage();
}
@@ -212,7 +224,10 @@
if (Freq==0.0) Env.v *= sqrt(0.5);
//fprintf(stderr,"Vol0 %8.3f\n", Vol0);
- len = Offset+Env.r+Env.c+Env.f+Pad;
+ Env.m += Env.r;
+ Env.f += Env.m;
+ Env.e += Env.f;
+ len = Env.e+Pad;
fnam1=NULL; fd1=-1;
//fprintf(stderr,"optind=%d argct=%d\n",optind,argct);
if (optind <= argct-2) {@@ -249,7 +264,6 @@
ct = len - st; if (ct>BSIZ) ct=BSIZ;
ReadN(fd1,ibuff,ct);
for (ibp=ibuff; ibp<ibuff+ct; ibp++,st++) {- int k;
double v = *ibp;
if (max<*ibp) max=*ibp;
@@ -256,15 +270,8 @@
else if (min>*ibp) min=*ibp;
v *= Vol0;
- k = st-Offset;
- if (k>=0 && k<Env.r) {- v += y*a(Env.r-k,Env.r)*Env.v;
- a_post(st);
- }else if (k-=Env.r, k>=0 && k<Env.c) {- v += y*Env.v;
- a_post(st);
- }else if (k-=Env.c, k>=0 && k<Env.f) {- v += y*a(k,Env.f)*Env.v;
+ if (st>=Env.r && st<Env.e) {+ v += y*env(&Env,st);
a_post(st);
}
--- a/test/ltest.pl
+++ b/test/ltest.pl
@@ -30,10 +30,15 @@
#my $effect='polyphase -cutoff 0.90';
#my $effect='filter 400-2000';
#my $effect='filter 400-2000 1024';
+
my ($rate0,$rate1)=(8000,22050); # sample rates
my $p=400; # silence before/after tonepulse
my $env="-e4000:16000:4000"; # attack, duration, drop
+#my ($rate0,$rate1)=(22050,8000); # sample rates
+#my $p=1102; # silence before/after tonepulse
+#my $env="-e11025:44100:11025"; # attack, duration, drop
+
# parse commandline arguments
my $updown = 0; # set to 1 for up/down rate-conversion test
my ($ding,$model,$t,$rms,$lim)=("./ding","./model","sw",11585.2, 0.5);@@ -72,8 +77,11 @@
mkdir("d",0775);my $f;
my %q;
-for ($f=0.00; $f<1.0; $f+=0.01) {+my $nyq = ($rate0<=$rate1)? 1.0:($rate1/$rate0);
+for ($f=0.00; $f<1.001; $f+=0.01) {my @mod;
+
+ #if ($f>0.995) { $f=0.999; } my $s=sprintf("%4.2f",$f);#print "$ding -v0.5 -o$p $env -p$p -f$s -d1.0 d/i$s.$t\n";
qx{$ding -v0.5 -o$p $env -p$p -f$s -d1.0 d/i$s.$t &>/dev/null};--- a/test/model.c
+++ b/test/model.c
@@ -163,8 +163,8 @@
}
c = eCenter(ibuff,len);
- Len1 = Env.c*Factor*0.6;
- c += Env.c*Factor*0.15;
+ Len1 = Env.c*Factor*0.6; /* 60% of original const-amplitude area */
+ c += Env.c*Factor*0.15; /* beginning after 30% */
{double a;
int b;
@@ -185,7 +185,6 @@
s = ip[n]; /* sigval at n */
- //if (fabs(n-Len1/2)<Len1/3) {f = 1; //a(n+del,Len1);
u = f*x; v = f*y;
sx += s*u;
@@ -194,7 +193,7 @@
h11 += u*u;
h12 += u*v;
h22 += v*v;
- //}
+
a_post(n);
}
//fprintf(stderr,"h12 %.10f\n", (double)h12/(sqrt(h11)*sqrt(h22)));
@@ -221,7 +220,7 @@
char *fnam1;
/* Parse the options */
- while ((optc = getopt(argct, argv, "d:o:e:p:f:h")) != -1) {+ while ((optc = getopt(argct, argv, "d:e:f:h")) != -1) {char *p;
switch(optc) {case 'd':