xref: /plan9/sys/src/games/mp3enc/pcm.c (revision 8f5875f3e9b20916b4c52ad4336922bc8653eb7b)
1 /* -*- mode: C; mode: fold -*- */
2 
3 /* $Id: pcm.c,v 1.7 2001/02/17 14:30:56 aleidinger Exp $ */
4 
5 /*
6  *  There are a lot of not tested return codes.
7  *  This is currently intention to make the code more readable
8  *  in the design phase.
9  *  Keine Ber�cksichtigung des Ein- und Ausschwingens des Downsamplefilters am
10  *  Anfang und am Ende wie bei meinem Resample-Programm
11  */
12 
13 #ifdef HAVE_CONFIG_H
14 # include <config.h>
15 #endif
16 
17 #include "bitstream.h"
18 #include "id3tag.h"
19 
20 #define FN(x)	do { static unsigned int cnt = 0; if (cnt < 100000) fprintf (stderr,"[%3u]>>>%s<<<\n",++cnt,(x)), fflush(stderr); } while (0)
21 
22 #ifdef KLEMM_44
23 
24 /*{{{ #includes                       */
25 
26 #include <assert.h>
27 #include <stddef.h>
28 #include <stdio.h>
29 #include <stdlib.h>
30 #include <limits.h>
31 #include <math.h>
32 #include <memory.h>
33 #include <unistd.h>
34 #include "pcm.h"
35 
36 /*}}}*/
37 /*{{{ bitstream object                */
38 
39 /*
40  *  About handling of octetstreams:
41  *
42  *  CHAR_BIT = 8: data points to a normal memory block which can be written to
43  *                disk via fwrite
44  *  CHAR_BIT > 8: On every char only the bits from 0...7 are used. Depending on the
45  *                host I/O it can be possible that you must first pack the data
46  *                stream before writing to disk/network/...
47  *  CHAR_BIT < 8: Not allowed by C
48  */
49 
octetstream_open(OUT size_t size)50 octetstream_t*  octetstream_open ( OUT size_t  size )
51 {
52     octetstream_t*  ret = (octetstream_t*) calloc ( 1, sizeof(octetstream_t) );
53     FN("octetstream_open");
54     ret -> data = (uint8_t*) calloc ( 1, size );
55     ret -> size = size;
56     return ret;
57 }
58 
octetstream_resize(INOUT octetstream_t * const os,OUT size_t size)59 int  octetstream_resize ( INOUT octetstream_t* const  os, OUT size_t  size )
60 {
61     FN("octetstream_resize");
62     if ( size > os->size ) {
63         os -> data = (uint8_t*) realloc ( os -> data, size );
64         memset ( os -> data + os -> size, 0, size - os -> size );
65         os -> size = size;
66     } else if ( size < os->size ) {
67 	os -> data = (uint8_t*) realloc ( os -> data, os->size = size );
68     }
69     return 0;
70 }
71 
octetstream_close(INOUT octetstream_t * const os)72 int  octetstream_close ( INOUT octetstream_t* const  os )
73 {
74     FN("octetstream_close");
75     if ( os == NULL )
76         return -1;
77     if ( os -> data != NULL ) {
78         free (os -> data);
79         os -> data = NULL;
80         free (os);
81         return 0;
82     } else {
83         free (os);
84         return -1;
85     }
86 }
87 
88 /*}}}*/
89 /*{{{ encode one frame                */
90 
lame_encode_frame(INOUT lame_t * lame,OUTTR sample_t ** inbuf,IN uint8_t * mp3buf,OUT size_t mp3buf_size)91 static inline int  lame_encode_frame (
92         INOUT lame_t*     lame,
93         OUTTR sample_t**  inbuf,
94         IN    uint8_t*    mp3buf,
95         OUT   size_t      mp3buf_size )
96 {
97     int  ret;
98 #if 1
99     int           i;
100     static int    j  = 0;
101     static FILE*  fp = NULL;
102     if ( fp == NULL )
103         fp = fopen ("pcm_data.txt", "w");
104     for ( i = 0; i < lame->frame_size; i++ )
105         fprintf ( fp, "%7d %11.4f %11.4f\n", j++, inbuf[0][i], inbuf[1][i] );
106     fprintf ( fp, "\n" );
107     fflush ( fp );
108 #endif
109 
110     FN("lame_encode_frame");
111 
112     switch ( lame -> coding ) {
113 #ifdef HAVE_MPEG_LAYER1
114     case coding_MPEG_Layer_1:
115         ret = lame_encode_mp1_frame ( lame->global_flags, inbuf[0], inbuf[1], mp3buf, mp3buf_size );
116         break;
117 #endif
118 #ifdef HAVE_MPEG_LAYER2
119     case coding_MPEG_Layer_2:
120         ret = lame_encode_mp2_frame ( lame->global_flags, inbuf[0], inbuf[1], mp3buf, mp3buf_size );
121         break;
122 #endif
123     case coding_MPEG_Layer_3:
124         ret = lame_encode_mp3_frame ( lame->global_flags, inbuf[0], inbuf[1], mp3buf, mp3buf_size );
125         break;
126 #ifdef HAVE_MPEG_PLUS
127     case coding_MPEG_plus:
128         ret = lame_encode_mpp_frame ( lame->global_flags, inbuf[0], inbuf[1], mp3buf, mp3buf_size );
129         break;
130 #endif
131 #ifdef HAVE_AAC
132     case coding_MPEG_AAC:
133         ret = lame_encode_aac_frame ( lame->global_flags, inbuf[0], inbuf[1], mp3buf, mp3buf_size );
134         break;
135 #endif
136 #ifdef HAVE_VORBIS
137     case coding_Ogg_Vorbis:
138         ret = lame_encode_ogg_frame ( lame->global_flags, inbuf[0], inbuf[1], mp3buf, mp3buf_size );
139         break;
140 #endif
141     default:
142         ret = -5;
143         break;
144     }
145 
146     if ( ret >= 0 ) {
147         lame->frame_count++;
148         if ( lame->analyzer_callback != NULL )
149             lame->analyzer_callback ( lame, lame->frame_size );
150     }
151     return ret;
152 }
153 
154 /*}}}*/
155 /*{{{ demultiplexing tools            */
156 
157 /*
158  * Now there are following the so called peek functions. They got a pointer and returning
159  * the data to this pointer as a sample. The size of the memory object can be from
160  * 8 up to 80 bits. One sample must be fully determined by it's bits, not by the history
161  * or by other channels or things like that. Also the input must have a integral byte size.
162  *
163  * That means:
164  *
165  *   - 4 or 6 bit PCM can't be decoded
166  *   - APCM can't be decoded
167  *   - ulaw/alaw *can* be decoded
168  *
169  * Note: In the future there will be a SMART define supported which only
170  *       support native endian shorts.
171  */
172 
173 static const int16_t  ulaw [256] = {
174     -32124, -31100, -30076, -29052, -28028, -27004, -25980, -24956,
175     -23932, -22908, -21884, -20860, -19836, -18812, -17788, -16764,
176     -15996, -15484, -14972, -14460, -13948, -13436, -12924, -12412,
177     -11900, -11388, -10876, -10364,  -9852,  -9340,  -8828,  -8316,
178      -7932,  -7676,  -7420,  -7164,  -6908,  -6652,  -6396,  -6140,
179      -5884,  -5628,  -5372,  -5116,  -4860,  -4604,  -4348,  -4092,
180      -3900,  -3772,  -3644,  -3516,  -3388,  -3260,  -3132,  -3004,
181      -2876,  -2748,  -2620,  -2492,  -2364,  -2236,  -2108,  -1980,
182      -1884,  -1820,  -1756,  -1692,  -1628,  -1564,  -1500,  -1436,
183      -1372,  -1308,  -1244,  -1180,  -1116,  -1052,   -988,   -924,
184       -876,   -844,   -812,   -780,   -748,   -716,   -684,   -652,
185       -620,   -588,   -556,   -524,   -492,   -460,   -428,   -396,
186       -372,   -356,   -340,   -324,   -308,   -292,   -276,   -260,
187       -244,   -228,   -212,   -196,   -180,   -164,   -148,   -132,
188       -120,   -112,   -104,    -96,    -88,    -80,    -72,    -64,
189        -56,    -48,    -40,    -32,    -24,    -16,     -8,      0,
190      32124,  31100,  30076,  29052,  28028,  27004,  25980,  24956,
191      23932,  22908,  21884,  20860,  19836,  18812,  17788,  16764,
192      15996,  15484,  14972,  14460,  13948,  13436,  12924,  12412,
193      11900,  11388,  10876,  10364,   9852,   9340,   8828,   8316,
194       7932,   7676,   7420,   7164,   6908,   6652,   6396,   6140,
195       5884,   5628,   5372,   5116,   4860,   4604,   4348,   4092,
196       3900,   3772,   3644,   3516,   3388,   3260,   3132,   3004,
197       2876,   2748,   2620,   2492,   2364,   2236,   2108,   1980,
198       1884,   1820,   1756,   1692,   1628,   1564,   1500,   1436,
199       1372,   1308,   1244,   1180,   1116,   1052,    988,    924,
200        876,    844,    812,    780,    748,    716,    684,    652,
201        620,    588,    556,    524,    492,    460,    428,    396,
202        372,    356,    340,    324,    308,    292,    276,    260,
203        244,    228,    212,    196,    180,    164,    148,    132,
204        120,    112,    104,     96,     88,     80,     72,     64,
205         56,     48,     40,     32,     24,     16,      8,      0,
206 };
207 
208 static const int16_t  alaw [256] = {
209      -5504,  -5248,  -6016,  -5760,  -4480,  -4224,  -4992,  -4736,
210      -7552,  -7296,  -8064,  -7808,  -6528,  -6272,  -7040,  -6784,
211      -2752,  -2624,  -3008,  -2880,  -2240,  -2112,  -2496,  -2368,
212      -3776,  -3648,  -4032,  -3904,  -3264,  -3136,  -3520,  -3392,
213     -22016, -20992, -24064, -23040, -17920, -16896, -19968, -18944,
214     -30208, -29184, -32256, -31232, -26112, -25088, -28160, -27136,
215     -11008, -10496, -12032, -11520,  -8960,  -8448,  -9984,  -9472,
216     -15104, -14592, -16128, -15616, -13056, -12544, -14080, -13568,
217       -344,   -328,   -376,   -360,   -280,   -264,   -312,   -296,
218       -472,   -456,   -504,   -488,   -408,   -392,   -440,   -424,
219        -88,    -72,   -120,   -104,    -24,     -8,    -56,    -40,
220       -216,   -200,   -248,   -232,   -152,   -136,   -184,   -168,
221      -1376,  -1312,  -1504,  -1440,  -1120,  -1056,  -1248,  -1184,
222      -1888,  -1824,  -2016,  -1952,  -1632,  -1568,  -1760,  -1696,
223       -688,   -656,   -752,   -720,   -560,   -528,   -624,   -592,
224       -944,   -912,  -1008,   -976,   -816,   -784,   -880,   -848,
225       5504,   5248,   6016,   5760,   4480,   4224,   4992,   4736,
226       7552,   7296,   8064,   7808,   6528,   6272,   7040,   6784,
227       2752,   2624,   3008,   2880,   2240,   2112,   2496,   2368,
228       3776,   3648,   4032,   3904,   3264,   3136,   3520,   3392,
229      22016,  20992,  24064,  23040,  17920,  16896,  19968,  18944,
230      30208,  29184,  32256,  31232,  26112,  25088,  28160,  27136,
231      11008,  10496,  12032,  11520,   8960,   8448,   9984,   9472,
232      15104,  14592,  16128,  15616,  13056,  12544,  14080,  13568,
233        344,    328,    376,    360,    280,    264,    312,    296,
234        472,    456,    504,    488,    408,    392,    440,    424,
235         88,     72,    120,    104,     24,      8,     56,     40,
236        216,    200,    248,    232,    152,    136,    184,    168,
237       1376,   1312,   1504,   1440,   1120,   1056,   1248,   1184,
238       1888,   1824,   2016,   1952,   1632,   1568,   1760,   1696,
239        688,    656,    752,    720,    560,    528,    624,    592,
240        944,    912,   1008,    976,    816,    784,    880,    848,
241 };
242 
243 /*
244  *  These macros get 1...4 bytes and calculating a value between -32768 and +32768.
245  *  The first argument of the macros is the MSB, the last the LSB.
246  *  These are used to decode integer PCM (if there is no faster code available).
247  */
248 
249 #define BYTES1(x1)          ((((int8_t)(x1)) << 8))
250 #define BYTES2(x1,x2)       ((((int8_t)(x1)) << 8) + (uint8_t)(x2))
251 #define BYTES3(x1,x2,x3)    ((((int8_t)(x1)) << 8) + (uint8_t)(x2) + 1/256.*(uint8_t)(x3))
252 #define BYTES4(x1,x2,x3,x4) ((((int8_t)(x1)) << 8) + (uint8_t)(x2) + 1/256.*(uint8_t)(x3) + 1/65536.*(uint8_t)(x4))
253 
254 /*
255  * The next two functions can read a 32 and 64 bit floating pointer number
256  * of an opposite byte order machine. This is done by byte swaping and using
257  * the normal way to read such a variable.
258  * Also note: peek function got their data from a variable called src, which is a pointer
259  * to const unsigned char.
260  */
261 
read_opposite_float32(OUT uint8_t * const src)262 static inline sample_t  read_opposite_float32 ( OUT uint8_t* const src )
263 {
264     uint8_t  tmp [4];
265     tmp[0] = src[3];
266     tmp[1] = src[2];
267     tmp[2] = src[1];
268     tmp[3] = src[0];
269     return *(const float32_t*)tmp;
270 }
271 
read_opposite_float64(OUT uint8_t * const src)272 static inline sample_t  read_opposite_float64 ( OUT uint8_t* const src )
273 {
274     uint8_t  tmp [8];
275     tmp[0] = src[7];
276     tmp[1] = src[6];
277     tmp[2] = src[5];
278     tmp[3] = src[4];
279     tmp[4] = src[3];
280     tmp[5] = src[2];
281     tmp[6] = src[1];
282     tmp[7] = src[0];
283     return *(const float64_t*)tmp;
284 }
285 
286 /*
287  *  The next two functions can read a 80 bit extended precision
288  *  IEEE-854 floating point number. This code also can read the numbers
289  *  on a machine not supporting 80 bit floats.
290  */
291 
read_float80_le(OUT uint8_t * const src)292 static sample_t  read_float80_le ( OUT uint8_t* const src )
293 {
294     int         sign     = src [9] & 128  ?  -1  :  +1;
295     long        exponent =(src [9] & 127) * 256 +
296                            src [8] - 16384 - 62;
297     long double mantisse = src [7] * 72057594037927936.L +
298                            src [6] * 281474976710656.L +
299                            src [5] * 1099511627776.L +
300                            src [4] * 4294967296.L +
301                            src [3] * 16777216.L +
302                            src [2] * 65536.L +
303                            src [1] * 256.L +
304                            src [0];
305     return sign * ldexp (mantisse, exponent);
306 }
307 
read_float80_be(OUT uint8_t * const src)308 static sample_t  read_float80_be ( OUT uint8_t* const src )
309 {
310     int         sign     = src [0] & 128  ?  -1  :  +1;
311     long        exponent =(src [0] & 127) * 256 +
312                            src [1] - 16384 - 62;
313     long double mantisse = src [2] * 72057594037927936.L +
314                            src [3] * 281474976710656.L +
315                            src [4] * 1099511627776.L +
316                            src [5] * 4294967296.L +
317                            src [6] * 16777216.L +
318                            src [7] * 65536.L +
319                            src [8] * 256.L +
320                            src [9];
321     return sign * ldexp (mantisse, exponent);
322 }
323 
read_opposite_float80(OUT uint8_t * const src)324 static inline sample_t  read_opposite_float80 ( OUT uint8_t* const src )
325 {
326     uint8_t  tmp [10];
327     tmp[0] = src[9];
328     tmp[1] = src[8];
329     tmp[2] = src[7];
330     tmp[3] = src[6];
331     tmp[4] = src[5];
332     tmp[5] = src[4];
333     tmp[6] = src[3];
334     tmp[7] = src[2];
335     tmp[8] = src[1];
336     tmp[9] = src[0];
337     return *(const float80_t*)tmp;
338 }
339 
340 
341 /*
342  *  Now we are building all resorting functions. The macro wants the name of the function
343  *  and a peek function or macro. The created function wants the following input:
344  *
345  *     destination: Where to store the result as sample_t vector? The distance between
346  *                  the elements is sizeof(sample_t)
347  *     source:      Where to get the bytes which should converted to destination
348  *     step size:   The distance between the input samples. This is at least the size
349  *                  of one input element, but can be more due to
350  *                    - alignment
351  *                    - interleaved multi-channel input
352  *     length:      The number of elements to process. Number must be positive.
353  */
354 
355 typedef void (*demux_t) ( IN sample_t* dst, OUT uint8_t* src, OUT ssize_t step, OUT size_t len );
356 
357 #define FUNCTION(name,expr)       \
358 static void  name (               \
359         IN  sample_t*  dst,       \
360 	OUT uint8_t*   src,       \
361 	OUT ssize_t    step,      \
362 	OUT size_t     len )      \
363 {                                 \
364     size_t  i = len;              \
365     do {                          \
366         *dst++ = (expr);          \
367         src   += (step);          \
368     } while (--i);                \
369 }
370 
371 /* Some of these function can be implemented byte order independent ... */
372 
373 FUNCTION ( copy_silence, 0 )
374 FUNCTION ( copy_u8     , BYTES1 (src[0]-128) )
375 FUNCTION ( copy_s8     , BYTES1 (src[0]    ) )
376 FUNCTION ( copy_alaw   , alaw[*src] )
377 FUNCTION ( copy_ulaw   , ulaw[*src] )
378 FUNCTION ( copy_s24_le , BYTES3 (src[2], src[1], src[0] ) )
379 FUNCTION ( copy_s24_be , BYTES3 (src[0], src[1], src[2] ) )
380 FUNCTION ( copy_u24_le , BYTES3 (src[2]-128, src[1], src[0] ) )
381 FUNCTION ( copy_u24_be , BYTES3 (src[0]-128, src[1], src[2] ) )
382 FUNCTION ( copy_f80_le , read_float80_le (src) )
383 FUNCTION ( copy_f80_be , read_float80_be (src) )
384 
385 /* ... and some are not or there are faster but byte order dependent possiblities */
386 
387 
388 #ifndef WORDS_BIGENDIAN
389 
390 /* little endian */
391 FUNCTION ( copy_s16_le , *(const int16_t*)src )
392 FUNCTION ( copy_s16_be , BYTES2 (src[0], src[1] ) )
393 FUNCTION ( copy_u16_le , *(const uint16_t*)src - 32768 )
394 FUNCTION ( copy_u16_be , BYTES2 (src[0]-128, src[1] ) )
395 FUNCTION ( copy_s32_le , 1/65536. * *(const int32_t*)src )
396 FUNCTION ( copy_s32_be , BYTES4 (src[0], src[1], src[2], src[3] ) )
397 FUNCTION ( copy_u32_le , 1/65536. * *(const uint32_t*)src - 32768. )
398 FUNCTION ( copy_u32_be , BYTES4 (src[0]-128, src[1], src[2], src[3] ) )
399 FUNCTION ( copy_f32_le , *(const float32_t*)src )
400 FUNCTION ( copy_f32_be , read_opposite_float32 (src) )
401 FUNCTION ( copy_f64_le , *(const float64_t*)src )
402 FUNCTION ( copy_f64_be , read_opposite_float64 (src) )
403 
404 #else
405 
406 /* big endian */
407 FUNCTION ( copy_s16_le , BYTES2 (src[1], src[0] ) )
408 FUNCTION ( copy_s16_be , *(const int16_t*)src )
409 FUNCTION ( copy_u16_le , BYTES2 (src[1]-128, src[0] ) )
410 FUNCTION ( copy_u16_be , *(const uint16_t*)src - 32768 )
411 FUNCTION ( copy_s32_le , BYTES4 (src[3], src[2], src[1], src[0] ) )
412 FUNCTION ( copy_s32_be , 1/65536. * *(const int32_t*)src )
413 FUNCTION ( copy_u32_le , BYTES4 (src[3]-128, src[2], src[1], src[0] ) )
414 FUNCTION ( copy_u32_be , 1/65536. * *(const uint32_t*)src - 32768. )
415 FUNCTION ( copy_f32_le , read_opposite_float32 (src) )
416 FUNCTION ( copy_f32_be , *(const float32_t*)src )
417 FUNCTION ( copy_f64_le , read_opposite_float64 (src) )
418 FUNCTION ( copy_f64_be , *(const float64_t*)src )
419 
420 #endif
421 
422 /*
423  *  The global collection of channel demultiplexer.
424  *  For every data type we have the size of one memory object and two
425  *  demultiplexer, one for big endians, one for little endians.
426  *  The demultiplexers
427  *    - 1st argument is the destination of the data converted to sample_t (and normalized to +/-32768)
428  *    - 2nd argument is a pointer to the source data, in any format
429  *    - 3rd argument is 1st argument of the table (multiplied by the channel count in the case of
430  *      interleaved data)
431  *    - 4th argument is the count of samples per channel to convert
432  */
433 
434 typedef struct {
435     const ssize_t  size;
436     const demux_t  demultiplexer_be;
437     const demux_t  demultiplexer_le;
438 } demux_info_t;
439 
440 const demux_info_t  demux_info [] = {
441     { -1, NULL        , NULL         }, /* 00: */
442     {  0, copy_silence, copy_silence }, /* 01: */
443     {  1, copy_u8     , copy_u8      }, /* 02: */
444     {  1, copy_s8     , copy_s8      }, /* 03: */
445     {  2, copy_u16_be , copy_u16_le  }, /* 04: */
446     {  2, copy_s16_be , copy_s16_le  }, /* 05: */
447     {  3, copy_u24_be , copy_u24_le  }, /* 06: */
448     {  3, copy_s24_be , copy_s24_le  }, /* 07: */
449     {  4, copy_u32_be , copy_u32_le  }, /* 08: */
450     {  4, copy_s32_be , copy_s32_le  }, /* 09: */
451     { -1, NULL        , NULL         }, /* 0A: */
452     { -1, NULL        , NULL         }, /* 0B: */
453     { -1, NULL        , NULL         }, /* 0C: */
454     { -1, NULL        , NULL         }, /* 0D: */
455     { -1, NULL        , NULL         }, /* 0E: */
456     { -1, NULL        , NULL         }, /* 0F: */
457     { -1, NULL        , NULL         }, /* 10: */
458     { -1, NULL        , NULL         }, /* 11: */
459     { -1, NULL        , NULL         }, /* 12: */
460     { -1, NULL        , NULL         }, /* 13: */
461     {  4, copy_f32_be , copy_f32_le  }, /* 14: */
462     { -1, NULL        , NULL         }, /* 15: */
463     { -1, NULL        , NULL         }, /* 16: */
464     { -1, NULL        , NULL         }, /* 17: */
465     {  8, copy_f64_be , copy_f64_le  }, /* 18: */
466     { -1, NULL        , NULL         }, /* 19: */
467     { 10, copy_f80_be , copy_f80_le  }, /* 1A: */
468     { -1, NULL        , NULL         }, /* 1B: */
469     { 12, copy_f80_be , copy_f80_le  }, /* 1C: */
470     { -1, NULL        , NULL         }, /* 1D: */
471     { -1, NULL        , NULL         }, /* 1E: */
472     { -1, NULL        , NULL         }, /* 1F: */
473     { 16, copy_f80_be , copy_f80_le  }, /* 20: */
474 
475     { sizeof(short)      , NULL        , NULL         }, /* 21: */
476     { sizeof(int)        , NULL        , NULL         }, /* 22: */
477     { sizeof(long)       , NULL        , NULL         }, /* 23: */
478     { sizeof(float)      , copy_f32_be , copy_f32_le  }, /* 24: */
479     { sizeof(double)     , copy_f64_be , copy_f64_le  }, /* 25: */
480     { sizeof(long double), NULL        , NULL         }, /* 26: */
481 
482     { -1, NULL        , NULL         }, /* 27: */
483     { -1, NULL        , NULL         }, /* 28: */
484     { -1, NULL        , NULL         }, /* 29: */
485     { -1, NULL        , NULL         }, /* 2A: */
486     { -1, NULL        , NULL         }, /* 2B: */
487     { -1, NULL        , NULL         }, /* 2C: */
488     { -1, NULL        , NULL         }, /* 2D: */
489     { -1, NULL        , NULL         }, /* 2E: */
490     { -1, NULL        , NULL         }, /* 2F: */
491     { -1, NULL        , NULL         }, /* 30: */
492     {  1, copy_alaw   , copy_alaw    }, /* 31: */
493     {  1, copy_ulaw   , copy_ulaw    }, /* 32: */
494     { -1, NULL        , NULL         }, /* 33: */
495     { -1, NULL        , NULL         }, /* 34: */
496     { -1, NULL        , NULL         }, /* 35: */
497     { -1, NULL        , NULL         }, /* 36: */
498     { -1, NULL        , NULL         }, /* 37: */
499     { -1, NULL        , NULL         }, /* 38: */
500     { -1, NULL        , NULL         }, /* 39: */
501     { -1, NULL        , NULL         }, /* 3A: */
502     { -1, NULL        , NULL         }, /* 3B: */
503     { -1, NULL        , NULL         }, /* 3C: */
504     { -1, NULL        , NULL         }, /* 3D: */
505     { -1, NULL        , NULL         }, /* 3E: */
506     { -1, NULL        , NULL         }, /* 3F: */
507 };
508 
509 /*
510  *  Selects the right demultiplexer from the attribute field
511  *  (currently endian and type information is used, rest is done in
512  *  lame_encode_pcm)
513  */
514 
select_demux(OUT uint32_t mode,IN demux_t * retf,IN ssize_t * size)515 static inline int  select_demux ( OUT uint32_t mode, IN demux_t* retf, IN ssize_t* size )
516 {
517     int                  big    = mode >> 24;
518     const demux_info_t*  tabptr = demux_info + ((mode >> 16) & 0x3F);
519 
520     FN("select_demux");
521 
522 #ifdef WORDS_BIGENDIAN
523     /* big endian */
524     big  = (big >> 1) ^ big ^ 1;        // 0=big, 1=little, 2=little, 3=big
525 #else
526     /* little endian */                 // 0=little, 1=big, 2=little, 3=big
527 #endif
528 
529     *size = tabptr -> size;
530     *retf = big & 1  ?  tabptr -> demultiplexer_be
531                      :  tabptr -> demultiplexer_le;
532     return *retf != NULL  ?  0  :  -1;
533 }
534 
535 /*}}}*/
536 /*{{{ amplify/resample stuff          */
537 
538 /*
539  *  routine to feed EXACTLY one frame (lame->frame_size) worth of data to the
540  *  encoding engine. All buffering, resampling, etc, handled by calling program.
541  */
542 
internal_lame_encoding_pcm(INOUT lame_t * const lame,INOUT octetstream_t * const os,OUT sample_t * const * const data,OUT size_t len)543 static inline int  internal_lame_encoding_pcm (
544         INOUT lame_t*           const lame,
545         INOUT octetstream_t*    const os,
546         OUT   sample_t* const * const data,
547         OUT   size_t                  len )
548 {
549     size_t           i;
550     double           ampl;
551     double           dampl;
552     int              ampl_on;
553     size_t           ch;
554     size_t           mf_needed;
555     int              ret;
556     size_t           n_in;
557     size_t           n_out;
558     size_t           remaining = len;
559     sample_t*        mfbuf [MAX_CHANNELS];
560     const sample_t*  pdata [MAX_CHANNELS];
561 
562     FN("internal_lame_encoding_pcm");
563 
564     /// this should be moved to lame_init_params();
565     // some sanity checks
566 #if ENCDELAY < MDCTDELAY
567 # error ENCDELAY is less than MDCTDELAY, see <encoder.h>
568 #endif
569 #if FFTOFFSET > BLKSIZE
570 # error FFTOFFSET is greater than BLKSIZE, see <encoder.h>
571 #endif
572 
573     mf_needed = BLKSIZE + lame->frame_size - FFTOFFSET;           // amount needed for FFT
574     mf_needed = MAX ( mf_needed, 286 + 576 * (1+lame->mode_gr) ); // amount needed for MDCT/filterbank
575     assert ( mf_needed <= MFSIZE );
576     ///
577 
578     os -> length = 0;
579     pdata [0]    = data [0];
580     pdata [1]    = data [1];
581     mfbuf [0]    = lame->mfbuf [0];
582     mfbuf [1]    = lame->mfbuf [1];
583 
584     ampl    = lame->last_ampl;
585     dampl   = 0.;
586     if ( lame->ampl != ampl ) {
587         if (remaining <= 1) {           // constant amplification
588             ampl_on = 1;
589         } else {                        // fading
590             ampl_on = 2;
591             dampl   = (lame->ampl - ampl) / remaining * lame->sampfreq_out / lame->sampfreq_in;
592         }
593         lame->last_ampl = lame->ampl;
594     } else if ( lame->ampl != 1. ) {    // constant amplification
595         ampl_on = 1;
596     } else {                            // no manipulation (fastest)
597         ampl_on = 0;
598     }
599 
600     while ( (int) remaining > 0 ) {
601 
602         FN("internal_lame_encoding_pcm 3");
603 
604         /* copy in new samples into mfbuf, with resampling if necessary */
605         if ( lame->resample_in != NULL )  {
606 
607             FN("internal_lame_encoding_pcm 10");
608 
609             for ( ch = 0; ch < lame->channels_out; ch++ ) {
610                 n_in  = remaining;
611                 n_out = lame->frame_size;
612 
613                 FN("internal_lame_encoding_pcm 12");
614 
615                 // resample filter virtually should have no delay !
616                 ret   = resample_buffer (
617                             lame->resample_in,
618                             mfbuf [ch] + lame->mf_size, &n_out,
619                             pdata [ch]                , &n_in,
620                             ch );
621                 if ( ret < 0 )
622                     return ret;
623                 pdata [ch] += n_in;
624             }
625 
626             FN("internal_lame_encoding_pcm 13");
627 
628         }
629         else {
630 
631             FN("internal_lame_encoding_pcm 14");
632 
633             n_in = n_out = MIN ( lame->frame_size, remaining );
634 
635             FN("internal_lame_encoding_pcm 15");
636 
637             for ( ch = 0; ch < lame->channels_out; ch++ ) {
638                 memcpy (  mfbuf [ch] + lame->mf_size,
639                           pdata [ch],
640                           n_out * sizeof (**mfbuf) );
641                 pdata [ch] += n_in;
642             }
643 
644             FN("internal_lame_encoding_pcm 16");
645 
646         }
647 
648         FN("internal_lame_encoding_pcm 4");
649 
650         switch ( ampl_on ) {
651         case 0:
652             break;
653         case 1:
654             for ( ch = 0; ch < lame->channels_out; ch++ )
655                 for ( i = 0; i < n_out; i++ )
656                     mfbuf [ch] [lame->mf_size + i] *= ampl;
657             break;
658         case 2:
659             for ( i = 0; i < n_out; i++, ampl += dampl )
660                 for ( ch = 0; ch < lame->channels_out; ch++ )
661                     mfbuf [ch] [lame->mf_size + i] *= ampl;
662             break;
663         default:
664             assert (0);
665             break;
666         }
667 
668         FN("internal_lame_encoding_pcm 4");
669 
670         fprintf ( stderr, "n_in=%d, n_out=%d, remaining=%d, remaining=%d\n", n_in, n_out, remaining, remaining-n_in );
671 
672         remaining                  -= n_in;
673         lame->mf_size              += n_out;
674         lame->mf_samples_to_encode += n_out;
675 
676         // encode ONE frame if enough data available
677         if ( lame->mf_size >= mf_needed ) {
678             // Enlarge octetstream buffer if (possibly) needed by (25% + 16K)
679             if ( os->size < 16384 + os->length )
680                 octetstream_resize ( os, os->size + os->size/4 + 16384 );
681             // Encode one frame
682             ret = lame_encode_frame ( lame, mfbuf,
683                                       os->data + os->length, os->size - os->length );
684 
685             if (ret < 0)
686                 return ret;
687             os->length += ret;
688 
689             FN("internal_lame_encoding_pcm 5");
690 
691             // shift out old samples
692             lame->mf_size              -= lame->frame_size;
693             lame->mf_samples_to_encode -= lame->frame_size;
694             for ( ch = 0; ch < lame->channels_out; ch++ )
695                 memmove ( mfbuf [ch] + 0, mfbuf [ch] + lame->frame_size, lame->mf_size * sizeof (**mfbuf) );
696 
697             FN("internal_lame_encoding_pcm 6");
698 
699         }
700     }
701     assert (remaining == 0);
702     return 0;
703 }
704 
705 /*}}}*/
706 /*{{{ demultiplexing stuff            */
707 
average(sample_t * dst,const sample_t * src1,const sample_t * src2,size_t len)708 static inline void  average ( sample_t* dst, const sample_t* src1, const sample_t* src2, size_t len )
709 {
710     FN("average");
711 
712     while (len--)
713         *dst++ = (*src1++ + *src2++) * 0.5;
714 }
715 
lame_encode_pcm(lame_t * const lame,octetstream_t * os,const void * pcm,size_t len,uint32_t flags)716 int  lame_encode_pcm (
717         lame_t* const   lame,
718         octetstream_t*  os,
719         const void*     pcm,
720         size_t          len,
721         uint32_t        flags )
722 {
723     sample_t*           data [2];
724     demux_t             retf;
725     ssize_t             size;
726     const uint8_t*      p;
727     const uint8_t* const*  q;
728     int                 ret;
729     size_t              channels_in = flags & 0xFFFF;
730 
731     FN("lame_encode_pcm");
732 
733     if (len == 0)
734         return 0;
735 
736     if (select_demux ( flags, &retf, &size ) < 0)
737         return -1;
738 
739     data[0] = (sample_t*) calloc (sizeof(sample_t), len);
740     data[1] = (sample_t*) calloc (sizeof(sample_t), len);
741 
742     switch (flags >> 28) {
743     case LAME_INTERLEAVED >> 28:
744         p = (const uint8_t*) pcm;
745         switch ( lame->channels_out ) {
746         case 1:
747             switch ( channels_in ) {
748             case 1:
749                 retf ( data[0], p+0*size, 1*size, len );
750                 break;
751             case 2:
752                 retf ( data[0], p+0*size, 2*size, len );
753                 retf ( data[1], p+1*size, 2*size, len );
754                 average ( data[0], data[0], data[1], len );
755                 memset ( data[1], 0, sizeof(sample_t)*len );
756                 break;
757             default:
758                 return -1;
759                 break;
760             }
761             break;
762         case 2:
763             switch ( channels_in ) {
764             case 1:
765                 retf ( data[0], p+0*size, 1*size, len );
766                 memcpy ( data[1], data[0], sizeof(sample_t)*len );
767                 break;
768             case 2:
769                 retf ( data[0], p+0*size, 2*size, len );
770                 retf ( data[1], p+1*size, 2*size, len );
771                 break;
772             default:
773                 return -1;
774                 break;
775             }
776             break;
777         default:
778             return -1;
779         }
780         break;
781 
782     case LAME_CHAINED >> 28:
783         p = (const uint8_t*) pcm;
784         switch ( lame->channels_out ) {
785         case 1:
786             switch ( channels_in ) {
787             case 1:
788                 retf ( data[0], p+0*size*len, size, len );
789                 break;
790             case 2:
791                 retf ( data[0], p+0*size*len, size, len );
792                 retf ( data[1], p+1*size*len, size, len );
793                 average ( data[0], data[0], data[1], len );
794                 memset ( data[1], 0, sizeof(sample_t)*len );
795                 break;
796             default:
797                 return -1;
798                 break;
799             }
800             break;
801         case 2:
802             switch ( channels_in ) {
803             case 1:
804                 retf ( data[0], p+0*size*len, size, len );
805                 memcpy ( data[1], data[0], sizeof(sample_t)*len );
806                 break;
807             case 2:
808                 retf ( data[0], p+0*size*len, size, len );
809                 retf ( data[1], p+1*size*len, size, len );
810                 break;
811             default:
812                 return -1;
813                 break;
814             }
815             break;
816         default:
817             return -1;
818         }
819         break;
820 
821     case LAME_INDIRECT >> 28:
822         q = (const uint8_t* const*) pcm;
823         switch ( lame->channels_out ) {
824         case 1:
825             switch ( channels_in ) {
826             case 1:
827                 retf ( data[0], q[0], size, len );
828                 break;
829             case 2:
830                 retf ( data[0], q[0], size, len );
831                 retf ( data[1], q[1], size, len );
832                 average ( data[0], data[0], data[1], len );
833                 memset ( data[1], 0, sizeof(sample_t)*len );
834                 break;
835             default:
836                 return -1;
837                 break;
838             }
839             break;
840         case 2:
841             switch ( channels_in ) {
842             case 1:
843                 retf ( data[0], q[0], size, len );
844                 memcpy ( data[1], data[0], sizeof(sample_t)*len );
845                 break;
846             case 2:
847                 retf ( data[0], q[0], size, len );
848                 retf ( data[1], q[1], size, len );
849                 break;
850             default:
851                 return -1;
852                 break;
853             }
854             break;
855         default:
856             return -1;
857         }
858         break;
859 
860     default:
861         return -1;
862     }
863 
864     ret = internal_lame_encoding_pcm ( lame, os, (sample_t const* const*) data, len );
865 
866     free ( data[0] );
867     free ( data[1] );
868     return ret;
869 }
870 
871 /*}}}*/
872 /*{{{ lame_encode_pcm_flush           */
873 
lame_encode_pcm_flush(lame_t * const lame,octetstream_t * const os)874 int  lame_encode_pcm_flush (
875         lame_t*        const  lame,
876         octetstream_t* const  os )
877 {
878     int  ret;
879     int  ret_cb;
880 
881     FN("lame_encode_pcm_flush");
882 
883     ret = lame_encode_pcm ( lame, os, NULL,
884                             lame->mf_samples_to_encode + lame->frame_size,
885 			    LAME_INTERLEAVED | LAME_SILENCE | 1 );
886 
887 
888     // ugly, not encapsulated
889     flush_bitstream (lame -> global_flags);       // mp3 related stuff. bit buffer might still contain some mp3 data
890     id3tag_write_v1 (lame -> global_flags);       // write a ID3 tag to the bitstream
891     ret_cb = copy_buffer (os->data, os->size - os->length, &lame->bs);
892     if (ret_cb < 0)
893         return ret_cb;
894     os->length += ret_cb;
895 
896 
897     return ret;
898 }
899 
900 /*
901  *  Data flow
902  *  ~~~~~~~~~
903  *  lame_encode_buffer (public)
904  *  lame_encode_buffer_interleaved (public)
905  *      - do some little preparations and calls lame_encode_pcm
906  *  lame_encode_pcm (public)
907  *      - demultiplexing
908  *      - type converting
909  *      - endian handling
910  *      - alignment handling
911  *      - channel number converting (currently averaging or duplicating)
912  *  internal_lame_encoding_pcm
913  *      - resampling
914  *      - level attenuator/amplifier, fader
915  *      - divide input into frames
916  *      - merge encoded data to one stream
917  *      - does a lot of things which should be done in the setup
918  *  lame_encode_frame
919  *      - calls the right next stage encoder
920  *      - counts the number of encoded frames
921  *      - analyzer call back
922  *  lame_encode_???_frame
923  *      - out of the scope of this file
924  */
925 
926 /*}}}*/
927 /*{{{ Legacy stuff                    */
928 
pointer2lame(void * const handle)929 static lame_t*  pointer2lame ( void* const handle )
930 {
931     lame_t*             lame = (lame_t*)handle;
932     lame_global_flags*  gfp  = (lame_global_flags*)handle;
933 
934     FN("pointer2lame");
935 
936     if ( lame == NULL )
937         return NULL;
938     if ( lame->Class_ID == LAME_ID )
939         return lame;
940 
941     if ( gfp->num_channels == 1  ||  gfp->num_channels == 2 ) {
942         lame = (lame_t*) (gfp->internal_flags);
943 
944         if ( lame == NULL )
945             return NULL;
946         if ( lame->Class_ID == LAME_ID )
947             return lame;
948     }
949     return NULL;
950 }
951 
952 
lame_encode_buffer(void * const gfp,const int16_t buffer_l[],const int16_t buffer_r[],const size_t nsamples,void * const mp3buf,const size_t mp3buf_size)953 int  lame_encode_buffer (
954         void* const    gfp,
955         const int16_t  buffer_l [],
956         const int16_t  buffer_r [],
957         const size_t   nsamples,
958         void* const    mp3buf,
959         const size_t   mp3buf_size )
960 {
961     const int16_t*  pcm [2];
962     octetstream_t*  os;
963     int             ret;
964     lame_t*         lame;
965 
966     FN("lame_encode_buffer");
967 
968     lame    = pointer2lame (gfp);
969     os      = octetstream_open (mp3buf_size);
970     pcm [0] = buffer_l;
971     pcm [1] = buffer_r;
972 
973     ret     = lame_encode_pcm ( lame, os, pcm, nsamples,
974                                 LAME_INDIRECT | LAME_NATIVE_ENDIAN | LAME_INT16 |
975                                 lame->channels_in );
976 
977     memcpy ( mp3buf, os->data, os->length );
978     if (ret == 0) ret = os->length;
979     octetstream_close (os);
980     return ret;
981 }
982 
983 
lame_encode_buffer_interleaved(void * const gfp,const int16_t buffer[],size_t nsamples,void * const mp3buf,const size_t mp3buf_size)984 int  lame_encode_buffer_interleaved (
985         void* const    gfp,
986         const int16_t  buffer [],
987         size_t         nsamples,
988         void* const    mp3buf,
989         const size_t   mp3buf_size )
990 {
991     octetstream_t*  os;
992     int             ret;
993     lame_t*         lame;
994 
995     FN("lame_encode_buffer_interleaved");
996 
997     lame    = pointer2lame (gfp);
998     os      = octetstream_open (mp3buf_size);
999 
1000     ret     = lame_encode_pcm ( lame, os, buffer, nsamples,
1001                                 LAME_INTERLEAVED | LAME_NATIVE_ENDIAN | LAME_INT16 |
1002                                 lame->channels_in );
1003 
1004     memcpy ( mp3buf, os->data, os->length );
1005     if (ret == 0) ret = os->length;
1006     octetstream_close (os);
1007     return ret;
1008 }
1009 
1010 
lame_encode_flush(void * const gfp,void * const mp3buf,const size_t mp3buf_size)1011 int  lame_encode_flush (
1012         void* const    gfp,
1013         void* const    mp3buf,
1014         const size_t   mp3buf_size )
1015 {
1016     octetstream_t*  os;
1017     int             ret;
1018     lame_t*         lame;
1019 
1020     FN("lame_encode_flush");
1021 
1022     lame    = pointer2lame (gfp);
1023     os      = octetstream_open (mp3buf_size);
1024 
1025     ret     = lame_encode_pcm_flush ( lame, os );
1026 
1027     memcpy ( mp3buf, os->data, os->length );
1028     if (ret == 0) ret = os->length;
1029     octetstream_close (os);
1030     return ret;
1031 }
1032 
1033 /*}}}*/
1034 /*{{{ sin/cos/sinc/rounding functions */
1035 
round_nearest(long double x)1036 static inline long  round_nearest ( long double x )
1037 {
1038     if ( x >= 0. )
1039         return (long)(x+0.5);
1040     else
1041         return (long)(x-0.5);
1042 }
1043 
1044 
round_down(long double x)1045 static inline long  round_down ( long double x )
1046 {
1047     if ( x >= 0. )
1048         return +(long)(+x);
1049     else
1050         return -(long)(-x);
1051 }
1052 
1053 
sinpi(long double x)1054 static inline long double  sinpi ( long double x )
1055 {
1056     x -= round_down (x) & ~1;
1057 
1058     switch ( (int)(4.*x) ) {
1059     default: assert (0);
1060     case  0: return +SIN ( M_PIl * (0.0+x) );
1061     case  1: return +COS ( M_PIl * (0.5-x) );
1062     case  2: return +COS ( M_PIl * (x-0.5) );
1063     case  3: return +SIN ( M_PIl * (1.0-x) );
1064     case  4: return -SIN ( M_PIl * (x-1.0) );
1065     case  5: return -COS ( M_PIl * (1.5-x) );
1066     case  6: return -COS ( M_PIl * (x-1.5) );
1067     case  7: return -SIN ( M_PIl * (2.0-x) );
1068     }
1069 }
1070 
1071 
cospi(long double x)1072 static inline long double  cospi ( long double x )
1073 {
1074     x -= round_down (x) & ~1;
1075 
1076     switch ( (int)(4.*x) ) {
1077     default: assert (0);
1078     case  0: return +COS ( M_PIl * (x-0.0) );
1079     case  1: return +SIN ( M_PIl * (0.5-x) );
1080     case  2: return -SIN ( M_PIl * (x-0.5) );
1081     case  3: return -COS ( M_PIl * (1.0-x) );
1082     case  4: return -COS ( M_PIl * (x-1.0) );
1083     case  5: return -SIN ( M_PIl * (1.5-x) );
1084     case  6: return +SIN ( M_PIl * (x-1.5) );
1085     case  7: return +COS ( M_PIl * (2.0-x) );
1086     }
1087 }
1088 
1089 
sinc(long double x)1090 static inline long double  sinc ( long double x )
1091 {
1092     if ( x == 0. )
1093         return 1.;
1094 
1095     if ( x <  0. )
1096         x = -x;
1097 
1098     return sinpi ( x ) / ( M_PIl * x );
1099 }
1100 
1101 /*}}}*/
1102 /*{{{ some window functions           */
1103 
hanning(double x)1104 static inline  double  hanning ( double x )
1105 {
1106     if ( fabs (x) >= 1 )
1107         return 0.;
1108 
1109     x = cospi (0.5 * x);
1110     return x * x;
1111 }
1112 
hamming(double x)1113 static inline  double  hamming ( double x )
1114 {
1115     if ( fabs (x) >= 1 )
1116         return 0.;
1117 
1118     x = cospi (x);
1119     return 0.54 + 0.46 * x;
1120 }
1121 
blackman(double x)1122 static inline  double  blackman ( double x )
1123 {
1124     if ( fabs (x) >= 1 )
1125         return 0.;
1126 
1127     x = cospi (x);
1128     return (0.16 * x + 0.50) * x + 0.34;        // using addition theorem of arc functions
1129 }
1130 
blackman1(double x)1131 static inline  double  blackman1 ( double x )
1132 {
1133     if ( fabs (x) >= 1 )
1134         return 0.;
1135 
1136     x += 1.;
1137     return 0.42 - 0.50*cospi(x) + 0.08*cospi(2*x);
1138 }
1139 
blackman2(double x)1140 static inline  double  blackman2 ( double x )
1141 {
1142     if ( fabs (x) >= 1 )
1143         return 0.;
1144 
1145     x += 1.;
1146     return 0.375 - 0.50*cospi(x) + 0.125*cospi(2*x);
1147 }
1148 
blackmanharris_nuttall(double x)1149 static inline  double  blackmanharris_nuttall ( double x )
1150 {
1151     if ( fabs (x) >= 1 )
1152         return 0.;
1153 
1154     x += 1.;
1155     return (10 - 15*cospi (x) + 6*cospi (2*x) - cospi (3*x) ) * (1./32);
1156 }
1157 
blackmanharris_min4(double x)1158 static inline  double  blackmanharris_min4 ( double x )
1159 {
1160     if ( fabs (x) >= 1 )
1161         return 0.;
1162 
1163     x += 1.;
1164     return 0.355768 - 0.487396*cospi (x) + 0.144232*cospi (2*x) - 0.012604*cospi (3*x);
1165 }
1166 
1167 /*
1168  * Blackman-Harris windows, which have the general equation:
1169  *
1170  *          w(n) = a0 - a1 cos(2pn/N) + a2 cos(4pn/N) - a3 cos(6pn/N).
1171  *
1172  * When Nuttall set the ak coefficients to the "minimum 4-term" values of a0 = 0.355768, a1 =
1173  * 0.487396, a2 = 0.144232, and a3 = 0.012604, he obtained the window and the frequency
1174  * response shown in Figure 8. The minimum-4-term window loses some frequency
1175  * resolution because it has a wide main lobe, but the first side lobe drops to -93 db.
1176  *
1177  * If you require fast side-lobe rolloff in an application, consider Nuttall's Equation (30). It
1178  * applies the following coefficients to the general equation above: a0 = 10/32, a1 = 15/32, a2 =
1179  * 6/32, and a3 = 1/32. Figure 8 also includes the window and frequency response for Nuttall's
1180  * Equation (30). This window has a first side-lobe level of -61 dB and a significant side-lobe
1181  * rolloff of -42 dB/octave.
1182  *
1183  * The first investigators of windowing functions determined that window functions should
1184  * have zero values at their boundaries and so should their successive derivatives. If a
1185  * window's kth derivative is zero at the boundaries, the peaks of the window's side lobes will
1186  * decay at a rate of 6(k+2) dB/octave. T&MW.
1187  */
1188 
1189 /*}}}*/
1190 /*{{{ scalar stuff                    */
1191 
1192 
1193 /**********************************************************
1194  *                                                        *
1195  *           Functions taken from scalar.nas              *
1196  *                                                        *
1197  **********************************************************/
1198 
1199 /*
1200  * The scalarxx versions with xx=04,08,12,16,20,24 are fixed size for xx element scalars
1201  * The scalar1n version is suitable for every non negative length
1202  * The scalar4n version is only suitable for positive lengths which are a multiple of 4
1203  *
1204  * The following are equivalent:
1205  *    scalar12 (p, q);
1206  *    scalar4n (p, q, 3);
1207  *    scalar1n (p, q, 12);
1208  */
1209 
scalar04_float32(const sample_t * p,const sample_t * q)1210 float_t  scalar04_float32 ( const sample_t* p, const sample_t* q )
1211 {
1212     return  p[0]*q[0] + p[1]*q[1] + p[2]*q[2] + p[3]*q[3];
1213 }
1214 
scalar08_float32(const sample_t * p,const sample_t * q)1215 float_t  scalar08_float32 ( const sample_t* p, const sample_t* q )
1216 {
1217     return  p[0]*q[0] + p[1]*q[1] + p[2]*q[2] + p[3]*q[3]
1218           + p[4]*q[4] + p[5]*q[5] + p[6]*q[6] + p[7]*q[7];
1219 }
1220 
scalar12_float32(const sample_t * p,const sample_t * q)1221 float_t  scalar12_float32 ( const sample_t* p, const sample_t* q )
1222 {
1223     return  p[0]*q[0] + p[1]*q[1] + p[ 2]*q[ 2] + p[ 3]*q[ 3]
1224           + p[4]*q[4] + p[5]*q[5] + p[ 6]*q[ 6] + p[ 7]*q[ 7]
1225           + p[8]*q[8] + p[9]*q[9] + p[10]*q[10] + p[11]*q[11];
1226 }
1227 
scalar16_float32(const sample_t * p,const sample_t * q)1228 float_t  scalar16_float32 ( const sample_t* p, const sample_t* q )
1229 {
1230     return  p[ 0]*q[ 0] + p[ 1]*q[ 1] + p[ 2]*q[ 2] + p[ 3]*q[ 3]
1231           + p[ 4]*q[ 4] + p[ 5]*q[ 5] + p[ 6]*q[ 6] + p[ 7]*q[ 7]
1232           + p[ 8]*q[ 8] + p[ 9]*q[ 9] + p[10]*q[10] + p[11]*q[11]
1233           + p[12]*q[12] + p[13]*q[13] + p[14]*q[14] + p[15]*q[15];
1234 }
1235 
scalar20_float32(const sample_t * p,const sample_t * q)1236 float_t  scalar20_float32 ( const sample_t* p, const sample_t* q )
1237 {
1238     return  p[ 0]*q[ 0] + p[ 1]*q[ 1] + p[ 2]*q[ 2] + p[ 3]*q[ 3]
1239           + p[ 4]*q[ 4] + p[ 5]*q[ 5] + p[ 6]*q[ 6] + p[ 7]*q[ 7]
1240           + p[ 8]*q[ 8] + p[ 9]*q[ 9] + p[10]*q[10] + p[11]*q[11]
1241           + p[12]*q[12] + p[13]*q[13] + p[14]*q[14] + p[15]*q[15]
1242           + p[16]*q[16] + p[17]*q[17] + p[18]*q[18] + p[19]*q[19];
1243 }
1244 
scalar24_float32(const sample_t * p,const sample_t * q)1245 float_t  scalar24_float32 ( const sample_t* p, const sample_t* q )
1246 {
1247     return  p[ 0]*q[ 0] + p[ 1]*q[ 1] + p[ 2]*q[ 2] + p[ 3]*q[ 3]
1248           + p[ 4]*q[ 4] + p[ 5]*q[ 5] + p[ 6]*q[ 6] + p[ 7]*q[ 7]
1249           + p[ 8]*q[ 8] + p[ 9]*q[ 9] + p[10]*q[10] + p[11]*q[11]
1250           + p[12]*q[12] + p[13]*q[13] + p[14]*q[14] + p[15]*q[15]
1251           + p[16]*q[16] + p[17]*q[17] + p[18]*q[18] + p[19]*q[19]
1252           + p[20]*q[20] + p[21]*q[21] + p[22]*q[22] + p[23]*q[23];
1253 }
1254 
scalar32_float32(const sample_t * p,const sample_t * q)1255 float_t  scalar32_float32 ( const sample_t* p, const sample_t* q )
1256 {
1257     return  p[ 0]*q[ 0] + p[ 1]*q[ 1] + p[ 2]*q[ 2] + p[ 3]*q[ 3]
1258           + p[ 4]*q[ 4] + p[ 5]*q[ 5] + p[ 6]*q[ 6] + p[ 7]*q[ 7]
1259           + p[ 8]*q[ 8] + p[ 9]*q[ 9] + p[10]*q[10] + p[11]*q[11]
1260           + p[12]*q[12] + p[13]*q[13] + p[14]*q[14] + p[15]*q[15]
1261           + p[16]*q[16] + p[17]*q[17] + p[18]*q[18] + p[19]*q[19]
1262           + p[20]*q[20] + p[21]*q[21] + p[22]*q[22] + p[23]*q[23]
1263           + p[24]*q[24] + p[25]*q[25] + p[26]*q[26] + p[27]*q[27]
1264           + p[28]*q[28] + p[29]*q[29] + p[30]*q[30] + p[31]*q[31];
1265 }
1266 
scalar4n_float32(const sample_t * p,const sample_t * q,size_t len)1267 float_t  scalar4n_float32 ( const sample_t* p, const sample_t* q, size_t len )
1268 {
1269     double sum = p[0]*q[0] + p[1]*q[1] + p[2]*q[2] + p[3]*q[3];
1270 
1271     while (--len) {
1272         p   += 4;
1273         q   += 4;
1274         sum += p[0]*q[0] + p[1]*q[1] + p[2]*q[2] + p[3]*q[3];
1275     }
1276     return sum;
1277 }
1278 
scalar1n_float32(const sample_t * p,const sample_t * q,size_t len)1279 float_t  scalar1n_float32 ( const sample_t* p, const sample_t* q, size_t len )
1280 {
1281     float_t sum;
1282 
1283     switch (len & 3) {
1284     case 0:
1285         sum  = 0.;
1286         break;
1287     case 1:
1288         sum  = p[0]*q[0];
1289         p   += 1;
1290         q   += 1;
1291         break;
1292     case 2:
1293         sum  = p[0]*q[0] + p[1]*q[1];
1294         p   += 2;
1295         q   += 2;
1296         break;
1297     default:
1298         sum  = p[0]*q[0] + p[1]*q[1] + p[2]*q[2];
1299         p   += 3;
1300         q   += 3;
1301         break;
1302     }
1303 
1304     for ( len >>= 2; len--; ) {
1305         sum += p[0]*q[0] + p[1]*q[1] + p[2]*q[2] + p[3]*q[3];
1306         p   += 4;
1307         q   += 4;
1308     }
1309     return sum;
1310 }
1311 
1312 
1313 scalar_t   scalar04 = scalar04_float32;
1314 scalar_t   scalar08 = scalar08_float32;
1315 scalar_t   scalar12 = scalar12_float32;
1316 scalar_t   scalar16 = scalar16_float32;
1317 scalar_t   scalar20 = scalar20_float32;
1318 scalar_t   scalar24 = scalar24_float32;
1319 scalar_t   scalar32 = scalar32_float32;
1320 scalarn_t  scalar4n = scalar4n_float32;
1321 scalarn_t  scalar1n = scalar1n_float32;
1322 
1323 
init_scalar_functions(OUT lame_t * const lame)1324 void  init_scalar_functions ( OUT lame_t* const lame )
1325 {
1326     if ( lame -> CPU_features.i387 ) {
1327         scalar04 = scalar04_float32_i387;
1328         scalar08 = scalar08_float32_i387;
1329         scalar12 = scalar12_float32_i387;
1330         scalar16 = scalar16_float32_i387;
1331         scalar20 = scalar20_float32_i387;
1332         scalar24 = scalar24_float32_i387;
1333         scalar32 = scalar32_float32_i387;
1334         scalar4n = scalar4n_float32_i387;
1335         scalar1n = scalar1n_float32_i387;
1336     }
1337 
1338     if ( lame -> CPU_features.AMD_3DNow ) {
1339         scalar04 = scalar04_float32_3DNow;
1340         scalar08 = scalar08_float32_3DNow;
1341         scalar12 = scalar12_float32_3DNow;
1342         scalar16 = scalar16_float32_3DNow;
1343         scalar20 = scalar20_float32_3DNow;
1344         scalar24 = scalar24_float32_3DNow;
1345         scalar32 = scalar32_float32_3DNow;
1346         scalar4n = scalar4n_float32_3DNow;
1347         scalar1n = scalar1n_float32_3DNow;
1348     }
1349 
1350     if ( lame -> CPU_features.SIMD ) {
1351         scalar04 = scalar04_float32_SIMD;
1352         scalar08 = scalar08_float32_SIMD;
1353         scalar12 = scalar12_float32_SIMD;
1354         scalar16 = scalar16_float32_SIMD;
1355         scalar20 = scalar20_float32_SIMD;
1356         scalar24 = scalar24_float32_SIMD;
1357         scalar32 = scalar32_float32_SIMD;
1358         scalar4n = scalar4n_float32_SIMD;
1359         scalar1n = scalar1n_float32_SIMD;
1360     }
1361 }
1362 
1363 /*}}}*/
1364 /*{{{ factorize                       */
1365 
1366 /*
1367  *  Tries to find the best integral ratio for two sampling frequencies
1368  *  It searches the best ratio a:b with b <= MAX_TABLES
1369  */
1370 
factorize(OUT long double f1,OUT long double f2,IN int * const x1,IN int * const x2)1371 static double  factorize (
1372         OUT long double  f1,
1373 	OUT long double  f2,
1374 	IN  int* const   x1,
1375 	IN  int* const   x2 )
1376 {
1377     unsigned     i;
1378     long         ltmp;
1379     long double  ftmp;
1380     double       minerror = 1.;
1381     double       abserror = 1.;
1382 
1383     assert ( f1 > 0. );
1384     assert ( f2 > 0. );
1385     assert ( x1 != NULL );
1386     assert ( x2 != NULL );
1387 
1388     for ( i = 1; i <= MAX_TABLES; i++ ) {
1389         ftmp  = f2 * i / f1;
1390         ltmp  = (long) ( ftmp + 0.5 );
1391         if ( fabs ( (ltmp-ftmp)/ftmp ) < minerror ) {
1392             *x1      = i;
1393             *x2      = (int)ltmp;
1394             abserror = (ltmp-ftmp) / ftmp;
1395             minerror = 0.9999847412109375 * fabs (abserror);
1396         }
1397     }
1398     return abserror;
1399 }
1400 
1401 /*
1402  *  Some programmers and some file format only supporting integral sampling frequencies
1403  *  This patch tries to find the right sampling frequency for some know rounding victims.
1404  *  Better is to support directly 64 or 80 bit FPU precision.
1405  */
1406 
unround_samplefrequency(OUT long double freq)1407 long double  unround_samplefrequency ( OUT long double freq )
1408 {
1409     if ( freq != (unsigned short)freq )
1410         return freq;
1411 
1412     switch ( (unsigned short)freq ) {
1413     case 48662: case 48663: return 1411200/ 29.L;  // SB 128 "48 kHz"         48662.06896551724137931034
1414     case 44055: case 44056: return 2863636/ 65.L;  // NTSC PCM/LD PCM         44055.93846153846153846153
1415     case 32072: case 32073: return 1411200/ 44.L;  // SB 128 "32 kHz"         32072.72727272727272727272
1416     case 31468: case 31469: return 2863636/ 91.L;  // NTSC Hi8 Digital Audio  31468.52747252747252747252
1417     case 23918: case 23919: return 1411200/ 59.L;  // SB 128 "24 kHz"         23918.64406779661016949152
1418     case 22254: case 22255: return  244800/ 11.L;  // MAC HQ                  22254.54545454545454545454
1419     case 16036: case 16037: return 1411200/ 88.L;  // SB 128 "16 kHz"         16036.36363636363636363636
1420     case 11959: case 11960: return 1411200/118.L;  // SB 128 "12 kHz"         11959.32203389830508474576
1421     case 11127: case 11128: return  122400/ 11.L;  // MAC LQ                  11127.27272727272727272727
1422     case  8018: case  8019: return 1411200/176.L;  // SB 128  "8 kHz"          8018.18181818181818181818
1423     case  8012: case  8013: return  312500/ 39.L;  // NeXT/Telco               8012.82051282051282051282
1424     case  5512: case  5513: return   44100/  8.L;  // 1/8 CD                   5512.50000000000000000000
1425     }
1426 
1427     return freq;
1428 }
1429 
1430 /*}}}*/
1431 /*{{{ resampling stuff                */
1432 
Calc_Coeffs(IN sample_t * Coeff,OUT int iNew,OUT int iOld,OUT long i,OUT long k,OUT double bandwidth)1433 static inline void  Calc_Coeffs (
1434         IN  sample_t*  Coeff,
1435         OUT int        iNew,
1436         OUT int        iOld,
1437         OUT long       i,
1438         OUT long       k,
1439         OUT double     bandwidth )   // 0.0 ... 1.0: used bandwidth from 0...fs_out/2
1440 {
1441     long double  w   = 1.L / (WINDOW_SIZE + 0.5);
1442     long double  sum = 0.;
1443     long double  tmp1;
1444     long double  tmp2;
1445     long         j;
1446 
1447     for ( j = 0; j <= 2*WINDOW_SIZE; j++ ) {
1448         tmp1 = ((k+j)*iNew - i*iOld) / (long double)iNew * bandwidth;
1449         tmp2 = ((k+j)*iNew - i*iOld) / (long double)iNew;
1450         sum += Coeff [j] = sinc (tmp1) * WINDOW (tmp2 * w);
1451     }
1452 
1453     for ( j = 0; j <= 2*WINDOW_SIZE; j++ )
1454         Coeff [j] /= sum;
1455 }
1456 
1457 
resample_open(OUT long double sampfreq_in,OUT long double sampfreq_out,OUT double lowpass_freq,OUT int quality)1458 resample_t*  resample_open (
1459         OUT long double  sampfreq_in,    // [Hz]
1460         OUT long double  sampfreq_out,   // [Hz]
1461         OUT double       lowpass_freq,   // [Hz] or <0 for auto mode
1462         OUT int          quality )       // Proposal: 0: default, 1: sample select, 2: linear interpolation
1463                                          //           4: 4-point interpolation, 32: 32-point interpolation
1464 {
1465     resample_t* const  ret = calloc ( sizeof(resample_t), 1 );
1466     long               i;
1467     sample_t*          p;
1468     double             err;
1469 
1470     FN("resample_open");
1471 
1472     if ( sampfreq_in == sampfreq_out )
1473         return NULL;
1474 
1475     err = factorize ( sampfreq_out, sampfreq_in, &(ret->scale_out), &(ret->scale_in) );
1476     fprintf ( stderr, "(%.6g => %.6g Hz, Ratio %d:%d",
1477               (double)sampfreq_in, (double)sampfreq_out, ret->scale_in, ret->scale_out );
1478     if ( fabs (err) > 5.e-10 )                                            // 16 msec/year
1479         fprintf ( stderr, ", Error %+.3f ppm", 1.e6*err );
1480     fprintf ( stderr, ")\n" );
1481     fflush  ( stderr );
1482 
1483     if ( ret->scale_out == ret->scale_in )
1484         return NULL;
1485 
1486     ret->Class_ID        = RESAMPLE_ID;
1487     ret->samplefreq_in   = sampfreq_in;
1488     ret->samplefreq_out  = sampfreq_out;
1489     ret->taps            = TAPS;
1490     ret->lowpass_freq    = lowpass_freq < 0.  ?  0.50 * MIN (sampfreq_in, sampfreq_out)  :  lowpass_freq;
1491     ret->in_old [0]      = calloc ( sizeof(sample_t*)     , 2*TAPS );
1492     ret->in_old [1]      = calloc ( sizeof(sample_t*)     , 2*TAPS );
1493     ret->src_step        = calloc ( sizeof(unsigned char*), (unsigned)ret->scale_out );
1494     ret->fir             = calloc ( sizeof(sample_t*)     , (unsigned)ret->scale_out );
1495     ret->firfree         = calloc ( sizeof(sample_t)      , TAPS * ret->scale_out + 64/sizeof(sample_t) );
1496     p                    = (sample_t*) ( ((ptrdiff_t)(ret->firfree) | 63) + 1 ); /* 64 byte alignment */
1497 
1498     for ( i = 0; i < ret->scale_out; i++, p += TAPS ) {
1499         ret->src_step [i] = round_nearest ( (i+1.L) * ret->scale_in / ret->scale_out )
1500                           - round_nearest ( (i+0.L) * ret->scale_in / ret->scale_out );
1501 
1502         Calc_Coeffs ( ret->fir [i] = p,
1503                       ret->scale_out, ret->scale_in,
1504                       i, round_nearest ( (double)i * ret->scale_in / ret->scale_out - WINDOW_SIZE ),
1505                       2.*ret->lowpass_freq / sampfreq_out );
1506     }
1507 
1508     return ret;
1509 }
1510 
1511 
1512 /* Should be done with more sanity checks */
1513 
resample_close(INOUT resample_t * const r)1514 int  resample_close ( INOUT resample_t* const r )
1515 {
1516     FN("resample_close");
1517 
1518     if ( r == NULL )
1519         return -1;
1520     free ( r->in_old [0] );
1521     free ( r->in_old [1] );
1522     free ( r->fir [0]    );
1523     free ( r->fir        );
1524     free ( r->src_step   );
1525     return 0;
1526 }
1527 
1528 
1529 // in the future some of the copies should be avoided by usage of mfbuf
1530 // also by the fill_buffer_resample() routine
1531 
1532 /*
1533  *  Current working scheme:
1534  *
1535  *  | old | new                data in in_buf, consist of old *and new data
1536  *  |0 1 2|3 4 5 6 ...
1537  *
1538  *   <--1-->
1539  *     <--2-->
1540  *       <--3-->
1541  *
1542  *        | new                data in in, only new data
1543  *        |0 1 2 3 4 5 6
1544  *         <--4-->
1545  *           <--5-->
1546  *             <--6-->
1547  */
1548 
resample_buffer(INOUT resample_t * const r,IN sample_t * const out,INOUT size_t * const out_req_len,OUT sample_t * const in,INOUT size_t * const in_avail_len,OUT size_t channel)1549 int  resample_buffer (                          // return code, 0 for success
1550         INOUT resample_t *const   r,            // internal structure
1551         IN    sample_t   *const   out,          // where to write the output data
1552         INOUT size_t     *const   out_req_len,  // requested output data len/really written output data len
1553         OUT   sample_t   *const   in,           // where are the input data?
1554         INOUT size_t     *const   in_avail_len, // available input data len/consumed input data len
1555         OUT   size_t              channel )     // number of the channel (needed for internal buffering)
1556 {
1557     sample_t*  p           = out;
1558     sample_t*  in_old      = r->in_old [channel];
1559     size_t     len         = *in_avail_len;
1560     size_t     desired_len = *out_req_len;
1561     size_t     i;
1562     size_t     k;
1563     int        l;
1564 
1565     FN("resample_buffer");
1566 
1567     // copy TAPS samples to the second half of in_old
1568     memcpy ( in_old + TAPS, in, MIN (len*sizeof(*in), TAPS*sizeof(*in)) );
1569 
1570     // calculating k and l
1571     l = r->fir_stepper [channel];
1572     k = r->inp_stepper [channel];
1573 
1574     // doing converting for all samples which can be done in 'in_old'
1575     for ( i = 0; i < desired_len  &&  k < TAPS; i++ ) {
1576         *p++ = scalar32 ( in_old + k, r->fir [l] );
1577         k   += r->src_step [l];
1578         if ( ++l >= r->scale_out ) l = 0;
1579     }
1580 
1581     // doing converting for all samples which can be done in 'in'
1582     k -= TAPS;
1583     for ( ; i < desired_len  &&  k < len - TAPS; i++ ) {
1584         *p++ = scalar32 ( in + k, r->fir [l] );
1585         k   += r->src_step [l];
1586         if ( ++l >= r->scale_out ) l = 0;
1587     }
1588 
1589     // writing back processed in and out
1590     *in_avail_len  = k - r->inp_stepper [channel];
1591     *out_req_len   = p - out;
1592 
1593     r->fir_stepper [channel] = l;
1594     r->inp_stepper [channel] = k - r->inp_stepper [channel] - len;
1595 
1596     // make of copy of the overlap
1597     if ( len >= TAPS )
1598         memcpy  ( in_old, in + len - TAPS, TAPS*sizeof(sample_t) );
1599     else
1600         memmove ( in_old, in_old + len   , TAPS*sizeof(sample_t) );
1601 
1602     return 0;
1603 }
1604 
1605 /*}}}*/
1606 
1607 #endif /* KLEMM_44 */
1608 
1609 /* end of pcm.c */
1610