![(please configure the [header_logo] section in trac.ini)](/repos/chrome/site/pj.jpg)
- Timestamp:
- Apr 7, 2007 2:53:15 PM (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
pjproject/branches/split-3rd-party/pjmedia/src/pjmedia/resample.c
r974 r1170 18 18 */ 19 19 20 /*21 * Based on:22 * resample-1.8.tar.gz from the23 * Digital Audio Resampling Home Page located at24 * http://www-ccrma.stanford.edu/~jos/resample/.25 *26 * SOFTWARE FOR SAMPLING-RATE CONVERSION AND FIR DIGITAL FILTER DESIGN27 *28 * Snippet from the resample.1 man page:29 *30 * HISTORY31 *32 * The first version of this software was written by Julius O. Smith III33 * <jos@ccrma.stanford.edu> at CCRMA <http://www-ccrma.stanford.edu> in34 * 1981. It was called SRCONV and was written in SAIL for PDP-1035 * compatible machines. The algorithm was first published in36 *37 * Smith, Julius O. and Phil Gossett. ``A Flexible Sampling-Rate38 * Conversion Method,'' Proceedings (2): 19.4.1-19.4.4, IEEE Conference39 * on Acoustics, Speech, and Signal Processing, San Diego, March 1984.40 *41 * An expanded tutorial based on this paper is available at the Digital42 * Audio Resampling Home Page given above.43 *44 * Circa 1988, the SRCONV program was translated from SAIL to C by45 * Christopher Lee Fraley working with Roger Dannenberg at CMU.46 *47 * Since then, the C version has been maintained by jos.48 *49 * Sndlib support was added 6/99 by John Gibson <jgg9c@virginia.edu>.50 *51 * The resample program is free software distributed in accordance52 * with the Lesser GNU Public License (LGPL). There is NO warranty; not53 * even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.54 */55 56 /* PJMEDIA modification:57 * - remove resample(), just use SrcUp, SrcUD, and SrcLinear directly.58 * - move FilterUp() and FilterUD() from filterkit.c59 * - move stddefs.h and resample.h to this file.60 * - const correctness.61 */62 20 #include <pjmedia/resample.h> 63 21 #include <pjmedia/errno.h> … … 66 24 #include <pj/pool.h> 67 25 26 #include "../../third_party/build/resample/resamplesubs.h" 68 27 69 28 #define THIS_FILE "resample.c" 70 29 71 72 /*73 * Taken from stddefs.h74 */75 #ifndef PI76 #define PI (3.14159265358979232846)77 #endif78 79 #ifndef PI280 #define PI2 (6.28318530717958465692)81 #endif82 83 #define D2R (0.01745329348) /* (2*pi)/360 */84 #define R2D (57.29577951) /* 360/(2*pi) */85 86 #ifndef MAX87 #define MAX(x,y) ((x)>(y) ?(x):(y))88 #endif89 #ifndef MIN90 #define MIN(x,y) ((x)<(y) ?(x):(y))91 #endif92 93 #ifndef ABS94 #define ABS(x) ((x)<0 ?(-(x)):(x))95 #endif96 97 #ifndef SGN98 #define SGN(x) ((x)<0 ?(-1):((x)==0?(0):(1)))99 #endif100 101 typedef char RES_BOOL;102 typedef short RES_HWORD;103 typedef int RES_WORD;104 typedef unsigned short RES_UHWORD;105 typedef unsigned int RES_UWORD;106 107 #define MAX_HWORD (32767)108 #define MIN_HWORD (-32768)109 110 #ifdef DEBUG111 #define INLINE112 #else113 #define INLINE inline114 #endif115 116 /*117 * Taken from resample.h118 *119 * The configuration constants below govern120 * the number of bits in the input sample and filter coefficients, the121 * number of bits to the right of the binary-point for fixed-point math, etc.122 *123 */124 125 /* Conversion constants */126 #define Nhc 8127 #define Na 7128 #define Np (Nhc+Na)129 #define Npc (1<<Nhc)130 #define Amask ((1<<Na)-1)131 #define Pmask ((1<<Np)-1)132 #define Nh 16133 #define Nb 16134 #define Nhxn 14135 #define Nhg (Nh-Nhxn)136 #define NLpScl 13137 138 /* Description of constants:139 *140 * Npc - is the number of look-up values available for the lowpass filter141 * between the beginning of its impulse response and the "cutoff time"142 * of the filter. The cutoff time is defined as the reciprocal of the143 * lowpass-filter cut off frequence in Hz. For example, if the144 * lowpass filter were a sinc function, Npc would be the index of the145 * impulse-response lookup-table corresponding to the first zero-146 * crossing of the sinc function. (The inverse first zero-crossing147 * time of a sinc function equals its nominal cutoff frequency in Hz.)148 * Npc must be a power of 2 due to the details of the current149 * implementation. The default value of 512 is sufficiently high that150 * using linear interpolation to fill in between the table entries151 * gives approximately 16-bit accuracy in filter coefficients.152 *153 * Nhc - is log base 2 of Npc.154 *155 * Na - is the number of bits devoted to linear interpolation of the156 * filter coefficients.157 *158 * Np - is Na + Nhc, the number of bits to the right of the binary point159 * in the integer "time" variable. To the left of the point, it indexes160 * the input array (X), and to the right, it is interpreted as a number161 * between 0 and 1 sample of the input X. Np must be less than 16 in162 * this implementation.163 *164 * Nh - is the number of bits in the filter coefficients. The sum of Nh and165 * the number of bits in the input data (typically 16) cannot exceed 32.166 * Thus Nh should be 16. The largest filter coefficient should nearly167 * fill 16 bits (32767).168 *169 * Nb - is the number of bits in the input data. The sum of Nb and Nh cannot170 * exceed 32.171 *172 * Nhxn - is the number of bits to right shift after multiplying each input173 * sample times a filter coefficient. It can be as great as Nh and as174 * small as 0. Nhxn = Nh-2 gives 2 guard bits in the multiply-add175 * accumulation. If Nhxn=0, the accumulation will soon overflow 32 bits.176 *177 * Nhg - is the number of guard bits in mpy-add accumulation (equal to Nh-Nhxn)178 *179 * NLpScl - is the number of bits allocated to the unity-gain normalization180 * factor. The output of the lowpass filter is multiplied by LpScl and181 * then right-shifted NLpScl bits. To avoid overflow, we must have182 * Nb+Nhg+NLpScl < 32.183 */184 185 186 #ifdef _MSC_VER187 # pragma warning(push, 3)188 //# pragma warning(disable: 4245) // Conversion from uint to ushort189 # pragma warning(disable: 4244) // Conversion from double to uint190 # pragma warning(disable: 4146) // unary minus operator applied to unsigned type, result still unsigned191 # pragma warning(disable: 4761) // integral size mismatch in argument; conversion supplied192 #endif193 194 #if defined(PJMEDIA_HAS_SMALL_FILTER) && PJMEDIA_HAS_SMALL_FILTER!=0195 # include "smallfilter.h"196 #else197 # define SMALL_FILTER_NMULT 0198 # define SMALL_FILTER_SCALE 0199 # define SMALL_FILTER_NWING 0200 # define SMALL_FILTER_IMP NULL201 # define SMALL_FILTER_IMPD NULL202 #endif203 204 #if defined(PJMEDIA_HAS_LARGE_FILTER) && PJMEDIA_HAS_LARGE_FILTER!=0205 # include "largefilter.h"206 #else207 # define LARGE_FILTER_NMULT 0208 # define LARGE_FILTER_SCALE 0209 # define LARGE_FILTER_NWING 0210 # define LARGE_FILTER_IMP NULL211 # define LARGE_FILTER_IMPD NULL212 #endif213 214 215 #undef INLINE216 #define INLINE217 #define HAVE_FILTER 0218 219 #ifndef NULL220 # define NULL 0221 #endif222 223 224 static INLINE RES_HWORD WordToHword(RES_WORD v, int scl)225 {226 RES_HWORD out;227 RES_WORD llsb = (1<<(scl-1));228 v += llsb; /* round */229 v >>= scl;230 if (v>MAX_HWORD) {231 v = MAX_HWORD;232 } else if (v < MIN_HWORD) {233 v = MIN_HWORD;234 }235 out = (RES_HWORD) v;236 return out;237 }238 239 /* Sampling rate conversion using linear interpolation for maximum speed.240 */241 static int242 SrcLinear(const RES_HWORD X[], RES_HWORD Y[], double pFactor, RES_UHWORD nx)243 {244 RES_HWORD iconst;245 RES_UWORD time = 0;246 const RES_HWORD *xp;247 RES_HWORD *Ystart, *Yend;248 RES_WORD v,x1,x2;249 250 double dt; /* Step through input signal */251 RES_UWORD dtb; /* Fixed-point version of Dt */252 RES_UWORD endTime; /* When time reaches EndTime, return to user */253 254 dt = 1.0/pFactor; /* Output sampling period */255 dtb = dt*(1<<Np) + 0.5; /* Fixed-point representation */256 257 Ystart = Y;258 Yend = Ystart + (unsigned)(nx * pFactor);259 endTime = time + (1<<Np)*(RES_WORD)nx;260 while (time < endTime)261 {262 iconst = (time) & Pmask;263 xp = &X[(time)>>Np]; /* Ptr to current input sample */264 x1 = *xp++;265 x2 = *xp;266 x1 *= ((1<<Np)-iconst);267 x2 *= iconst;268 v = x1 + x2;269 *Y++ = WordToHword(v,Np); /* Deposit output */270 time += dtb; /* Move to next sample by time increment */271 }272 return (Y - Ystart); /* Return number of output samples */273 }274 275 static RES_WORD FilterUp(const RES_HWORD Imp[], const RES_HWORD ImpD[],276 RES_UHWORD Nwing, RES_BOOL Interp,277 const RES_HWORD *Xp, RES_HWORD Ph, RES_HWORD Inc)278 {279 const RES_HWORD *Hp;280 const RES_HWORD *Hdp = NULL;281 const RES_HWORD *End;282 RES_HWORD a = 0;283 RES_WORD v, t;284 285 v=0;286 Hp = &Imp[Ph>>Na];287 End = &Imp[Nwing];288 if (Interp) {289 Hdp = &ImpD[Ph>>Na];290 a = Ph & Amask;291 }292 if (Inc == 1) /* If doing right wing... */293 { /* ...drop extra coeff, so when Ph is */294 End--; /* 0.5, we don't do too many mult's */295 if (Ph == 0) /* If the phase is zero... */296 { /* ...then we've already skipped the */297 Hp += Npc; /* first sample, so we must also */298 Hdp += Npc; /* skip ahead in Imp[] and ImpD[] */299 }300 }301 if (Interp)302 while (Hp < End) {303 t = *Hp; /* Get filter coeff */304 t += (((RES_WORD)*Hdp)*a)>>Na; /* t is now interp'd filter coeff */305 Hdp += Npc; /* Filter coeff differences step */306 t *= *Xp; /* Mult coeff by input sample */307 if (t & (1<<(Nhxn-1))) /* Round, if needed */308 t += (1<<(Nhxn-1));309 t >>= Nhxn; /* Leave some guard bits, but come back some */310 v += t; /* The filter output */311 Hp += Npc; /* Filter coeff step */312 313 Xp += Inc; /* Input signal step. NO CHECK ON BOUNDS */314 }315 else316 while (Hp < End) {317 t = *Hp; /* Get filter coeff */318 t *= *Xp; /* Mult coeff by input sample */319 if (t & (1<<(Nhxn-1))) /* Round, if needed */320 t += (1<<(Nhxn-1));321 t >>= Nhxn; /* Leave some guard bits, but come back some */322 v += t; /* The filter output */323 Hp += Npc; /* Filter coeff step */324 Xp += Inc; /* Input signal step. NO CHECK ON BOUNDS */325 }326 return(v);327 }328 329 330 static RES_WORD FilterUD(const RES_HWORD Imp[], const RES_HWORD ImpD[],331 RES_UHWORD Nwing, RES_BOOL Interp,332 const RES_HWORD *Xp, RES_HWORD Ph, RES_HWORD Inc, RES_UHWORD dhb)333 {334 RES_HWORD a;335 const RES_HWORD *Hp, *Hdp, *End;336 RES_WORD v, t;337 RES_UWORD Ho;338 339 v=0;340 Ho = (Ph*(RES_UWORD)dhb)>>Np;341 End = &Imp[Nwing];342 if (Inc == 1) /* If doing right wing... */343 { /* ...drop extra coeff, so when Ph is */344 End--; /* 0.5, we don't do too many mult's */345 if (Ph == 0) /* If the phase is zero... */346 Ho += dhb; /* ...then we've already skipped the */347 } /* first sample, so we must also */348 /* skip ahead in Imp[] and ImpD[] */349 if (Interp)350 while ((Hp = &Imp[Ho>>Na]) < End) {351 t = *Hp; /* Get IR sample */352 Hdp = &ImpD[Ho>>Na]; /* get interp (lower Na) bits from diff table*/353 a = Ho & Amask; /* a is logically between 0 and 1 */354 t += (((RES_WORD)*Hdp)*a)>>Na; /* t is now interp'd filter coeff */355 t *= *Xp; /* Mult coeff by input sample */356 if (t & 1<<(Nhxn-1)) /* Round, if needed */357 t += 1<<(Nhxn-1);358 t >>= Nhxn; /* Leave some guard bits, but come back some */359 v += t; /* The filter output */360 Ho += dhb; /* IR step */361 Xp += Inc; /* Input signal step. NO CHECK ON BOUNDS */362 }363 else364 while ((Hp = &Imp[Ho>>Na]) < End) {365 t = *Hp; /* Get IR sample */366 t *= *Xp; /* Mult coeff by input sample */367 if (t & 1<<(Nhxn-1)) /* Round, if needed */368 t += 1<<(Nhxn-1);369 t >>= Nhxn; /* Leave some guard bits, but come back some */370 v += t; /* The filter output */371 Ho += dhb; /* IR step */372 Xp += Inc; /* Input signal step. NO CHECK ON BOUNDS */373 }374 return(v);375 }376 377 /* Sampling rate up-conversion only subroutine;378 * Slightly faster than down-conversion;379 */380 static int SrcUp(const RES_HWORD X[], RES_HWORD Y[], double pFactor,381 RES_UHWORD nx, RES_UHWORD pNwing, RES_UHWORD pLpScl,382 const RES_HWORD pImp[], const RES_HWORD pImpD[], RES_BOOL Interp)383 {384 const RES_HWORD *xp;385 RES_HWORD *Ystart, *Yend;386 RES_WORD v;387 388 double dt; /* Step through input signal */389 RES_UWORD dtb; /* Fixed-point version of Dt */390 RES_UWORD time = 0;391 RES_UWORD endTime; /* When time reaches EndTime, return to user */392 393 dt = 1.0/pFactor; /* Output sampling period */394 dtb = dt*(1<<Np) + 0.5; /* Fixed-point representation */395 396 Ystart = Y;397 Yend = Ystart + (unsigned)(nx * pFactor);398 endTime = time + (1<<Np)*(RES_WORD)nx;399 while (time < endTime)400 {401 xp = &X[time>>Np]; /* Ptr to current input sample */402 /* Perform left-wing inner product */403 v = 0;404 v = FilterUp(pImp, pImpD, pNwing, Interp, xp, (RES_HWORD)(time&Pmask),-1);405 406 /* Perform right-wing inner product */407 v += FilterUp(pImp, pImpD, pNwing, Interp, xp+1, (RES_HWORD)((-time)&Pmask),1);408 409 v >>= Nhg; /* Make guard bits */410 v *= pLpScl; /* Normalize for unity filter gain */411 *Y++ = WordToHword(v,NLpScl); /* strip guard bits, deposit output */412 time += dtb; /* Move to next sample by time increment */413 }414 return (Y - Ystart); /* Return the number of output samples */415 }416 417 418 /* Sampling rate conversion subroutine */419 420 static int SrcUD(const RES_HWORD X[], RES_HWORD Y[], double pFactor,421 RES_UHWORD nx, RES_UHWORD pNwing, RES_UHWORD pLpScl,422 const RES_HWORD pImp[], const RES_HWORD pImpD[], RES_BOOL Interp)423 {424 const RES_HWORD *xp;425 RES_HWORD *Ystart, *Yend;426 RES_WORD v;427 428 double dh; /* Step through filter impulse response */429 double dt; /* Step through input signal */430 RES_UWORD time = 0;431 RES_UWORD endTime; /* When time reaches EndTime, return to user */432 RES_UWORD dhb, dtb; /* Fixed-point versions of Dh,Dt */433 434 dt = 1.0/pFactor; /* Output sampling period */435 dtb = dt*(1<<Np) + 0.5; /* Fixed-point representation */436 437 dh = MIN(Npc, pFactor*Npc); /* Filter sampling period */438 dhb = dh*(1<<Na) + 0.5; /* Fixed-point representation */439 440 Ystart = Y;441 Yend = Ystart + (unsigned)(nx * pFactor);442 endTime = time + (1<<Np)*(RES_WORD)nx;443 while (time < endTime)444 {445 xp = &X[time>>Np]; /* Ptr to current input sample */446 v = FilterUD(pImp, pImpD, pNwing, Interp, xp, (RES_HWORD)(time&Pmask),447 -1, dhb); /* Perform left-wing inner product */448 v += FilterUD(pImp, pImpD, pNwing, Interp, xp+1, (RES_HWORD)((-time)&Pmask),449 1, dhb); /* Perform right-wing inner product */450 v >>= Nhg; /* Make guard bits */451 v *= pLpScl; /* Normalize for unity filter gain */452 *Y++ = WordToHword(v,NLpScl); /* strip guard bits, deposit output */453 time += dtb; /* Move to next sample by time increment */454 }455 return (Y - Ystart); /* Return the number of output samples */456 }457 458 459 /* ***************************************************************************460 *461 * PJMEDIA RESAMPLE462 *463 * ***************************************************************************464 */465 30 466 31 struct pjmedia_resample … … 473 38 pj_int16_t *buffer; /* Input buffer. */ 474 39 }; 40 41 42 #ifndef MAX 43 # define MAX(a,b) ((a) >= (b) ? (a) : (b)) 44 #endif 475 45 476 46 … … 539 109 */ 540 110 if (large_filter) 541 resample->xoff = (LARGE_FILTER_NMULT + 1) / 2.0 * 542 MAX(1.0, 1.0/resample->factor); 111 resample->xoff = (unsigned) 112 ((resample_LARGE_FILTER_NMULT + 1) / 2.0 * 113 MAX(1.0, 1.0/resample->factor)); 543 114 else 544 resample->xoff = (SMALL_FILTER_NMULT + 1) / 2.0 * 545 MAX(1.0, 1.0/resample->factor); 115 resample->xoff = (unsigned) 116 ((resample_SMALL_FILTER_NMULT + 1) / 2.0 * 117 MAX(1.0, 1.0/resample->factor)); 546 118 547 119 … … 653 225 SrcUp(resample->buffer + resample->xoff, output, 654 226 resample->factor, resample->frame_size, 655 LARGE_FILTER_NWING,LARGE_FILTER_SCALE,656 LARGE_FILTER_IMP,LARGE_FILTER_IMPD,227 resample_LARGE_FILTER_NWING, resample_LARGE_FILTER_SCALE, 228 resample_LARGE_FILTER_IMP, resample_LARGE_FILTER_IMPD, 657 229 PJ_TRUE); 658 230 } else { 659 231 SrcUp(resample->buffer + resample->xoff, output, 660 232 resample->factor, resample->frame_size, 661 SMALL_FILTER_NWING,SMALL_FILTER_SCALE,662 SMALL_FILTER_IMP,SMALL_FILTER_IMPD,233 resample_SMALL_FILTER_NWING, resample_SMALL_FILTER_SCALE, 234 resample_SMALL_FILTER_IMP, resample_SMALL_FILTER_IMPD, 663 235 PJ_TRUE); 664 236 } … … 670 242 SrcUD( resample->buffer + resample->xoff, output, 671 243 resample->factor, resample->frame_size, 672 LARGE_FILTER_NWING,673 LARGE_FILTER_SCALE * resample->factor + 0.5,674 LARGE_FILTER_IMP,LARGE_FILTER_IMPD,244 resample_LARGE_FILTER_NWING, 245 resample_LARGE_FILTER_SCALE * resample->factor + 0.5, 246 resample_LARGE_FILTER_IMP, resample_LARGE_FILTER_IMPD, 675 247 PJ_TRUE); 676 248 … … 679 251 SrcUD( resample->buffer + resample->xoff, output, 680 252 resample->factor, resample->frame_size, 681 SMALL_FILTER_NWING,682 SMALL_FILTER_SCALE * resample->factor + 0.5,683 SMALL_FILTER_IMP,SMALL_FILTER_IMPD,253 resample_SMALL_FILTER_NWING, 254 resample_SMALL_FILTER_SCALE * resample->factor + 0.5, 255 resample_SMALL_FILTER_IMP, resample_SMALL_FILTER_IMPD, 684 256 PJ_TRUE); 685 257
Note: See TracChangeset
for help on using the changeset viewer.
