- Timestamp:
- Apr 9, 2007 7:06:08 AM (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
pjproject/branches/split-3rd-party/pjmedia/src/pjmedia/resample_resample.c
r1174 r1177 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> 21 63 22 #include <pjmedia/errno.h> 64 23 #include <pj/assert.h> … … 67 26 68 27 28 #if defined(PJMEDIA_HAS_LIBRESAMPLE) && PJMEDIA_HAS_LIBRESAMPLE != 0 29 30 #include <third_party/resample/include/resamplesubs.h> 31 69 32 #define THIS_FILE "resample.c" 70 33 71 34 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 35 466 36 struct pjmedia_resample … … 504 74 } 505 75 506 #if !defined(PJMEDIA_HAS_LARGE_FILTER) || PJMEDIA_HAS_LARGE_FILTER==0507 /*508 * If large filter is excluded in the build, then prevent application509 * from using it.510 */511 if (high_quality && large_filter) {512 large_filter = PJ_FALSE;513 PJ_LOG(5,(THIS_FILE,514 "Resample uses small filter because large filter is "515 "disabled"));516 }517 #endif518 519 #if !defined(PJMEDIA_HAS_SMALL_FILTER) || PJMEDIA_HAS_SMALL_FILTER==0520 /*521 * If small filter is excluded in the build and application wants to522 * use it, then drop to linear conversion.523 */524 if (high_quality && large_filter == 0) {525 high_quality = PJ_FALSE;526 PJ_LOG(4,(THIS_FILE,527 "Resample uses linear because small filter is disabled"));528 }529 #endif530 531 76 resample->factor = rate_out * 1.0 / rate_in; 532 77 resample->large_filter = large_filter; … … 541 86 * resample->xoff = large_filter ? 32 : 6; 542 87 */ 543 if (large_filter) 544 resample->xoff = (LARGE_FILTER_NMULT + 1) / 2.0 * 545 MAX(1.0, 1.0/resample->factor); 546 else 547 resample->xoff = (SMALL_FILTER_NMULT + 1) / 2.0 * 548 MAX(1.0, 1.0/resample->factor); 549 88 resample->xoff = res_GetXOFF(resample->factor, (char)large_filter); 550 89 551 90 size = (samples_per_frame + 2*resample->xoff) * sizeof(pj_int16_t); … … 650 189 dst_buf = resample->buffer + resample->xoff*2; 651 190 pjmedia_copy_samples(dst_buf, input, resample->frame_size); 652 653 if (resample->factor >= 1) { 654 655 if (resample->large_filter) { 656 SrcUp(resample->buffer + resample->xoff, output, 657 resample->factor, resample->frame_size, 658 LARGE_FILTER_NWING, LARGE_FILTER_SCALE, 659 LARGE_FILTER_IMP, LARGE_FILTER_IMPD, 660 PJ_TRUE); 661 } else { 662 SrcUp(resample->buffer + resample->xoff, output, 663 resample->factor, resample->frame_size, 664 SMALL_FILTER_NWING, SMALL_FILTER_SCALE, 665 SMALL_FILTER_IMP, SMALL_FILTER_IMPD, 666 PJ_TRUE); 667 } 668 669 } else { 670 671 if (resample->large_filter) { 672 673 SrcUD( resample->buffer + resample->xoff, output, 674 resample->factor, resample->frame_size, 675 LARGE_FILTER_NWING, 676 LARGE_FILTER_SCALE * resample->factor + 0.5, 677 LARGE_FILTER_IMP, LARGE_FILTER_IMPD, 678 PJ_TRUE); 679 680 } else { 681 682 SrcUD( resample->buffer + resample->xoff, output, 683 resample->factor, resample->frame_size, 684 SMALL_FILTER_NWING, 685 SMALL_FILTER_SCALE * resample->factor + 0.5, 686 SMALL_FILTER_IMP, SMALL_FILTER_IMPD, 687 PJ_TRUE); 688 689 } 690 691 } 191 192 res_Resample(resample->buffer + resample->xoff, output, 193 resample->factor, (pj_uint16_t)resample->frame_size, 194 (char)resample->large_filter, (char)PJ_TRUE); 692 195 693 196 dst_buf = resample->buffer; … … 696 199 697 200 } else { 698 SrcLinear( input, output, resample->factor, resample->frame_size); 201 res_SrcLinear( input, output, resample->factor, 202 (pj_uint16_t)resample->frame_size); 699 203 } 700 204 } … … 710 214 PJ_UNUSED_ARG(resample); 711 215 } 216 217 218 #else /* PJMEDIA_HAS_LIBRESAMPLE */ 219 220 int pjmedia_resample_resample_excluded; 221 222 #endif /* PJMEDIA_HAS_LIBRESAMPLE */ 223
Note: See TracChangeset
for help on using the changeset viewer.