2 * lame utility library source file
4 * Copyright (c) 1999 Albert L Faber
5 * Copyright (c) 2000-2005 Alexander Leidinger
7 * This library is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU Library General Public
9 * License as published by the Free Software Foundation; either
10 * version 2 of the License, or (at your option) any later version.
12 * This library is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Library General Public License for more details.
17 * You should have received a copy of the GNU Library General Public
18 * License along with this library; if not, write to the
19 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
20 * Boston, MA 02111-1307, USA.
23 /* $Id: util.c,v 1.154 2011/10/18 21:51:20 robert Exp $ */
36 #if defined(__FreeBSD__) && !defined(__alpha__)
37 # include <machine/floatingpoint.h>
41 /***********************************************************************
43 * Global Function Definitions
45 ***********************************************************************/
46 /*empty and close mallocs in gfc */
49 free_id3tag(lame_internal_flags * const gfc)
51 if (gfc->tag_spec.title != 0) {
52 free(gfc->tag_spec.title);
53 gfc->tag_spec.title = 0;
55 if (gfc->tag_spec.artist != 0) {
56 free(gfc->tag_spec.artist);
57 gfc->tag_spec.artist = 0;
59 if (gfc->tag_spec.album != 0) {
60 free(gfc->tag_spec.album);
61 gfc->tag_spec.album = 0;
63 if (gfc->tag_spec.comment != 0) {
64 free(gfc->tag_spec.comment);
65 gfc->tag_spec.comment = 0;
68 if (gfc->tag_spec.albumart != 0) {
69 free(gfc->tag_spec.albumart);
70 gfc->tag_spec.albumart = 0;
71 gfc->tag_spec.albumart_size = 0;
72 gfc->tag_spec.albumart_mimetype = MIMETYPE_NONE;
74 if (gfc->tag_spec.values != 0) {
76 for (i = 0; i < gfc->tag_spec.num_values; ++i) {
77 free(gfc->tag_spec.values[i]);
79 free(gfc->tag_spec.values);
80 gfc->tag_spec.values = 0;
81 gfc->tag_spec.num_values = 0;
83 if (gfc->tag_spec.v2_head != 0) {
84 FrameDataNode *node = gfc->tag_spec.v2_head;
86 void *p = node->dsc.ptr.b;
87 void *q = node->txt.ptr.b;
94 gfc->tag_spec.v2_head = 0;
95 gfc->tag_spec.v2_tail = 0;
101 free_global_data(lame_internal_flags * gfc)
103 if (gfc && gfc->cd_psy) {
104 if (gfc->cd_psy->l.s3) {
105 /* XXX allocated in psymodel_init() */
106 free(gfc->cd_psy->l.s3);
108 if (gfc->cd_psy->s.s3) {
109 /* XXX allocated in psymodel_init() */
110 free(gfc->cd_psy->s.s3);
119 freegfc(lame_internal_flags * const gfc)
120 { /* bit stream structure */
124 for (i = 0; i <= 2 * BPC; i++)
125 if (gfc->sv_enc.blackfilt[i] != NULL) {
126 free(gfc->sv_enc.blackfilt[i]);
127 gfc->sv_enc.blackfilt[i] = NULL;
129 if (gfc->sv_enc.inbuf_old[0]) {
130 free(gfc->sv_enc.inbuf_old[0]);
131 gfc->sv_enc.inbuf_old[0] = NULL;
133 if (gfc->sv_enc.inbuf_old[1]) {
134 free(gfc->sv_enc.inbuf_old[1]);
135 gfc->sv_enc.inbuf_old[1] = NULL;
138 if (gfc->bs.buf != NULL) {
143 if (gfc->VBR_seek_table.bag) {
144 free(gfc->VBR_seek_table.bag);
145 gfc->VBR_seek_table.bag = NULL;
146 gfc->VBR_seek_table.size = 0;
151 if (gfc->sv_rpg.rgdata) {
152 free(gfc->sv_rpg.rgdata);
154 if (gfc->sv_enc.in_buffer_0) {
155 free(gfc->sv_enc.in_buffer_0);
157 if (gfc->sv_enc.in_buffer_1) {
158 free(gfc->sv_enc.in_buffer_1);
162 #ifdef DECODE_ON_THE_FLY
164 hip_decode_exit(gfc->hip);
169 free_global_data(gfc);
174 /*those ATH formulas are returning
175 their minimum value for input = -1*/
178 ATHformula_GB(FLOAT f, FLOAT value, FLOAT f_min, FLOAT f_max)
180 /* from Painter & Spanias
181 modified by Gabriel Bouvigne to better fit the reality
182 ath = 3.640 * pow(f,-0.8)
183 - 6.800 * exp(-0.6*pow(f-3.4,2.0))
184 + 6.000 * exp(-0.15*pow(f-8.7,2.0))
185 + 0.6* 0.001 * pow(f,4.0);
188 In the past LAME was using the Painter &Spanias formula.
189 But we had some recurrent problems with HF content.
190 We measured real ATH values, and found the older formula
191 to be inacurate in the higher part. So we made this new
192 formula and this solved most of HF problematic testcases.
193 The tradeoff is that in VBR mode it increases a lot the
197 /*this curve can be udjusted according to the VBR scale:
198 it adjusts from something close to Painter & Spanias
199 on V9 up to Bouvigne's formula for V0. This way the VBR
200 bitrate is more balanced according to the -V value.*/
204 /* the following Hack allows to ask for the lowest value */
208 f /= 1000; /* convert to khz */
212 ath = 3.640 * pow(f, -0.8)
213 - 6.800 * exp(-0.6 * pow(f - 3.4, 2.0))
214 + 6.000 * exp(-0.15 * pow(f - 8.7, 2.0))
215 + (0.6 + 0.04 * value) * 0.001 * pow(f, 4.0);
222 ATHformula(SessionConfig_t const *cfg, FLOAT f)
225 switch (cfg->ATHtype) {
227 ath = ATHformula_GB(f, 9, 0.1f, 24.0f);
230 ath = ATHformula_GB(f, -1, 0.1f, 24.0f); /*over sensitive, should probably be removed */
233 ath = ATHformula_GB(f, 0, 0.1f, 24.0f);
236 ath = ATHformula_GB(f, 1, 0.1f, 24.0f) + 6; /*modification of GB formula by Roel */
239 ath = ATHformula_GB(f, cfg->ATHcurve, 0.1f, 24.0f);
242 ath = ATHformula_GB(f, cfg->ATHcurve, 3.41f, 16.1f);
245 ath = ATHformula_GB(f, 0, 0.1f, 24.0f);
251 /* see for example "Zwicker: Psychoakustik, 1982; ISBN 3-540-11401-7 */
253 freq2bark(FLOAT freq)
255 /* input: freq in hz output: barks */
259 return 13.0 * atan(.76 * freq) + 3.5 * atan(freq * freq / (7.5 * 7.5));
263 extern FLOAT freq2cbw(FLOAT freq);
265 /* see for example "Zwicker: Psychoakustik, 1982; ISBN 3-540-11401-7 */
269 /* input: freq in hz output: critical band width */
271 return 25 + 75 * pow(1 + 1.4 * (freq * freq), 0.69);
279 #define ABS(A) (((A)>0) ? (A) : -(A))
282 FindNearestBitrate(int bRate, /* legal rates from 8 to 320 */
283 int version, int samplerate)
284 { /* MPEG-1 or MPEG-2 LSF */
288 if (samplerate < 16000)
291 bitrate = bitrate_table[version][1];
293 for (i = 2; i <= 14; i++) {
294 if (bitrate_table[version][i] > 0) {
295 if (ABS(bitrate_table[version][i] - bRate) < ABS(bitrate - bRate))
296 bitrate = bitrate_table[version][i];
307 #define Min(A, B) ((A) < (B) ? (A) : (B))
310 #define Max(A, B) ((A) > (B) ? (A) : (B))
314 /* Used to find table index when
315 * we need bitrate-based values
316 * determined using tables
320 * Gabriel Bouvigne 2002-11-03
323 nearestBitrateFullIndex(uint16_t bitrate)
325 /* borrowed from DM abr presets */
327 const int full_bitrate_table[] =
328 { 8, 16, 24, 32, 40, 48, 56, 64, 80, 96, 112, 128, 160, 192, 224, 256, 320 };
331 int lower_range = 0, lower_range_kbps = 0, upper_range = 0, upper_range_kbps = 0;
337 /* We assume specified bitrate will be 320kbps */
338 upper_range_kbps = full_bitrate_table[16];
340 lower_range_kbps = full_bitrate_table[16];
343 /* Determine which significant bitrates the value specified falls between,
344 * if loop ends without breaking then we were correct above that the value was 320
346 for (b = 0; b < 16; b++) {
347 if ((Max(bitrate, full_bitrate_table[b + 1])) != bitrate) {
348 upper_range_kbps = full_bitrate_table[b + 1];
350 lower_range_kbps = full_bitrate_table[b];
352 break; /* We found upper range */
356 /* Determine which range the value specified is closer to */
357 if ((upper_range_kbps - bitrate) > (bitrate - lower_range_kbps)) {
367 /* map frequency to a valid MP3 sample frequency
369 * Robert Hegemann 2000-07-01
372 map2MP3Frequency(int freq)
395 BitrateIndex(int bRate, /* legal rates from 32 to 448 kbps */
396 int version, /* MPEG-1 or MPEG-2/2.5 LSF */
398 { /* convert bitrate in kbps to index */
400 if (samplerate < 16000)
402 for (i = 0; i <= 14; i++) {
403 if (bitrate_table[version][i] > 0) {
404 if (bitrate_table[version][i] == bRate) {
412 /* convert samp freq in Hz to index */
415 SmpFrqIndex(int sample_freq, int *const version)
417 switch (sample_freq) {
452 /*****************************************************************************
454 * End of bit_stream.c package
456 *****************************************************************************/
467 /* resampling via FIR filter, blackman window */
469 blackman(FLOAT x, FLOAT fcn, int l)
471 /* This algorithm from:
472 SIGNAL PROCESSING ALGORITHMS IN FORTRAN AND C
473 S.D. Stearns and R.A. David, Prentice-Hall, 1992
476 FLOAT const wcn = (PI * fcn);
485 bkwn = 0.42 - 0.5 * cos(2 * x * PI) + 0.08 * cos(4 * x * PI);
489 return (bkwn * sin(l * wcn * x2) / (PI * l * x2));
497 /* gcd - greatest common divisor */
498 /* Joint work of Euclid and M. Hendry */
503 /* assert ( i > 0 && j > 0 ); */
504 return j ? gcd(j, i % j) : i;
510 fill_buffer_resample(lame_internal_flags * gfc,
512 int desired_len, sample_t const *inbuf, int len, int *num_used, int ch)
514 SessionConfig_t const *const cfg = &gfc->cfg;
515 EncStateVar_t *esv = &gfc->sv_enc;
516 double resample_ratio = (double)cfg->samplerate_in / (double)cfg->samplerate_out;
518 FLOAT offset, xvalue;
523 int bpc; /* number of convolution functions to pre-compute */
524 bpc = cfg->samplerate_out / gcd(cfg->samplerate_out, cfg->samplerate_in);
528 intratio = (fabs(resample_ratio - floor(.5 + resample_ratio)) < .0001);
529 fcn = 1.00 / resample_ratio;
532 filter_l = 31; /* must be odd */
533 filter_l += intratio; /* unless resample_ratio=int, it must be even */
536 BLACKSIZE = filter_l + 1; /* size of data needed for FIR */
538 if (gfc->fill_buffer_resample_init == 0) {
539 esv->inbuf_old[0] = calloc(BLACKSIZE, sizeof(esv->inbuf_old[0][0]));
540 esv->inbuf_old[1] = calloc(BLACKSIZE, sizeof(esv->inbuf_old[0][0]));
541 for (i = 0; i <= 2 * bpc; ++i)
542 esv->blackfilt[i] = calloc(BLACKSIZE, sizeof(esv->blackfilt[0][0]));
547 /* precompute blackman filter coefficients */
548 for (j = 0; j <= 2 * bpc; j++) {
550 offset = (j - bpc) / (2. * bpc);
551 for (i = 0; i <= filter_l; i++)
552 sum += esv->blackfilt[j][i] = blackman(i - offset, fcn, filter_l);
553 for (i = 0; i <= filter_l; i++)
554 esv->blackfilt[j][i] /= sum;
556 gfc->fill_buffer_resample_init = 1;
559 inbuf_old = esv->inbuf_old[ch];
561 /* time of j'th element in inbuf = itime + j/ifreq; */
562 /* time of k'th element in outbuf = j/ofreq */
563 for (k = 0; k < desired_len; k++) {
564 double time0 = k * resample_ratio; /* time of k'th output sample */
567 j = floor(time0 - esv->itime[ch]);
569 /* check if we need more input data */
570 if ((filter_l + j - filter_l / 2) >= len)
573 /* blackman filter. by default, window centered at j+.5(filter_l%2) */
574 /* but we want a window centered at time0. */
575 offset = (time0 - esv->itime[ch] - (j + .5 * (filter_l % 2)));
576 assert(fabs(offset) <= .501);
578 /* find the closest precomputed window for this offset: */
579 joff = floor((offset * 2 * bpc) + bpc + .5);
582 for (i = 0; i <= filter_l; ++i) {
583 int const j2 = i + j - filter_l / 2;
586 assert(j2 + BLACKSIZE >= 0);
587 y = (j2 < 0) ? inbuf_old[BLACKSIZE + j2] : inbuf[j2];
589 xvalue += y * esv->blackfilt[joff][i];
591 xvalue += y * blackman(i - offset, fcn, filter_l); /* very slow! */
598 /* k = number of samples added to outbuf */
599 /* last k sample used data from [j-filter_l/2,j+filter_l-filter_l/2] */
601 /* how many samples of input data were used: */
602 *num_used = Min(len, filter_l + j - filter_l / 2);
604 /* adjust our input time counter. Incriment by the number of samples used,
605 * then normalize so that next output sample is at time 0, next
606 * input buffer is at time itime[ch] */
607 esv->itime[ch] += *num_used - k * resample_ratio;
609 /* save the last BLACKSIZE samples into the inbuf_old buffer */
610 if (*num_used >= BLACKSIZE) {
611 for (i = 0; i < BLACKSIZE; i++)
612 inbuf_old[i] = inbuf[*num_used + i - BLACKSIZE];
615 /* shift in *num_used samples into inbuf_old */
616 int const n_shift = BLACKSIZE - *num_used; /* number of samples to shift */
618 /* shift n_shift samples by *num_used, to make room for the
619 * num_used new samples */
620 for (i = 0; i < n_shift; ++i)
621 inbuf_old[i] = inbuf_old[i + *num_used];
623 /* shift in the *num_used samples */
624 for (j = 0; i < BLACKSIZE; ++i, ++j)
625 inbuf_old[i] = inbuf[j];
627 assert(j == *num_used);
629 return k; /* return the number samples created at the new samplerate */
633 isResamplingNecessary(SessionConfig_t const* cfg)
635 int const l = cfg->samplerate_out * 0.9995f;
636 int const h = cfg->samplerate_out * 1.0005f;
637 return (cfg->samplerate_in < l) || (h < cfg->samplerate_in) ? 1 : 0;
640 /* copy in new samples from in_buffer into mfbuf, with resampling
641 if necessary. n_in = number of samples from the input buffer that
642 were used. n_out = number of samples copied into mfbuf */
645 fill_buffer(lame_internal_flags * gfc,
646 sample_t * const mfbuf[2], sample_t const * const in_buffer[2], int nsamples, int *n_in, int *n_out)
648 SessionConfig_t const *const cfg = &gfc->cfg;
649 int mf_size = gfc->sv_enc.mf_size;
650 int framesize = 576 * cfg->mode_gr;
652 int nch = cfg->channels_out;
654 /* copy in new samples into mfbuf, with resampling if necessary */
655 if (isResamplingNecessary(cfg)) {
658 fill_buffer_resample(gfc, mfbuf[ch]+mf_size,
659 framesize, in_buffer[ch], nsamples, n_in, ch);
660 } while (++ch < nch);
664 nout = Min(framesize, nsamples);
666 memcpy(mfbuf[ch]+mf_size,in_buffer[ch],sizeof(sample_t) * nout);
667 } while (++ch < nch);
679 /***********************************************************************
683 ***********************************************************************/
686 lame_report_def(const char *format, va_list args)
688 (void) vfprintf(stderr, format, args);
689 fflush(stderr); /* an debug function should flush immediately */
693 lame_report_fnc(lame_report_function print_f, const char *format, ...)
697 va_start(args, format);
698 print_f(format, args);
705 lame_debugf(const lame_internal_flags* gfc, const char *format, ...)
707 if (gfc && gfc->report_dbg) {
709 va_start(args, format);
710 gfc->report_dbg(format, args);
717 lame_msgf(const lame_internal_flags* gfc, const char *format, ...)
719 if (gfc && gfc->report_msg) {
721 va_start(args, format);
722 gfc->report_msg(format, args);
729 lame_errorf(const lame_internal_flags* gfc, const char *format, ...)
731 if (gfc && gfc->report_err) {
733 va_start(args, format);
734 gfc->report_err(format, args);
741 /***********************************************************************
743 * routines to detect CPU specific features like 3DNow, MMX, SSE
745 * donated by Frank Klemm
746 * added Robert Hegemann 2000-10-10
748 ***********************************************************************/
751 extern int has_MMX_nasm(void);
752 extern int has_3DNow_nasm(void);
753 extern int has_SSE_nasm(void);
754 extern int has_SSE2_nasm(void);
761 return has_MMX_nasm();
763 return 0; /* don't know, assume not */
771 return has_3DNow_nasm();
773 return 0; /* don't know, assume not */
781 return has_SSE_nasm();
783 #if defined( _M_X64 ) || defined( MIN_ARCH_SSE )
786 return 0; /* don't know, assume not */
795 return has_SSE2_nasm();
797 #if defined( _M_X64 ) || defined( MIN_ARCH_SSE )
800 return 0; /* don't know, assume not */
808 /* extremly system dependent stuff, move to a lib to make the code readable */
809 /*==========================================================================*/
814 * Disable floating point exceptions
820 #if defined(__FreeBSD__) && !defined(__alpha__)
822 /* seet floating point mask to the Linux default */
825 /* if bit is set, we get SIGFPE on that error! */
826 fpsetmask(mask & ~(FP_X_INV | FP_X_DZ));
827 /* DEBUGF("FreeBSD mask is 0x%x\n",mask); */
831 #if defined(__riscos__) && !defined(ABORTFP)
832 /* Disable FPE's under RISC OS */
833 /* if bit is set, we disable trapping that error! */
834 /* _FPE_IVO : invalid operation */
835 /* _FPE_DVZ : divide by zero */
836 /* _FPE_OFL : overflow */
837 /* _FPE_UFL : underflow */
838 /* _FPE_INX : inexact */
839 DisableFPETraps(_FPE_IVO | _FPE_DVZ | _FPE_OFL);
844 * The default is to ignore FPE's, unless compiled with -DABORTFP
845 * so add code below to ENABLE FPE's.
849 #if defined(_MSC_VER)
853 the following fix seems to be a workaround for a problem in the
854 parent process calling LAME. It would be better to fix the broken
855 application => code disabled.
858 /* set affinity to a single CPU. Fix for EAC/lame on SMP systems from
859 "Todd Richmond" <todd.richmond@openwave.com> */
862 SetProcessAffinityMask(GetCurrentProcess(), si.dwActiveProcessorMask);
866 mask = _controlfp(0, 0);
867 mask &= ~(_EM_OVERFLOW | _EM_UNDERFLOW | _EM_ZERODIVIDE | _EM_INVALID);
868 mask = _controlfp(mask, _MCW_EM);
870 #elif defined(__CYGWIN__)
871 # define _FPU_GETCW(cw) __asm__ ("fnstcw %0" : "=m" (*&cw))
872 # define _FPU_SETCW(cw) __asm__ ("fldcw %0" : : "m" (*&cw))
874 # define _EM_INEXACT 0x00000020 /* inexact (precision) */
875 # define _EM_UNDERFLOW 0x00000010 /* underflow */
876 # define _EM_OVERFLOW 0x00000008 /* overflow */
877 # define _EM_ZERODIVIDE 0x00000004 /* zero divide */
878 # define _EM_INVALID 0x00000001 /* invalid */
882 /* Set the FPU control word to abort on most FPEs */
883 mask &= ~(_EM_OVERFLOW | _EM_ZERODIVIDE | _EM_INVALID);
886 # elif defined(__linux__)
889 # include <fpu_control.h>
891 # define _FPU_GETCW(cw) __asm__ ("fnstcw %0" : "=m" (*&cw))
894 # define _FPU_SETCW(cw) __asm__ ("fldcw %0" : : "m" (*&cw))
898 * Set the Linux mask to abort on most FPE's
899 * if bit is set, we _mask_ SIGFPE on that error!
900 * mask &= ~( _FPU_MASK_IM | _FPU_MASK_ZM | _FPU_MASK_OM | _FPU_MASK_UM );
905 mask &= ~(_FPU_MASK_IM | _FPU_MASK_ZM | _FPU_MASK_OM);
917 /***********************************************************************
919 * Fast Log Approximation for log2, used to approximate every other log
921 * maximum absolute error for log10 is around 10-6
922 * maximum *relative* error can be high when x is almost 1 because error/log10(x) tends toward x/e
924 * use it if typical RESULT values are > 1e-5 (for example if x>1.00001 or x<0.99999)
925 * or if the relative precision in the domain around 1 is not important (result in 1 is exact and 0)
927 ***********************************************************************/
930 #define LOG2_SIZE (512)
931 #define LOG2_SIZE_L2 (9)
933 static ieee754_float32_t log_table[LOG2_SIZE + 1];
943 /* Range for log2(x) over [1,2[ is [0,1[ */
944 assert((1 << LOG2_SIZE_L2) == LOG2_SIZE);
947 for (j = 0; j < LOG2_SIZE + 1; j++)
948 log_table[j] = log(1.0f + j / (ieee754_float32_t) LOG2_SIZE) / log(2.0f);
956 fast_log2(ieee754_float32_t x)
958 ieee754_float32_t log2val, partial;
965 mantisse = fi.i & 0x7fffff;
966 log2val = ((fi.i >> 23) & 0xFF) - 0x7f;
967 partial = (mantisse & ((1 << (23 - LOG2_SIZE_L2)) - 1));
968 partial *= 1.0f / ((1 << (23 - LOG2_SIZE_L2)));
971 mantisse >>= (23 - LOG2_SIZE_L2);
973 /* log2val += log_table[mantisse]; without interpolation the results are not good */
974 log2val += log_table[mantisse] * (1.0f - partial) + log_table[mantisse + 1] * partial;
979 #else /* Don't use FAST_LOG */