celp_math: Replace duplicate ff_dot_productf() by ff_scalarproduct_c()
[ffmpeg.git] / libavcodec / wmavoice.c
1 /*
2  * Windows Media Audio Voice decoder.
3  * Copyright (c) 2009 Ronald S. Bultje
4  *
5  * This file is part of Libav.
6  *
7  * Libav is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public
9  * License as published by the Free Software Foundation; either
10  * version 2.1 of the License, or (at your option) any later version.
11  *
12  * Libav 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  * Lesser General Public License for more details.
16  *
17  * You should have received a copy of the GNU Lesser General Public
18  * License along with Libav; if not, write to the Free Software
19  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
20  */
21
22 /**
23  * @file
24  * @brief Windows Media Audio Voice compatible decoder
25  * @author Ronald S. Bultje <rsbultje@gmail.com>
26  */
27
28 #define UNCHECKED_BITSTREAM_READER 1
29
30 #include <math.h>
31
32 #include "dsputil.h"
33 #include "avcodec.h"
34 #include "get_bits.h"
35 #include "put_bits.h"
36 #include "wmavoice_data.h"
37 #include "celp_filters.h"
38 #include "acelp_vectors.h"
39 #include "acelp_filters.h"
40 #include "lsp.h"
41 #include "libavutil/lzo.h"
42 #include "dct.h"
43 #include "rdft.h"
44 #include "sinewin.h"
45
46 #define MAX_BLOCKS           8   ///< maximum number of blocks per frame
47 #define MAX_LSPS             16  ///< maximum filter order
48 #define MAX_LSPS_ALIGN16     16  ///< same as #MAX_LSPS; needs to be multiple
49                                  ///< of 16 for ASM input buffer alignment
50 #define MAX_FRAMES           3   ///< maximum number of frames per superframe
51 #define MAX_FRAMESIZE        160 ///< maximum number of samples per frame
52 #define MAX_SIGNAL_HISTORY   416 ///< maximum excitation signal history
53 #define MAX_SFRAMESIZE       (MAX_FRAMESIZE * MAX_FRAMES)
54                                  ///< maximum number of samples per superframe
55 #define SFRAME_CACHE_MAXSIZE 256 ///< maximum cache size for frame data that
56                                  ///< was split over two packets
57 #define VLC_NBITS            6   ///< number of bits to read per VLC iteration
58
59 /**
60  * Frame type VLC coding.
61  */
62 static VLC frame_type_vlc;
63
64 /**
65  * Adaptive codebook types.
66  */
67 enum {
68     ACB_TYPE_NONE       = 0, ///< no adaptive codebook (only hardcoded fixed)
69     ACB_TYPE_ASYMMETRIC = 1, ///< adaptive codebook with per-frame pitch, which
70                              ///< we interpolate to get a per-sample pitch.
71                              ///< Signal is generated using an asymmetric sinc
72                              ///< window function
73                              ///< @note see #wmavoice_ipol1_coeffs
74     ACB_TYPE_HAMMING    = 2  ///< Per-block pitch with signal generation using
75                              ///< a Hamming sinc window function
76                              ///< @note see #wmavoice_ipol2_coeffs
77 };
78
79 /**
80  * Fixed codebook types.
81  */
82 enum {
83     FCB_TYPE_SILENCE    = 0, ///< comfort noise during silence
84                              ///< generated from a hardcoded (fixed) codebook
85                              ///< with per-frame (low) gain values
86     FCB_TYPE_HARDCODED  = 1, ///< hardcoded (fixed) codebook with per-block
87                              ///< gain values
88     FCB_TYPE_AW_PULSES  = 2, ///< Pitch-adaptive window (AW) pulse signals,
89                              ///< used in particular for low-bitrate streams
90     FCB_TYPE_EXC_PULSES = 3, ///< Innovation (fixed) codebook pulse sets in
91                              ///< combinations of either single pulses or
92                              ///< pulse pairs
93 };
94
95 /**
96  * Description of frame types.
97  */
98 static const struct frame_type_desc {
99     uint8_t n_blocks;     ///< amount of blocks per frame (each block
100                           ///< (contains 160/#n_blocks samples)
101     uint8_t log_n_blocks; ///< log2(#n_blocks)
102     uint8_t acb_type;     ///< Adaptive codebook type (ACB_TYPE_*)
103     uint8_t fcb_type;     ///< Fixed codebook type (FCB_TYPE_*)
104     uint8_t dbl_pulses;   ///< how many pulse vectors have pulse pairs
105                           ///< (rather than just one single pulse)
106                           ///< only if #fcb_type == #FCB_TYPE_EXC_PULSES
107     uint16_t frame_size;  ///< the amount of bits that make up the block
108                           ///< data (per frame)
109 } frame_descs[17] = {
110     { 1, 0, ACB_TYPE_NONE,       FCB_TYPE_SILENCE,    0,   0 },
111     { 2, 1, ACB_TYPE_NONE,       FCB_TYPE_HARDCODED,  0,  28 },
112     { 2, 1, ACB_TYPE_ASYMMETRIC, FCB_TYPE_AW_PULSES,  0,  46 },
113     { 2, 1, ACB_TYPE_ASYMMETRIC, FCB_TYPE_EXC_PULSES, 2,  80 },
114     { 2, 1, ACB_TYPE_ASYMMETRIC, FCB_TYPE_EXC_PULSES, 5, 104 },
115     { 4, 2, ACB_TYPE_ASYMMETRIC, FCB_TYPE_EXC_PULSES, 0, 108 },
116     { 4, 2, ACB_TYPE_ASYMMETRIC, FCB_TYPE_EXC_PULSES, 2, 132 },
117     { 4, 2, ACB_TYPE_ASYMMETRIC, FCB_TYPE_EXC_PULSES, 5, 168 },
118     { 2, 1, ACB_TYPE_HAMMING,    FCB_TYPE_EXC_PULSES, 0,  64 },
119     { 2, 1, ACB_TYPE_HAMMING,    FCB_TYPE_EXC_PULSES, 2,  80 },
120     { 2, 1, ACB_TYPE_HAMMING,    FCB_TYPE_EXC_PULSES, 5, 104 },
121     { 4, 2, ACB_TYPE_HAMMING,    FCB_TYPE_EXC_PULSES, 0, 108 },
122     { 4, 2, ACB_TYPE_HAMMING,    FCB_TYPE_EXC_PULSES, 2, 132 },
123     { 4, 2, ACB_TYPE_HAMMING,    FCB_TYPE_EXC_PULSES, 5, 168 },
124     { 8, 3, ACB_TYPE_HAMMING,    FCB_TYPE_EXC_PULSES, 0, 176 },
125     { 8, 3, ACB_TYPE_HAMMING,    FCB_TYPE_EXC_PULSES, 2, 208 },
126     { 8, 3, ACB_TYPE_HAMMING,    FCB_TYPE_EXC_PULSES, 5, 256 }
127 };
128
129 /**
130  * WMA Voice decoding context.
131  */
132 typedef struct {
133     /**
134      * @name Global values specified in the stream header / extradata or used all over.
135      * @{
136      */
137     AVFrame frame;
138     GetBitContext gb;             ///< packet bitreader. During decoder init,
139                                   ///< it contains the extradata from the
140                                   ///< demuxer. During decoding, it contains
141                                   ///< packet data.
142     int8_t vbm_tree[25];          ///< converts VLC codes to frame type
143
144     int spillover_bitsize;        ///< number of bits used to specify
145                                   ///< #spillover_nbits in the packet header
146                                   ///< = ceil(log2(ctx->block_align << 3))
147     int history_nsamples;         ///< number of samples in history for signal
148                                   ///< prediction (through ACB)
149
150     /* postfilter specific values */
151     int do_apf;                   ///< whether to apply the averaged
152                                   ///< projection filter (APF)
153     int denoise_strength;         ///< strength of denoising in Wiener filter
154                                   ///< [0-11]
155     int denoise_tilt_corr;        ///< Whether to apply tilt correction to the
156                                   ///< Wiener filter coefficients (postfilter)
157     int dc_level;                 ///< Predicted amount of DC noise, based
158                                   ///< on which a DC removal filter is used
159
160     int lsps;                     ///< number of LSPs per frame [10 or 16]
161     int lsp_q_mode;               ///< defines quantizer defaults [0, 1]
162     int lsp_def_mode;             ///< defines different sets of LSP defaults
163                                   ///< [0, 1]
164     int frame_lsp_bitsize;        ///< size (in bits) of LSPs, when encoded
165                                   ///< per-frame (independent coding)
166     int sframe_lsp_bitsize;       ///< size (in bits) of LSPs, when encoded
167                                   ///< per superframe (residual coding)
168
169     int min_pitch_val;            ///< base value for pitch parsing code
170     int max_pitch_val;            ///< max value + 1 for pitch parsing
171     int pitch_nbits;              ///< number of bits used to specify the
172                                   ///< pitch value in the frame header
173     int block_pitch_nbits;        ///< number of bits used to specify the
174                                   ///< first block's pitch value
175     int block_pitch_range;        ///< range of the block pitch
176     int block_delta_pitch_nbits;  ///< number of bits used to specify the
177                                   ///< delta pitch between this and the last
178                                   ///< block's pitch value, used in all but
179                                   ///< first block
180     int block_delta_pitch_hrange; ///< 1/2 range of the delta (full range is
181                                   ///< from -this to +this-1)
182     uint16_t block_conv_table[4]; ///< boundaries for block pitch unit/scale
183                                   ///< conversion
184
185     /**
186      * @}
187      *
188      * @name Packet values specified in the packet header or related to a packet.
189      *
190      * A packet is considered to be a single unit of data provided to this
191      * decoder by the demuxer.
192      * @{
193      */
194     int spillover_nbits;          ///< number of bits of the previous packet's
195                                   ///< last superframe preceding this
196                                   ///< packet's first full superframe (useful
197                                   ///< for re-synchronization also)
198     int has_residual_lsps;        ///< if set, superframes contain one set of
199                                   ///< LSPs that cover all frames, encoded as
200                                   ///< independent and residual LSPs; if not
201                                   ///< set, each frame contains its own, fully
202                                   ///< independent, LSPs
203     int skip_bits_next;           ///< number of bits to skip at the next call
204                                   ///< to #wmavoice_decode_packet() (since
205                                   ///< they're part of the previous superframe)
206
207     uint8_t sframe_cache[SFRAME_CACHE_MAXSIZE + FF_INPUT_BUFFER_PADDING_SIZE];
208                                   ///< cache for superframe data split over
209                                   ///< multiple packets
210     int sframe_cache_size;        ///< set to >0 if we have data from an
211                                   ///< (incomplete) superframe from a previous
212                                   ///< packet that spilled over in the current
213                                   ///< packet; specifies the amount of bits in
214                                   ///< #sframe_cache
215     PutBitContext pb;             ///< bitstream writer for #sframe_cache
216
217     /**
218      * @}
219      *
220      * @name Frame and superframe values
221      * Superframe and frame data - these can change from frame to frame,
222      * although some of them do in that case serve as a cache / history for
223      * the next frame or superframe.
224      * @{
225      */
226     double prev_lsps[MAX_LSPS];   ///< LSPs of the last frame of the previous
227                                   ///< superframe
228     int last_pitch_val;           ///< pitch value of the previous frame
229     int last_acb_type;            ///< frame type [0-2] of the previous frame
230     int pitch_diff_sh16;          ///< ((cur_pitch_val - #last_pitch_val)
231                                   ///< << 16) / #MAX_FRAMESIZE
232     float silence_gain;           ///< set for use in blocks if #ACB_TYPE_NONE
233
234     int aw_idx_is_ext;            ///< whether the AW index was encoded in
235                                   ///< 8 bits (instead of 6)
236     int aw_pulse_range;           ///< the range over which #aw_pulse_set1()
237                                   ///< can apply the pulse, relative to the
238                                   ///< value in aw_first_pulse_off. The exact
239                                   ///< position of the first AW-pulse is within
240                                   ///< [pulse_off, pulse_off + this], and
241                                   ///< depends on bitstream values; [16 or 24]
242     int aw_n_pulses[2];           ///< number of AW-pulses in each block; note
243                                   ///< that this number can be negative (in
244                                   ///< which case it basically means "zero")
245     int aw_first_pulse_off[2];    ///< index of first sample to which to
246                                   ///< apply AW-pulses, or -0xff if unset
247     int aw_next_pulse_off_cache;  ///< the position (relative to start of the
248                                   ///< second block) at which pulses should
249                                   ///< start to be positioned, serves as a
250                                   ///< cache for pitch-adaptive window pulses
251                                   ///< between blocks
252
253     int frame_cntr;               ///< current frame index [0 - 0xFFFE]; is
254                                   ///< only used for comfort noise in #pRNG()
255     float gain_pred_err[6];       ///< cache for gain prediction
256     float excitation_history[MAX_SIGNAL_HISTORY];
257                                   ///< cache of the signal of previous
258                                   ///< superframes, used as a history for
259                                   ///< signal generation
260     float synth_history[MAX_LSPS]; ///< see #excitation_history
261     /**
262      * @}
263      *
264      * @name Postfilter values
265      *
266      * Variables used for postfilter implementation, mostly history for
267      * smoothing and so on, and context variables for FFT/iFFT.
268      * @{
269      */
270     RDFTContext rdft, irdft;      ///< contexts for FFT-calculation in the
271                                   ///< postfilter (for denoise filter)
272     DCTContext dct, dst;          ///< contexts for phase shift (in Hilbert
273                                   ///< transform, part of postfilter)
274     float sin[511], cos[511];     ///< 8-bit cosine/sine windows over [-pi,pi]
275                                   ///< range
276     float postfilter_agc;         ///< gain control memory, used in
277                                   ///< #adaptive_gain_control()
278     float dcf_mem[2];             ///< DC filter history
279     float zero_exc_pf[MAX_SIGNAL_HISTORY + MAX_SFRAMESIZE];
280                                   ///< zero filter output (i.e. excitation)
281                                   ///< by postfilter
282     float denoise_filter_cache[MAX_FRAMESIZE];
283     int   denoise_filter_cache_size; ///< samples in #denoise_filter_cache
284     DECLARE_ALIGNED(32, float, tilted_lpcs_pf)[0x80];
285                                   ///< aligned buffer for LPC tilting
286     DECLARE_ALIGNED(32, float, denoise_coeffs_pf)[0x80];
287                                   ///< aligned buffer for denoise coefficients
288     DECLARE_ALIGNED(32, float, synth_filter_out_buf)[0x80 + MAX_LSPS_ALIGN16];
289                                   ///< aligned buffer for postfilter speech
290                                   ///< synthesis
291     /**
292      * @}
293      */
294 } WMAVoiceContext;
295
296 /**
297  * Set up the variable bit mode (VBM) tree from container extradata.
298  * @param gb bit I/O context.
299  *           The bit context (s->gb) should be loaded with byte 23-46 of the
300  *           container extradata (i.e. the ones containing the VBM tree).
301  * @param vbm_tree pointer to array to which the decoded VBM tree will be
302  *                 written.
303  * @return 0 on success, <0 on error.
304  */
305 static av_cold int decode_vbmtree(GetBitContext *gb, int8_t vbm_tree[25])
306 {
307     static const uint8_t bits[] = {
308          2,  2,  2,  4,  4,  4,
309          6,  6,  6,  8,  8,  8,
310         10, 10, 10, 12, 12, 12,
311         14, 14, 14, 14
312     };
313     static const uint16_t codes[] = {
314           0x0000, 0x0001, 0x0002,        //              00/01/10
315           0x000c, 0x000d, 0x000e,        //           11+00/01/10
316           0x003c, 0x003d, 0x003e,        //         1111+00/01/10
317           0x00fc, 0x00fd, 0x00fe,        //       111111+00/01/10
318           0x03fc, 0x03fd, 0x03fe,        //     11111111+00/01/10
319           0x0ffc, 0x0ffd, 0x0ffe,        //   1111111111+00/01/10
320           0x3ffc, 0x3ffd, 0x3ffe, 0x3fff // 111111111111+xx
321     };
322     int cntr[8] = { 0 }, n, res;
323
324     memset(vbm_tree, 0xff, sizeof(vbm_tree[0]) * 25);
325     for (n = 0; n < 17; n++) {
326         res = get_bits(gb, 3);
327         if (cntr[res] > 3) // should be >= 3 + (res == 7))
328             return -1;
329         vbm_tree[res * 3 + cntr[res]++] = n;
330     }
331     INIT_VLC_STATIC(&frame_type_vlc, VLC_NBITS, sizeof(bits),
332                     bits, 1, 1, codes, 2, 2, 132);
333     return 0;
334 }
335
336 /**
337  * Set up decoder with parameters from demuxer (extradata etc.).
338  */
339 static av_cold int wmavoice_decode_init(AVCodecContext *ctx)
340 {
341     int n, flags, pitch_range, lsp16_flag;
342     WMAVoiceContext *s = ctx->priv_data;
343
344     /**
345      * Extradata layout:
346      * - byte  0-18: WMAPro-in-WMAVoice extradata (see wmaprodec.c),
347      * - byte 19-22: flags field (annoyingly in LE; see below for known
348      *               values),
349      * - byte 23-46: variable bitmode tree (really just 17 * 3 bits,
350      *               rest is 0).
351      */
352     if (ctx->extradata_size != 46) {
353         av_log(ctx, AV_LOG_ERROR,
354                "Invalid extradata size %d (should be 46)\n",
355                ctx->extradata_size);
356         return -1;
357     }
358     flags                = AV_RL32(ctx->extradata + 18);
359     s->spillover_bitsize = 3 + av_ceil_log2(ctx->block_align);
360     s->do_apf            =    flags & 0x1;
361     if (s->do_apf) {
362         ff_rdft_init(&s->rdft,  7, DFT_R2C);
363         ff_rdft_init(&s->irdft, 7, IDFT_C2R);
364         ff_dct_init(&s->dct,  6, DCT_I);
365         ff_dct_init(&s->dst,  6, DST_I);
366
367         ff_sine_window_init(s->cos, 256);
368         memcpy(&s->sin[255], s->cos, 256 * sizeof(s->cos[0]));
369         for (n = 0; n < 255; n++) {
370             s->sin[n]       = -s->sin[510 - n];
371             s->cos[510 - n] =  s->cos[n];
372         }
373     }
374     s->denoise_strength  =   (flags >> 2) & 0xF;
375     if (s->denoise_strength >= 12) {
376         av_log(ctx, AV_LOG_ERROR,
377                "Invalid denoise filter strength %d (max=11)\n",
378                s->denoise_strength);
379         return -1;
380     }
381     s->denoise_tilt_corr = !!(flags & 0x40);
382     s->dc_level          =   (flags >> 7) & 0xF;
383     s->lsp_q_mode        = !!(flags & 0x2000);
384     s->lsp_def_mode      = !!(flags & 0x4000);
385     lsp16_flag           =    flags & 0x1000;
386     if (lsp16_flag) {
387         s->lsps               = 16;
388         s->frame_lsp_bitsize  = 34;
389         s->sframe_lsp_bitsize = 60;
390     } else {
391         s->lsps               = 10;
392         s->frame_lsp_bitsize  = 24;
393         s->sframe_lsp_bitsize = 48;
394     }
395     for (n = 0; n < s->lsps; n++)
396         s->prev_lsps[n] = M_PI * (n + 1.0) / (s->lsps + 1.0);
397
398     init_get_bits(&s->gb, ctx->extradata + 22, (ctx->extradata_size - 22) << 3);
399     if (decode_vbmtree(&s->gb, s->vbm_tree) < 0) {
400         av_log(ctx, AV_LOG_ERROR, "Invalid VBM tree; broken extradata?\n");
401         return -1;
402     }
403
404     s->min_pitch_val    = ((ctx->sample_rate << 8)      /  400 + 50) >> 8;
405     s->max_pitch_val    = ((ctx->sample_rate << 8) * 37 / 2000 + 50) >> 8;
406     pitch_range         = s->max_pitch_val - s->min_pitch_val;
407     if (pitch_range <= 0) {
408         av_log(ctx, AV_LOG_ERROR, "Invalid pitch range; broken extradata?\n");
409         return -1;
410     }
411     s->pitch_nbits      = av_ceil_log2(pitch_range);
412     s->last_pitch_val   = 40;
413     s->last_acb_type    = ACB_TYPE_NONE;
414     s->history_nsamples = s->max_pitch_val + 8;
415
416     if (s->min_pitch_val < 1 || s->history_nsamples > MAX_SIGNAL_HISTORY) {
417         int min_sr = ((((1 << 8) - 50) * 400) + 0xFF) >> 8,
418             max_sr = ((((MAX_SIGNAL_HISTORY - 8) << 8) + 205) * 2000 / 37) >> 8;
419
420         av_log(ctx, AV_LOG_ERROR,
421                "Unsupported samplerate %d (min=%d, max=%d)\n",
422                ctx->sample_rate, min_sr, max_sr); // 322-22097 Hz
423
424         return -1;
425     }
426
427     s->block_conv_table[0]      = s->min_pitch_val;
428     s->block_conv_table[1]      = (pitch_range * 25) >> 6;
429     s->block_conv_table[2]      = (pitch_range * 44) >> 6;
430     s->block_conv_table[3]      = s->max_pitch_val - 1;
431     s->block_delta_pitch_hrange = (pitch_range >> 3) & ~0xF;
432     if (s->block_delta_pitch_hrange <= 0) {
433         av_log(ctx, AV_LOG_ERROR, "Invalid delta pitch hrange; broken extradata?\n");
434         return -1;
435     }
436     s->block_delta_pitch_nbits  = 1 + av_ceil_log2(s->block_delta_pitch_hrange);
437     s->block_pitch_range        = s->block_conv_table[2] +
438                                   s->block_conv_table[3] + 1 +
439                                   2 * (s->block_conv_table[1] - 2 * s->min_pitch_val);
440     s->block_pitch_nbits        = av_ceil_log2(s->block_pitch_range);
441
442     ctx->sample_fmt             = AV_SAMPLE_FMT_FLT;
443
444     avcodec_get_frame_defaults(&s->frame);
445     ctx->coded_frame = &s->frame;
446
447     return 0;
448 }
449
450 /**
451  * @name Postfilter functions
452  * Postfilter functions (gain control, wiener denoise filter, DC filter,
453  * kalman smoothening, plus surrounding code to wrap it)
454  * @{
455  */
456 /**
457  * Adaptive gain control (as used in postfilter).
458  *
459  * Identical to #ff_adaptive_gain_control() in acelp_vectors.c, except
460  * that the energy here is calculated using sum(abs(...)), whereas the
461  * other codecs (e.g. AMR-NB, SIPRO) use sqrt(dotproduct(...)).
462  *
463  * @param out output buffer for filtered samples
464  * @param in input buffer containing the samples as they are after the
465  *           postfilter steps so far
466  * @param speech_synth input buffer containing speech synth before postfilter
467  * @param size input buffer size
468  * @param alpha exponential filter factor
469  * @param gain_mem pointer to filter memory (single float)
470  */
471 static void adaptive_gain_control(float *out, const float *in,
472                                   const float *speech_synth,
473                                   int size, float alpha, float *gain_mem)
474 {
475     int i;
476     float speech_energy = 0.0, postfilter_energy = 0.0, gain_scale_factor;
477     float mem = *gain_mem;
478
479     for (i = 0; i < size; i++) {
480         speech_energy     += fabsf(speech_synth[i]);
481         postfilter_energy += fabsf(in[i]);
482     }
483     gain_scale_factor = (1.0 - alpha) * speech_energy / postfilter_energy;
484
485     for (i = 0; i < size; i++) {
486         mem = alpha * mem + gain_scale_factor;
487         out[i] = in[i] * mem;
488     }
489
490     *gain_mem = mem;
491 }
492
493 /**
494  * Kalman smoothing function.
495  *
496  * This function looks back pitch +/- 3 samples back into history to find
497  * the best fitting curve (that one giving the optimal gain of the two
498  * signals, i.e. the highest dot product between the two), and then
499  * uses that signal history to smoothen the output of the speech synthesis
500  * filter.
501  *
502  * @param s WMA Voice decoding context
503  * @param pitch pitch of the speech signal
504  * @param in input speech signal
505  * @param out output pointer for smoothened signal
506  * @param size input/output buffer size
507  *
508  * @returns -1 if no smoothening took place, e.g. because no optimal
509  *          fit could be found, or 0 on success.
510  */
511 static int kalman_smoothen(WMAVoiceContext *s, int pitch,
512                            const float *in, float *out, int size)
513 {
514     int n;
515     float optimal_gain = 0, dot;
516     const float *ptr = &in[-FFMAX(s->min_pitch_val, pitch - 3)],
517                 *end = &in[-FFMIN(s->max_pitch_val, pitch + 3)],
518                 *best_hist_ptr;
519
520     /* find best fitting point in history */
521     do {
522         dot = ff_scalarproduct_float_c(in, ptr, size);
523         if (dot > optimal_gain) {
524             optimal_gain  = dot;
525             best_hist_ptr = ptr;
526         }
527     } while (--ptr >= end);
528
529     if (optimal_gain <= 0)
530         return -1;
531     dot = ff_scalarproduct_float_c(best_hist_ptr, best_hist_ptr, size);
532     if (dot <= 0) // would be 1.0
533         return -1;
534
535     if (optimal_gain <= dot) {
536         dot = dot / (dot + 0.6 * optimal_gain); // 0.625-1.000
537     } else
538         dot = 0.625;
539
540     /* actual smoothing */
541     for (n = 0; n < size; n++)
542         out[n] = best_hist_ptr[n] + dot * (in[n] - best_hist_ptr[n]);
543
544     return 0;
545 }
546
547 /**
548  * Get the tilt factor of a formant filter from its transfer function
549  * @see #tilt_factor() in amrnbdec.c, which does essentially the same,
550  *      but somehow (??) it does a speech synthesis filter in the
551  *      middle, which is missing here
552  *
553  * @param lpcs LPC coefficients
554  * @param n_lpcs Size of LPC buffer
555  * @returns the tilt factor
556  */
557 static float tilt_factor(const float *lpcs, int n_lpcs)
558 {
559     float rh0, rh1;
560
561     rh0 = 1.0     + ff_scalarproduct_float_c(lpcs,  lpcs,    n_lpcs);
562     rh1 = lpcs[0] + ff_scalarproduct_float_c(lpcs, &lpcs[1], n_lpcs - 1);
563
564     return rh1 / rh0;
565 }
566
567 /**
568  * Derive denoise filter coefficients (in real domain) from the LPCs.
569  */
570 static void calc_input_response(WMAVoiceContext *s, float *lpcs,
571                                 int fcb_type, float *coeffs, int remainder)
572 {
573     float last_coeff, min = 15.0, max = -15.0;
574     float irange, angle_mul, gain_mul, range, sq;
575     int n, idx;
576
577     /* Create frequency power spectrum of speech input (i.e. RDFT of LPCs) */
578     s->rdft.rdft_calc(&s->rdft, lpcs);
579 #define log_range(var, assign) do { \
580         float tmp = log10f(assign);  var = tmp; \
581         max       = FFMAX(max, tmp); min = FFMIN(min, tmp); \
582     } while (0)
583     log_range(last_coeff,  lpcs[1]         * lpcs[1]);
584     for (n = 1; n < 64; n++)
585         log_range(lpcs[n], lpcs[n * 2]     * lpcs[n * 2] +
586                            lpcs[n * 2 + 1] * lpcs[n * 2 + 1]);
587     log_range(lpcs[0],     lpcs[0]         * lpcs[0]);
588 #undef log_range
589     range    = max - min;
590     lpcs[64] = last_coeff;
591
592     /* Now, use this spectrum to pick out these frequencies with higher
593      * (relative) power/energy (which we then take to be "not noise"),
594      * and set up a table (still in lpc[]) of (relative) gains per frequency.
595      * These frequencies will be maintained, while others ("noise") will be
596      * decreased in the filter output. */
597     irange    = 64.0 / range; // so irange*(max-value) is in the range [0, 63]
598     gain_mul  = range * (fcb_type == FCB_TYPE_HARDCODED ? (5.0 / 13.0) :
599                                                           (5.0 / 14.7));
600     angle_mul = gain_mul * (8.0 * M_LN10 / M_PI);
601     for (n = 0; n <= 64; n++) {
602         float pwr;
603
604         idx = FFMAX(0, lrint((max - lpcs[n]) * irange) - 1);
605         pwr = wmavoice_denoise_power_table[s->denoise_strength][idx];
606         lpcs[n] = angle_mul * pwr;
607
608         /* 70.57 =~ 1/log10(1.0331663) */
609         idx = (pwr * gain_mul - 0.0295) * 70.570526123;
610         if (idx > 127) { // fallback if index falls outside table range
611             coeffs[n] = wmavoice_energy_table[127] *
612                         powf(1.0331663, idx - 127);
613         } else
614             coeffs[n] = wmavoice_energy_table[FFMAX(0, idx)];
615     }
616
617     /* calculate the Hilbert transform of the gains, which we do (since this
618      * is a sinus input) by doing a phase shift (in theory, H(sin())=cos()).
619      * Hilbert_Transform(RDFT(x)) = Laplace_Transform(x), which calculates the
620      * "moment" of the LPCs in this filter. */
621     s->dct.dct_calc(&s->dct, lpcs);
622     s->dst.dct_calc(&s->dst, lpcs);
623
624     /* Split out the coefficient indexes into phase/magnitude pairs */
625     idx = 255 + av_clip(lpcs[64],               -255, 255);
626     coeffs[0]  = coeffs[0]  * s->cos[idx];
627     idx = 255 + av_clip(lpcs[64] - 2 * lpcs[63], -255, 255);
628     last_coeff = coeffs[64] * s->cos[idx];
629     for (n = 63;; n--) {
630         idx = 255 + av_clip(-lpcs[64] - 2 * lpcs[n - 1], -255, 255);
631         coeffs[n * 2 + 1] = coeffs[n] * s->sin[idx];
632         coeffs[n * 2]     = coeffs[n] * s->cos[idx];
633
634         if (!--n) break;
635
636         idx = 255 + av_clip( lpcs[64] - 2 * lpcs[n - 1], -255, 255);
637         coeffs[n * 2 + 1] = coeffs[n] * s->sin[idx];
638         coeffs[n * 2]     = coeffs[n] * s->cos[idx];
639     }
640     coeffs[1] = last_coeff;
641
642     /* move into real domain */
643     s->irdft.rdft_calc(&s->irdft, coeffs);
644
645     /* tilt correction and normalize scale */
646     memset(&coeffs[remainder], 0, sizeof(coeffs[0]) * (128 - remainder));
647     if (s->denoise_tilt_corr) {
648         float tilt_mem = 0;
649
650         coeffs[remainder - 1] = 0;
651         ff_tilt_compensation(&tilt_mem,
652                              -1.8 * tilt_factor(coeffs, remainder - 1),
653                              coeffs, remainder);
654     }
655     sq = (1.0 / 64.0) * sqrtf(1 / ff_scalarproduct_float_c(coeffs, coeffs, remainder));
656     for (n = 0; n < remainder; n++)
657         coeffs[n] *= sq;
658 }
659
660 /**
661  * This function applies a Wiener filter on the (noisy) speech signal as
662  * a means to denoise it.
663  *
664  * - take RDFT of LPCs to get the power spectrum of the noise + speech;
665  * - using this power spectrum, calculate (for each frequency) the Wiener
666  *    filter gain, which depends on the frequency power and desired level
667  *    of noise subtraction (when set too high, this leads to artifacts)
668  *    We can do this symmetrically over the X-axis (so 0-4kHz is the inverse
669  *    of 4-8kHz);
670  * - by doing a phase shift, calculate the Hilbert transform of this array
671  *    of per-frequency filter-gains to get the filtering coefficients;
672  * - smoothen/normalize/de-tilt these filter coefficients as desired;
673  * - take RDFT of noisy sound, apply the coefficients and take its IRDFT
674  *    to get the denoised speech signal;
675  * - the leftover (i.e. output of the IRDFT on denoised speech data beyond
676  *    the frame boundary) are saved and applied to subsequent frames by an
677  *    overlap-add method (otherwise you get clicking-artifacts).
678  *
679  * @param s WMA Voice decoding context
680  * @param fcb_type Frame (codebook) type
681  * @param synth_pf input: the noisy speech signal, output: denoised speech
682  *                 data; should be 16-byte aligned (for ASM purposes)
683  * @param size size of the speech data
684  * @param lpcs LPCs used to synthesize this frame's speech data
685  */
686 static void wiener_denoise(WMAVoiceContext *s, int fcb_type,
687                            float *synth_pf, int size,
688                            const float *lpcs)
689 {
690     int remainder, lim, n;
691
692     if (fcb_type != FCB_TYPE_SILENCE) {
693         float *tilted_lpcs = s->tilted_lpcs_pf,
694               *coeffs = s->denoise_coeffs_pf, tilt_mem = 0;
695
696         tilted_lpcs[0]           = 1.0;
697         memcpy(&tilted_lpcs[1], lpcs, sizeof(lpcs[0]) * s->lsps);
698         memset(&tilted_lpcs[s->lsps + 1], 0,
699                sizeof(tilted_lpcs[0]) * (128 - s->lsps - 1));
700         ff_tilt_compensation(&tilt_mem, 0.7 * tilt_factor(lpcs, s->lsps),
701                              tilted_lpcs, s->lsps + 2);
702
703         /* The IRDFT output (127 samples for 7-bit filter) beyond the frame
704          * size is applied to the next frame. All input beyond this is zero,
705          * and thus all output beyond this will go towards zero, hence we can
706          * limit to min(size-1, 127-size) as a performance consideration. */
707         remainder = FFMIN(127 - size, size - 1);
708         calc_input_response(s, tilted_lpcs, fcb_type, coeffs, remainder);
709
710         /* apply coefficients (in frequency spectrum domain), i.e. complex
711          * number multiplication */
712         memset(&synth_pf[size], 0, sizeof(synth_pf[0]) * (128 - size));
713         s->rdft.rdft_calc(&s->rdft, synth_pf);
714         s->rdft.rdft_calc(&s->rdft, coeffs);
715         synth_pf[0] *= coeffs[0];
716         synth_pf[1] *= coeffs[1];
717         for (n = 1; n < 64; n++) {
718             float v1 = synth_pf[n * 2], v2 = synth_pf[n * 2 + 1];
719             synth_pf[n * 2]     = v1 * coeffs[n * 2] - v2 * coeffs[n * 2 + 1];
720             synth_pf[n * 2 + 1] = v2 * coeffs[n * 2] + v1 * coeffs[n * 2 + 1];
721         }
722         s->irdft.rdft_calc(&s->irdft, synth_pf);
723     }
724
725     /* merge filter output with the history of previous runs */
726     if (s->denoise_filter_cache_size) {
727         lim = FFMIN(s->denoise_filter_cache_size, size);
728         for (n = 0; n < lim; n++)
729             synth_pf[n] += s->denoise_filter_cache[n];
730         s->denoise_filter_cache_size -= lim;
731         memmove(s->denoise_filter_cache, &s->denoise_filter_cache[size],
732                 sizeof(s->denoise_filter_cache[0]) * s->denoise_filter_cache_size);
733     }
734
735     /* move remainder of filter output into a cache for future runs */
736     if (fcb_type != FCB_TYPE_SILENCE) {
737         lim = FFMIN(remainder, s->denoise_filter_cache_size);
738         for (n = 0; n < lim; n++)
739             s->denoise_filter_cache[n] += synth_pf[size + n];
740         if (lim < remainder) {
741             memcpy(&s->denoise_filter_cache[lim], &synth_pf[size + lim],
742                    sizeof(s->denoise_filter_cache[0]) * (remainder - lim));
743             s->denoise_filter_cache_size = remainder;
744         }
745     }
746 }
747
748 /**
749  * Averaging projection filter, the postfilter used in WMAVoice.
750  *
751  * This uses the following steps:
752  * - A zero-synthesis filter (generate excitation from synth signal)
753  * - Kalman smoothing on excitation, based on pitch
754  * - Re-synthesized smoothened output
755  * - Iterative Wiener denoise filter
756  * - Adaptive gain filter
757  * - DC filter
758  *
759  * @param s WMAVoice decoding context
760  * @param synth Speech synthesis output (before postfilter)
761  * @param samples Output buffer for filtered samples
762  * @param size Buffer size of synth & samples
763  * @param lpcs Generated LPCs used for speech synthesis
764  * @param zero_exc_pf destination for zero synthesis filter (16-byte aligned)
765  * @param fcb_type Frame type (silence, hardcoded, AW-pulses or FCB-pulses)
766  * @param pitch Pitch of the input signal
767  */
768 static void postfilter(WMAVoiceContext *s, const float *synth,
769                        float *samples,    int size,
770                        const float *lpcs, float *zero_exc_pf,
771                        int fcb_type,      int pitch)
772 {
773     float synth_filter_in_buf[MAX_FRAMESIZE / 2],
774           *synth_pf = &s->synth_filter_out_buf[MAX_LSPS_ALIGN16],
775           *synth_filter_in = zero_exc_pf;
776
777     assert(size <= MAX_FRAMESIZE / 2);
778
779     /* generate excitation from input signal */
780     ff_celp_lp_zero_synthesis_filterf(zero_exc_pf, lpcs, synth, size, s->lsps);
781
782     if (fcb_type >= FCB_TYPE_AW_PULSES &&
783         !kalman_smoothen(s, pitch, zero_exc_pf, synth_filter_in_buf, size))
784         synth_filter_in = synth_filter_in_buf;
785
786     /* re-synthesize speech after smoothening, and keep history */
787     ff_celp_lp_synthesis_filterf(synth_pf, lpcs,
788                                  synth_filter_in, size, s->lsps);
789     memcpy(&synth_pf[-s->lsps], &synth_pf[size - s->lsps],
790            sizeof(synth_pf[0]) * s->lsps);
791
792     wiener_denoise(s, fcb_type, synth_pf, size, lpcs);
793
794     adaptive_gain_control(samples, synth_pf, synth, size, 0.99,
795                           &s->postfilter_agc);
796
797     if (s->dc_level > 8) {
798         /* remove ultra-low frequency DC noise / highpass filter;
799          * coefficients are identical to those used in SIPR decoding,
800          * and very closely resemble those used in AMR-NB decoding. */
801         ff_acelp_apply_order_2_transfer_function(samples, samples,
802             (const float[2]) { -1.99997,      1.0 },
803             (const float[2]) { -1.9330735188, 0.93589198496 },
804             0.93980580475, s->dcf_mem, size);
805     }
806 }
807 /**
808  * @}
809  */
810
811 /**
812  * Dequantize LSPs
813  * @param lsps output pointer to the array that will hold the LSPs
814  * @param num number of LSPs to be dequantized
815  * @param values quantized values, contains n_stages values
816  * @param sizes range (i.e. max value) of each quantized value
817  * @param n_stages number of dequantization runs
818  * @param table dequantization table to be used
819  * @param mul_q LSF multiplier
820  * @param base_q base (lowest) LSF values
821  */
822 static void dequant_lsps(double *lsps, int num,
823                          const uint16_t *values,
824                          const uint16_t *sizes,
825                          int n_stages, const uint8_t *table,
826                          const double *mul_q,
827                          const double *base_q)
828 {
829     int n, m;
830
831     memset(lsps, 0, num * sizeof(*lsps));
832     for (n = 0; n < n_stages; n++) {
833         const uint8_t *t_off = &table[values[n] * num];
834         double base = base_q[n], mul = mul_q[n];
835
836         for (m = 0; m < num; m++)
837             lsps[m] += base + mul * t_off[m];
838
839         table += sizes[n] * num;
840     }
841 }
842
843 /**
844  * @name LSP dequantization routines
845  * LSP dequantization routines, for 10/16LSPs and independent/residual coding.
846  * @note we assume enough bits are available, caller should check.
847  * lsp10i() consumes 24 bits; lsp10r() consumes an additional 24 bits;
848  * lsp16i() consumes 34 bits; lsp16r() consumes an additional 26 bits.
849  * @{
850  */
851 /**
852  * Parse 10 independently-coded LSPs.
853  */
854 static void dequant_lsp10i(GetBitContext *gb, double *lsps)
855 {
856     static const uint16_t vec_sizes[4] = { 256, 64, 32, 32 };
857     static const double mul_lsf[4] = {
858         5.2187144800e-3,    1.4626986422e-3,
859         9.6179549166e-4,    1.1325736225e-3
860     };
861     static const double base_lsf[4] = {
862         M_PI * -2.15522e-1, M_PI * -6.1646e-2,
863         M_PI * -3.3486e-2,  M_PI * -5.7408e-2
864     };
865     uint16_t v[4];
866
867     v[0] = get_bits(gb, 8);
868     v[1] = get_bits(gb, 6);
869     v[2] = get_bits(gb, 5);
870     v[3] = get_bits(gb, 5);
871
872     dequant_lsps(lsps, 10, v, vec_sizes, 4, wmavoice_dq_lsp10i,
873                  mul_lsf, base_lsf);
874 }
875
876 /**
877  * Parse 10 independently-coded LSPs, and then derive the tables to
878  * generate LSPs for the other frames from them (residual coding).
879  */
880 static void dequant_lsp10r(GetBitContext *gb,
881                            double *i_lsps, const double *old,
882                            double *a1, double *a2, int q_mode)
883 {
884     static const uint16_t vec_sizes[3] = { 128, 64, 64 };
885     static const double mul_lsf[3] = {
886         2.5807601174e-3,    1.2354460219e-3,   1.1763821673e-3
887     };
888     static const double base_lsf[3] = {
889         M_PI * -1.07448e-1, M_PI * -5.2706e-2, M_PI * -5.1634e-2
890     };
891     const float (*ipol_tab)[2][10] = q_mode ?
892         wmavoice_lsp10_intercoeff_b : wmavoice_lsp10_intercoeff_a;
893     uint16_t interpol, v[3];
894     int n;
895
896     dequant_lsp10i(gb, i_lsps);
897
898     interpol = get_bits(gb, 5);
899     v[0]     = get_bits(gb, 7);
900     v[1]     = get_bits(gb, 6);
901     v[2]     = get_bits(gb, 6);
902
903     for (n = 0; n < 10; n++) {
904         double delta = old[n] - i_lsps[n];
905         a1[n]        = ipol_tab[interpol][0][n] * delta + i_lsps[n];
906         a1[10 + n]   = ipol_tab[interpol][1][n] * delta + i_lsps[n];
907     }
908
909     dequant_lsps(a2, 20, v, vec_sizes, 3, wmavoice_dq_lsp10r,
910                  mul_lsf, base_lsf);
911 }
912
913 /**
914  * Parse 16 independently-coded LSPs.
915  */
916 static void dequant_lsp16i(GetBitContext *gb, double *lsps)
917 {
918     static const uint16_t vec_sizes[5] = { 256, 64, 128, 64, 128 };
919     static const double mul_lsf[5] = {
920         3.3439586280e-3,    6.9908173703e-4,
921         3.3216608306e-3,    1.0334960326e-3,
922         3.1899104283e-3
923     };
924     static const double base_lsf[5] = {
925         M_PI * -1.27576e-1, M_PI * -2.4292e-2,
926         M_PI * -1.28094e-1, M_PI * -3.2128e-2,
927         M_PI * -1.29816e-1
928     };
929     uint16_t v[5];
930
931     v[0] = get_bits(gb, 8);
932     v[1] = get_bits(gb, 6);
933     v[2] = get_bits(gb, 7);
934     v[3] = get_bits(gb, 6);
935     v[4] = get_bits(gb, 7);
936
937     dequant_lsps( lsps,     5,  v,     vec_sizes,    2,
938                  wmavoice_dq_lsp16i1,  mul_lsf,     base_lsf);
939     dequant_lsps(&lsps[5],  5, &v[2], &vec_sizes[2], 2,
940                  wmavoice_dq_lsp16i2, &mul_lsf[2], &base_lsf[2]);
941     dequant_lsps(&lsps[10], 6, &v[4], &vec_sizes[4], 1,
942                  wmavoice_dq_lsp16i3, &mul_lsf[4], &base_lsf[4]);
943 }
944
945 /**
946  * Parse 16 independently-coded LSPs, and then derive the tables to
947  * generate LSPs for the other frames from them (residual coding).
948  */
949 static void dequant_lsp16r(GetBitContext *gb,
950                            double *i_lsps, const double *old,
951                            double *a1, double *a2, int q_mode)
952 {
953     static const uint16_t vec_sizes[3] = { 128, 128, 128 };
954     static const double mul_lsf[3] = {
955         1.2232979501e-3,   1.4062241527e-3,   1.6114744851e-3
956     };
957     static const double base_lsf[3] = {
958         M_PI * -5.5830e-2, M_PI * -5.2908e-2, M_PI * -5.4776e-2
959     };
960     const float (*ipol_tab)[2][16] = q_mode ?
961         wmavoice_lsp16_intercoeff_b : wmavoice_lsp16_intercoeff_a;
962     uint16_t interpol, v[3];
963     int n;
964
965     dequant_lsp16i(gb, i_lsps);
966
967     interpol = get_bits(gb, 5);
968     v[0]     = get_bits(gb, 7);
969     v[1]     = get_bits(gb, 7);
970     v[2]     = get_bits(gb, 7);
971
972     for (n = 0; n < 16; n++) {
973         double delta = old[n] - i_lsps[n];
974         a1[n]        = ipol_tab[interpol][0][n] * delta + i_lsps[n];
975         a1[16 + n]   = ipol_tab[interpol][1][n] * delta + i_lsps[n];
976     }
977
978     dequant_lsps( a2,     10,  v,     vec_sizes,    1,
979                  wmavoice_dq_lsp16r1,  mul_lsf,     base_lsf);
980     dequant_lsps(&a2[10], 10, &v[1], &vec_sizes[1], 1,
981                  wmavoice_dq_lsp16r2, &mul_lsf[1], &base_lsf[1]);
982     dequant_lsps(&a2[20], 12, &v[2], &vec_sizes[2], 1,
983                  wmavoice_dq_lsp16r3, &mul_lsf[2], &base_lsf[2]);
984 }
985
986 /**
987  * @}
988  * @name Pitch-adaptive window coding functions
989  * The next few functions are for pitch-adaptive window coding.
990  * @{
991  */
992 /**
993  * Parse the offset of the first pitch-adaptive window pulses, and
994  * the distribution of pulses between the two blocks in this frame.
995  * @param s WMA Voice decoding context private data
996  * @param gb bit I/O context
997  * @param pitch pitch for each block in this frame
998  */
999 static void aw_parse_coords(WMAVoiceContext *s, GetBitContext *gb,
1000                             const int *pitch)
1001 {
1002     static const int16_t start_offset[94] = {
1003         -11,  -9,  -7,  -5,  -3,  -1,   1,   3,   5,   7,   9,  11,
1004          13,  15,  18,  17,  19,  20,  21,  22,  23,  24,  25,  26,
1005          27,  28,  29,  30,  31,  32,  33,  35,  37,  39,  41,  43,
1006          45,  47,  49,  51,  53,  55,  57,  59,  61,  63,  65,  67,
1007          69,  71,  73,  75,  77,  79,  81,  83,  85,  87,  89,  91,
1008          93,  95,  97,  99, 101, 103, 105, 107, 109, 111, 113, 115,
1009         117, 119, 121, 123, 125, 127, 129, 131, 133, 135, 137, 139,
1010         141, 143, 145, 147, 149, 151, 153, 155, 157, 159
1011     };
1012     int bits, offset;
1013
1014     /* position of pulse */
1015     s->aw_idx_is_ext = 0;
1016     if ((bits = get_bits(gb, 6)) >= 54) {
1017         s->aw_idx_is_ext = 1;
1018         bits += (bits - 54) * 3 + get_bits(gb, 2);
1019     }
1020
1021     /* for a repeated pulse at pulse_off with a pitch_lag of pitch[], count
1022      * the distribution of the pulses in each block contained in this frame. */
1023     s->aw_pulse_range        = FFMIN(pitch[0], pitch[1]) > 32 ? 24 : 16;
1024     for (offset = start_offset[bits]; offset < 0; offset += pitch[0]) ;
1025     s->aw_n_pulses[0]        = (pitch[0] - 1 + MAX_FRAMESIZE / 2 - offset) / pitch[0];
1026     s->aw_first_pulse_off[0] = offset - s->aw_pulse_range / 2;
1027     offset                  += s->aw_n_pulses[0] * pitch[0];
1028     s->aw_n_pulses[1]        = (pitch[1] - 1 + MAX_FRAMESIZE - offset) / pitch[1];
1029     s->aw_first_pulse_off[1] = offset - (MAX_FRAMESIZE + s->aw_pulse_range) / 2;
1030
1031     /* if continuing from a position before the block, reset position to
1032      * start of block (when corrected for the range over which it can be
1033      * spread in aw_pulse_set1()). */
1034     if (start_offset[bits] < MAX_FRAMESIZE / 2) {
1035         while (s->aw_first_pulse_off[1] - pitch[1] + s->aw_pulse_range > 0)
1036             s->aw_first_pulse_off[1] -= pitch[1];
1037         if (start_offset[bits] < 0)
1038             while (s->aw_first_pulse_off[0] - pitch[0] + s->aw_pulse_range > 0)
1039                 s->aw_first_pulse_off[0] -= pitch[0];
1040     }
1041 }
1042
1043 /**
1044  * Apply second set of pitch-adaptive window pulses.
1045  * @param s WMA Voice decoding context private data
1046  * @param gb bit I/O context
1047  * @param block_idx block index in frame [0, 1]
1048  * @param fcb structure containing fixed codebook vector info
1049  */
1050 static void aw_pulse_set2(WMAVoiceContext *s, GetBitContext *gb,
1051                           int block_idx, AMRFixed *fcb)
1052 {
1053     uint16_t use_mask_mem[9]; // only 5 are used, rest is padding
1054     uint16_t *use_mask = use_mask_mem + 2;
1055     /* in this function, idx is the index in the 80-bit (+ padding) use_mask
1056      * bit-array. Since use_mask consists of 16-bit values, the lower 4 bits
1057      * of idx are the position of the bit within a particular item in the
1058      * array (0 being the most significant bit, and 15 being the least
1059      * significant bit), and the remainder (>> 4) is the index in the
1060      * use_mask[]-array. This is faster and uses less memory than using a
1061      * 80-byte/80-int array. */
1062     int pulse_off = s->aw_first_pulse_off[block_idx],
1063         pulse_start, n, idx, range, aidx, start_off = 0;
1064
1065     /* set offset of first pulse to within this block */
1066     if (s->aw_n_pulses[block_idx] > 0)
1067         while (pulse_off + s->aw_pulse_range < 1)
1068             pulse_off += fcb->pitch_lag;
1069
1070     /* find range per pulse */
1071     if (s->aw_n_pulses[0] > 0) {
1072         if (block_idx == 0) {
1073             range = 32;
1074         } else /* block_idx = 1 */ {
1075             range = 8;
1076             if (s->aw_n_pulses[block_idx] > 0)
1077                 pulse_off = s->aw_next_pulse_off_cache;
1078         }
1079     } else
1080         range = 16;
1081     pulse_start = s->aw_n_pulses[block_idx] > 0 ? pulse_off - range / 2 : 0;
1082
1083     /* aw_pulse_set1() already applies pulses around pulse_off (to be exactly,
1084      * in the range of [pulse_off, pulse_off + s->aw_pulse_range], and thus
1085      * we exclude that range from being pulsed again in this function. */
1086     memset(&use_mask[-2], 0, 2 * sizeof(use_mask[0]));
1087     memset( use_mask,   -1, 5 * sizeof(use_mask[0]));
1088     memset(&use_mask[5], 0, 2 * sizeof(use_mask[0]));
1089     if (s->aw_n_pulses[block_idx] > 0)
1090         for (idx = pulse_off; idx < MAX_FRAMESIZE / 2; idx += fcb->pitch_lag) {
1091             int excl_range         = s->aw_pulse_range; // always 16 or 24
1092             uint16_t *use_mask_ptr = &use_mask[idx >> 4];
1093             int first_sh           = 16 - (idx & 15);
1094             *use_mask_ptr++       &= 0xFFFFu << first_sh;
1095             excl_range            -= first_sh;
1096             if (excl_range >= 16) {
1097                 *use_mask_ptr++    = 0;
1098                 *use_mask_ptr     &= 0xFFFF >> (excl_range - 16);
1099             } else
1100                 *use_mask_ptr     &= 0xFFFF >> excl_range;
1101         }
1102
1103     /* find the 'aidx'th offset that is not excluded */
1104     aidx = get_bits(gb, s->aw_n_pulses[0] > 0 ? 5 - 2 * block_idx : 4);
1105     for (n = 0; n <= aidx; pulse_start++) {
1106         for (idx = pulse_start; idx < 0; idx += fcb->pitch_lag) ;
1107         if (idx >= MAX_FRAMESIZE / 2) { // find from zero
1108             if (use_mask[0])      idx = 0x0F;
1109             else if (use_mask[1]) idx = 0x1F;
1110             else if (use_mask[2]) idx = 0x2F;
1111             else if (use_mask[3]) idx = 0x3F;
1112             else if (use_mask[4]) idx = 0x4F;
1113             else                  return;
1114             idx -= av_log2_16bit(use_mask[idx >> 4]);
1115         }
1116         if (use_mask[idx >> 4] & (0x8000 >> (idx & 15))) {
1117             use_mask[idx >> 4] &= ~(0x8000 >> (idx & 15));
1118             n++;
1119             start_off = idx;
1120         }
1121     }
1122
1123     fcb->x[fcb->n] = start_off;
1124     fcb->y[fcb->n] = get_bits1(gb) ? -1.0 : 1.0;
1125     fcb->n++;
1126
1127     /* set offset for next block, relative to start of that block */
1128     n = (MAX_FRAMESIZE / 2 - start_off) % fcb->pitch_lag;
1129     s->aw_next_pulse_off_cache = n ? fcb->pitch_lag - n : 0;
1130 }
1131
1132 /**
1133  * Apply first set of pitch-adaptive window pulses.
1134  * @param s WMA Voice decoding context private data
1135  * @param gb bit I/O context
1136  * @param block_idx block index in frame [0, 1]
1137  * @param fcb storage location for fixed codebook pulse info
1138  */
1139 static void aw_pulse_set1(WMAVoiceContext *s, GetBitContext *gb,
1140                           int block_idx, AMRFixed *fcb)
1141 {
1142     int val = get_bits(gb, 12 - 2 * (s->aw_idx_is_ext && !block_idx));
1143     float v;
1144
1145     if (s->aw_n_pulses[block_idx] > 0) {
1146         int n, v_mask, i_mask, sh, n_pulses;
1147
1148         if (s->aw_pulse_range == 24) { // 3 pulses, 1:sign + 3:index each
1149             n_pulses = 3;
1150             v_mask   = 8;
1151             i_mask   = 7;
1152             sh       = 4;
1153         } else { // 4 pulses, 1:sign + 2:index each
1154             n_pulses = 4;
1155             v_mask   = 4;
1156             i_mask   = 3;
1157             sh       = 3;
1158         }
1159
1160         for (n = n_pulses - 1; n >= 0; n--, val >>= sh) {
1161             fcb->y[fcb->n] = (val & v_mask) ? -1.0 : 1.0;
1162             fcb->x[fcb->n] = (val & i_mask) * n_pulses + n +
1163                                  s->aw_first_pulse_off[block_idx];
1164             while (fcb->x[fcb->n] < 0)
1165                 fcb->x[fcb->n] += fcb->pitch_lag;
1166             if (fcb->x[fcb->n] < MAX_FRAMESIZE / 2)
1167                 fcb->n++;
1168         }
1169     } else {
1170         int num2 = (val & 0x1FF) >> 1, delta, idx;
1171
1172         if (num2 < 1 * 79)      { delta = 1; idx = num2 + 1; }
1173         else if (num2 < 2 * 78) { delta = 3; idx = num2 + 1 - 1 * 77; }
1174         else if (num2 < 3 * 77) { delta = 5; idx = num2 + 1 - 2 * 76; }
1175         else                    { delta = 7; idx = num2 + 1 - 3 * 75; }
1176         v = (val & 0x200) ? -1.0 : 1.0;
1177
1178         fcb->no_repeat_mask |= 3 << fcb->n;
1179         fcb->x[fcb->n]       = idx - delta;
1180         fcb->y[fcb->n]       = v;
1181         fcb->x[fcb->n + 1]   = idx;
1182         fcb->y[fcb->n + 1]   = (val & 1) ? -v : v;
1183         fcb->n              += 2;
1184     }
1185 }
1186
1187 /**
1188  * @}
1189  *
1190  * Generate a random number from frame_cntr and block_idx, which will lief
1191  * in the range [0, 1000 - block_size] (so it can be used as an index in a
1192  * table of size 1000 of which you want to read block_size entries).
1193  *
1194  * @param frame_cntr current frame number
1195  * @param block_num current block index
1196  * @param block_size amount of entries we want to read from a table
1197  *                   that has 1000 entries
1198  * @return a (non-)random number in the [0, 1000 - block_size] range.
1199  */
1200 static int pRNG(int frame_cntr, int block_num, int block_size)
1201 {
1202     /* array to simplify the calculation of z:
1203      * y = (x % 9) * 5 + 6;
1204      * z = (49995 * x) / y;
1205      * Since y only has 9 values, we can remove the division by using a
1206      * LUT and using FASTDIV-style divisions. For each of the 9 values
1207      * of y, we can rewrite z as:
1208      * z = x * (49995 / y) + x * ((49995 % y) / y)
1209      * In this table, each col represents one possible value of y, the
1210      * first number is 49995 / y, and the second is the FASTDIV variant
1211      * of 49995 % y / y. */
1212     static const unsigned int div_tbl[9][2] = {
1213         { 8332,  3 * 715827883U }, // y =  6
1214         { 4545,  0 * 390451573U }, // y = 11
1215         { 3124, 11 * 268435456U }, // y = 16
1216         { 2380, 15 * 204522253U }, // y = 21
1217         { 1922, 23 * 165191050U }, // y = 26
1218         { 1612, 23 * 138547333U }, // y = 31
1219         { 1388, 27 * 119304648U }, // y = 36
1220         { 1219, 16 * 104755300U }, // y = 41
1221         { 1086, 39 *  93368855U }  // y = 46
1222     };
1223     unsigned int z, y, x = MUL16(block_num, 1877) + frame_cntr;
1224     if (x >= 0xFFFF) x -= 0xFFFF;   // max value of x is 8*1877+0xFFFE=0x13AA6,
1225                                     // so this is effectively a modulo (%)
1226     y = x - 9 * MULH(477218589, x); // x % 9
1227     z = (uint16_t) (x * div_tbl[y][0] + UMULH(x, div_tbl[y][1]));
1228                                     // z = x * 49995 / (y * 5 + 6)
1229     return z % (1000 - block_size);
1230 }
1231
1232 /**
1233  * Parse hardcoded signal for a single block.
1234  * @note see #synth_block().
1235  */
1236 static void synth_block_hardcoded(WMAVoiceContext *s, GetBitContext *gb,
1237                                  int block_idx, int size,
1238                                  const struct frame_type_desc *frame_desc,
1239                                  float *excitation)
1240 {
1241     float gain;
1242     int n, r_idx;
1243
1244     assert(size <= MAX_FRAMESIZE);
1245
1246     /* Set the offset from which we start reading wmavoice_std_codebook */
1247     if (frame_desc->fcb_type == FCB_TYPE_SILENCE) {
1248         r_idx = pRNG(s->frame_cntr, block_idx, size);
1249         gain  = s->silence_gain;
1250     } else /* FCB_TYPE_HARDCODED */ {
1251         r_idx = get_bits(gb, 8);
1252         gain  = wmavoice_gain_universal[get_bits(gb, 6)];
1253     }
1254
1255     /* Clear gain prediction parameters */
1256     memset(s->gain_pred_err, 0, sizeof(s->gain_pred_err));
1257
1258     /* Apply gain to hardcoded codebook and use that as excitation signal */
1259     for (n = 0; n < size; n++)
1260         excitation[n] = wmavoice_std_codebook[r_idx + n] * gain;
1261 }
1262
1263 /**
1264  * Parse FCB/ACB signal for a single block.
1265  * @note see #synth_block().
1266  */
1267 static void synth_block_fcb_acb(WMAVoiceContext *s, GetBitContext *gb,
1268                                 int block_idx, int size,
1269                                 int block_pitch_sh2,
1270                                 const struct frame_type_desc *frame_desc,
1271                                 float *excitation)
1272 {
1273     static const float gain_coeff[6] = {
1274         0.8169, -0.06545, 0.1726, 0.0185, -0.0359, 0.0458
1275     };
1276     float pulses[MAX_FRAMESIZE / 2], pred_err, acb_gain, fcb_gain;
1277     int n, idx, gain_weight;
1278     AMRFixed fcb;
1279
1280     assert(size <= MAX_FRAMESIZE / 2);
1281     memset(pulses, 0, sizeof(*pulses) * size);
1282
1283     fcb.pitch_lag      = block_pitch_sh2 >> 2;
1284     fcb.pitch_fac      = 1.0;
1285     fcb.no_repeat_mask = 0;
1286     fcb.n              = 0;
1287
1288     /* For the other frame types, this is where we apply the innovation
1289      * (fixed) codebook pulses of the speech signal. */
1290     if (frame_desc->fcb_type == FCB_TYPE_AW_PULSES) {
1291         aw_pulse_set1(s, gb, block_idx, &fcb);
1292         aw_pulse_set2(s, gb, block_idx, &fcb);
1293     } else /* FCB_TYPE_EXC_PULSES */ {
1294         int offset_nbits = 5 - frame_desc->log_n_blocks;
1295
1296         fcb.no_repeat_mask = -1;
1297         /* similar to ff_decode_10_pulses_35bits(), but with single pulses
1298          * (instead of double) for a subset of pulses */
1299         for (n = 0; n < 5; n++) {
1300             float sign;
1301             int pos1, pos2;
1302
1303             sign           = get_bits1(gb) ? 1.0 : -1.0;
1304             pos1           = get_bits(gb, offset_nbits);
1305             fcb.x[fcb.n]   = n + 5 * pos1;
1306             fcb.y[fcb.n++] = sign;
1307             if (n < frame_desc->dbl_pulses) {
1308                 pos2           = get_bits(gb, offset_nbits);
1309                 fcb.x[fcb.n]   = n + 5 * pos2;
1310                 fcb.y[fcb.n++] = (pos1 < pos2) ? -sign : sign;
1311             }
1312         }
1313     }
1314     ff_set_fixed_vector(pulses, &fcb, 1.0, size);
1315
1316     /* Calculate gain for adaptive & fixed codebook signal.
1317      * see ff_amr_set_fixed_gain(). */
1318     idx = get_bits(gb, 7);
1319     fcb_gain = expf(ff_scalarproduct_float_c(s->gain_pred_err, gain_coeff, 6) -
1320                     5.2409161640 + wmavoice_gain_codebook_fcb[idx]);
1321     acb_gain = wmavoice_gain_codebook_acb[idx];
1322     pred_err = av_clipf(wmavoice_gain_codebook_fcb[idx],
1323                         -2.9957322736 /* log(0.05) */,
1324                          1.6094379124 /* log(5.0)  */);
1325
1326     gain_weight = 8 >> frame_desc->log_n_blocks;
1327     memmove(&s->gain_pred_err[gain_weight], s->gain_pred_err,
1328             sizeof(*s->gain_pred_err) * (6 - gain_weight));
1329     for (n = 0; n < gain_weight; n++)
1330         s->gain_pred_err[n] = pred_err;
1331
1332     /* Calculation of adaptive codebook */
1333     if (frame_desc->acb_type == ACB_TYPE_ASYMMETRIC) {
1334         int len;
1335         for (n = 0; n < size; n += len) {
1336             int next_idx_sh16;
1337             int abs_idx    = block_idx * size + n;
1338             int pitch_sh16 = (s->last_pitch_val << 16) +
1339                              s->pitch_diff_sh16 * abs_idx;
1340             int pitch      = (pitch_sh16 + 0x6FFF) >> 16;
1341             int idx_sh16   = ((pitch << 16) - pitch_sh16) * 8 + 0x58000;
1342             idx            = idx_sh16 >> 16;
1343             if (s->pitch_diff_sh16) {
1344                 if (s->pitch_diff_sh16 > 0) {
1345                     next_idx_sh16 = (idx_sh16) &~ 0xFFFF;
1346                 } else
1347                     next_idx_sh16 = (idx_sh16 + 0x10000) &~ 0xFFFF;
1348                 len = av_clip((idx_sh16 - next_idx_sh16) / s->pitch_diff_sh16 / 8,
1349                               1, size - n);
1350             } else
1351                 len = size;
1352
1353             ff_acelp_interpolatef(&excitation[n], &excitation[n - pitch],
1354                                   wmavoice_ipol1_coeffs, 17,
1355                                   idx, 9, len);
1356         }
1357     } else /* ACB_TYPE_HAMMING */ {
1358         int block_pitch = block_pitch_sh2 >> 2;
1359         idx             = block_pitch_sh2 & 3;
1360         if (idx) {
1361             ff_acelp_interpolatef(excitation, &excitation[-block_pitch],
1362                                   wmavoice_ipol2_coeffs, 4,
1363                                   idx, 8, size);
1364         } else
1365             av_memcpy_backptr((uint8_t *) excitation, sizeof(float) * block_pitch,
1366                               sizeof(float) * size);
1367     }
1368
1369     /* Interpolate ACB/FCB and use as excitation signal */
1370     ff_weighted_vector_sumf(excitation, excitation, pulses,
1371                             acb_gain, fcb_gain, size);
1372 }
1373
1374 /**
1375  * Parse data in a single block.
1376  * @note we assume enough bits are available, caller should check.
1377  *
1378  * @param s WMA Voice decoding context private data
1379  * @param gb bit I/O context
1380  * @param block_idx index of the to-be-read block
1381  * @param size amount of samples to be read in this block
1382  * @param block_pitch_sh2 pitch for this block << 2
1383  * @param lsps LSPs for (the end of) this frame
1384  * @param prev_lsps LSPs for the last frame
1385  * @param frame_desc frame type descriptor
1386  * @param excitation target memory for the ACB+FCB interpolated signal
1387  * @param synth target memory for the speech synthesis filter output
1388  * @return 0 on success, <0 on error.
1389  */
1390 static void synth_block(WMAVoiceContext *s, GetBitContext *gb,
1391                         int block_idx, int size,
1392                         int block_pitch_sh2,
1393                         const double *lsps, const double *prev_lsps,
1394                         const struct frame_type_desc *frame_desc,
1395                         float *excitation, float *synth)
1396 {
1397     double i_lsps[MAX_LSPS];
1398     float lpcs[MAX_LSPS];
1399     float fac;
1400     int n;
1401
1402     if (frame_desc->acb_type == ACB_TYPE_NONE)
1403         synth_block_hardcoded(s, gb, block_idx, size, frame_desc, excitation);
1404     else
1405         synth_block_fcb_acb(s, gb, block_idx, size, block_pitch_sh2,
1406                             frame_desc, excitation);
1407
1408     /* convert interpolated LSPs to LPCs */
1409     fac = (block_idx + 0.5) / frame_desc->n_blocks;
1410     for (n = 0; n < s->lsps; n++) // LSF -> LSP
1411         i_lsps[n] = cos(prev_lsps[n] + fac * (lsps[n] - prev_lsps[n]));
1412     ff_acelp_lspd2lpc(i_lsps, lpcs, s->lsps >> 1);
1413
1414     /* Speech synthesis */
1415     ff_celp_lp_synthesis_filterf(synth, lpcs, excitation, size, s->lsps);
1416 }
1417
1418 /**
1419  * Synthesize output samples for a single frame.
1420  * @note we assume enough bits are available, caller should check.
1421  *
1422  * @param ctx WMA Voice decoder context
1423  * @param gb bit I/O context (s->gb or one for cross-packet superframes)
1424  * @param frame_idx Frame number within superframe [0-2]
1425  * @param samples pointer to output sample buffer, has space for at least 160
1426  *                samples
1427  * @param lsps LSP array
1428  * @param prev_lsps array of previous frame's LSPs
1429  * @param excitation target buffer for excitation signal
1430  * @param synth target buffer for synthesized speech data
1431  * @return 0 on success, <0 on error.
1432  */
1433 static int synth_frame(AVCodecContext *ctx, GetBitContext *gb, int frame_idx,
1434                        float *samples,
1435                        const double *lsps, const double *prev_lsps,
1436                        float *excitation, float *synth)
1437 {
1438     WMAVoiceContext *s = ctx->priv_data;
1439     int n, n_blocks_x2, log_n_blocks_x2, cur_pitch_val;
1440     int pitch[MAX_BLOCKS], last_block_pitch;
1441
1442     /* Parse frame type ("frame header"), see frame_descs */
1443     int bd_idx = s->vbm_tree[get_vlc2(gb, frame_type_vlc.table, 6, 3)], block_nsamples;
1444
1445     if (bd_idx < 0) {
1446         av_log(ctx, AV_LOG_ERROR,
1447                "Invalid frame type VLC code, skipping\n");
1448         return -1;
1449     }
1450
1451     block_nsamples = MAX_FRAMESIZE / frame_descs[bd_idx].n_blocks;
1452
1453     /* Pitch calculation for ACB_TYPE_ASYMMETRIC ("pitch-per-frame") */
1454     if (frame_descs[bd_idx].acb_type == ACB_TYPE_ASYMMETRIC) {
1455         /* Pitch is provided per frame, which is interpreted as the pitch of
1456          * the last sample of the last block of this frame. We can interpolate
1457          * the pitch of other blocks (and even pitch-per-sample) by gradually
1458          * incrementing/decrementing prev_frame_pitch to cur_pitch_val. */
1459         n_blocks_x2      = frame_descs[bd_idx].n_blocks << 1;
1460         log_n_blocks_x2  = frame_descs[bd_idx].log_n_blocks + 1;
1461         cur_pitch_val    = s->min_pitch_val + get_bits(gb, s->pitch_nbits);
1462         cur_pitch_val    = FFMIN(cur_pitch_val, s->max_pitch_val - 1);
1463         if (s->last_acb_type == ACB_TYPE_NONE ||
1464             20 * abs(cur_pitch_val - s->last_pitch_val) >
1465                 (cur_pitch_val + s->last_pitch_val))
1466             s->last_pitch_val = cur_pitch_val;
1467
1468         /* pitch per block */
1469         for (n = 0; n < frame_descs[bd_idx].n_blocks; n++) {
1470             int fac = n * 2 + 1;
1471
1472             pitch[n] = (MUL16(fac,                 cur_pitch_val) +
1473                         MUL16((n_blocks_x2 - fac), s->last_pitch_val) +
1474                         frame_descs[bd_idx].n_blocks) >> log_n_blocks_x2;
1475         }
1476
1477         /* "pitch-diff-per-sample" for calculation of pitch per sample */
1478         s->pitch_diff_sh16 =
1479             ((cur_pitch_val - s->last_pitch_val) << 16) / MAX_FRAMESIZE;
1480     }
1481
1482     /* Global gain (if silence) and pitch-adaptive window coordinates */
1483     switch (frame_descs[bd_idx].fcb_type) {
1484     case FCB_TYPE_SILENCE:
1485         s->silence_gain = wmavoice_gain_silence[get_bits(gb, 8)];
1486         break;
1487     case FCB_TYPE_AW_PULSES:
1488         aw_parse_coords(s, gb, pitch);
1489         break;
1490     }
1491
1492     for (n = 0; n < frame_descs[bd_idx].n_blocks; n++) {
1493         int bl_pitch_sh2;
1494
1495         /* Pitch calculation for ACB_TYPE_HAMMING ("pitch-per-block") */
1496         switch (frame_descs[bd_idx].acb_type) {
1497         case ACB_TYPE_HAMMING: {
1498             /* Pitch is given per block. Per-block pitches are encoded as an
1499              * absolute value for the first block, and then delta values
1500              * relative to this value) for all subsequent blocks. The scale of
1501              * this pitch value is semi-logaritmic compared to its use in the
1502              * decoder, so we convert it to normal scale also. */
1503             int block_pitch,
1504                 t1 = (s->block_conv_table[1] - s->block_conv_table[0]) << 2,
1505                 t2 = (s->block_conv_table[2] - s->block_conv_table[1]) << 1,
1506                 t3 =  s->block_conv_table[3] - s->block_conv_table[2] + 1;
1507
1508             if (n == 0) {
1509                 block_pitch = get_bits(gb, s->block_pitch_nbits);
1510             } else
1511                 block_pitch = last_block_pitch - s->block_delta_pitch_hrange +
1512                                  get_bits(gb, s->block_delta_pitch_nbits);
1513             /* Convert last_ so that any next delta is within _range */
1514             last_block_pitch = av_clip(block_pitch,
1515                                        s->block_delta_pitch_hrange,
1516                                        s->block_pitch_range -
1517                                            s->block_delta_pitch_hrange);
1518
1519             /* Convert semi-log-style scale back to normal scale */
1520             if (block_pitch < t1) {
1521                 bl_pitch_sh2 = (s->block_conv_table[0] << 2) + block_pitch;
1522             } else {
1523                 block_pitch -= t1;
1524                 if (block_pitch < t2) {
1525                     bl_pitch_sh2 =
1526                         (s->block_conv_table[1] << 2) + (block_pitch << 1);
1527                 } else {
1528                     block_pitch -= t2;
1529                     if (block_pitch < t3) {
1530                         bl_pitch_sh2 =
1531                             (s->block_conv_table[2] + block_pitch) << 2;
1532                     } else
1533                         bl_pitch_sh2 = s->block_conv_table[3] << 2;
1534                 }
1535             }
1536             pitch[n] = bl_pitch_sh2 >> 2;
1537             break;
1538         }
1539
1540         case ACB_TYPE_ASYMMETRIC: {
1541             bl_pitch_sh2 = pitch[n] << 2;
1542             break;
1543         }
1544
1545         default: // ACB_TYPE_NONE has no pitch
1546             bl_pitch_sh2 = 0;
1547             break;
1548         }
1549
1550         synth_block(s, gb, n, block_nsamples, bl_pitch_sh2,
1551                     lsps, prev_lsps, &frame_descs[bd_idx],
1552                     &excitation[n * block_nsamples],
1553                     &synth[n * block_nsamples]);
1554     }
1555
1556     /* Averaging projection filter, if applicable. Else, just copy samples
1557      * from synthesis buffer */
1558     if (s->do_apf) {
1559         double i_lsps[MAX_LSPS];
1560         float lpcs[MAX_LSPS];
1561
1562         for (n = 0; n < s->lsps; n++) // LSF -> LSP
1563             i_lsps[n] = cos(0.5 * (prev_lsps[n] + lsps[n]));
1564         ff_acelp_lspd2lpc(i_lsps, lpcs, s->lsps >> 1);
1565         postfilter(s, synth, samples, 80, lpcs,
1566                    &s->zero_exc_pf[s->history_nsamples + MAX_FRAMESIZE * frame_idx],
1567                    frame_descs[bd_idx].fcb_type, pitch[0]);
1568
1569         for (n = 0; n < s->lsps; n++) // LSF -> LSP
1570             i_lsps[n] = cos(lsps[n]);
1571         ff_acelp_lspd2lpc(i_lsps, lpcs, s->lsps >> 1);
1572         postfilter(s, &synth[80], &samples[80], 80, lpcs,
1573                    &s->zero_exc_pf[s->history_nsamples + MAX_FRAMESIZE * frame_idx + 80],
1574                    frame_descs[bd_idx].fcb_type, pitch[0]);
1575     } else
1576         memcpy(samples, synth, 160 * sizeof(synth[0]));
1577
1578     /* Cache values for next frame */
1579     s->frame_cntr++;
1580     if (s->frame_cntr >= 0xFFFF) s->frame_cntr -= 0xFFFF; // i.e. modulo (%)
1581     s->last_acb_type = frame_descs[bd_idx].acb_type;
1582     switch (frame_descs[bd_idx].acb_type) {
1583     case ACB_TYPE_NONE:
1584         s->last_pitch_val = 0;
1585         break;
1586     case ACB_TYPE_ASYMMETRIC:
1587         s->last_pitch_val = cur_pitch_val;
1588         break;
1589     case ACB_TYPE_HAMMING:
1590         s->last_pitch_val = pitch[frame_descs[bd_idx].n_blocks - 1];
1591         break;
1592     }
1593
1594     return 0;
1595 }
1596
1597 /**
1598  * Ensure minimum value for first item, maximum value for last value,
1599  * proper spacing between each value and proper ordering.
1600  *
1601  * @param lsps array of LSPs
1602  * @param num size of LSP array
1603  *
1604  * @note basically a double version of #ff_acelp_reorder_lsf(), might be
1605  *       useful to put in a generic location later on. Parts are also
1606  *       present in #ff_set_min_dist_lsf() + #ff_sort_nearly_sorted_floats(),
1607  *       which is in float.
1608  */
1609 static void stabilize_lsps(double *lsps, int num)
1610 {
1611     int n, m, l;
1612
1613     /* set minimum value for first, maximum value for last and minimum
1614      * spacing between LSF values.
1615      * Very similar to ff_set_min_dist_lsf(), but in double. */
1616     lsps[0]       = FFMAX(lsps[0],       0.0015 * M_PI);
1617     for (n = 1; n < num; n++)
1618         lsps[n]   = FFMAX(lsps[n],       lsps[n - 1] + 0.0125 * M_PI);
1619     lsps[num - 1] = FFMIN(lsps[num - 1], 0.9985 * M_PI);
1620
1621     /* reorder (looks like one-time / non-recursed bubblesort).
1622      * Very similar to ff_sort_nearly_sorted_floats(), but in double. */
1623     for (n = 1; n < num; n++) {
1624         if (lsps[n] < lsps[n - 1]) {
1625             for (m = 1; m < num; m++) {
1626                 double tmp = lsps[m];
1627                 for (l = m - 1; l >= 0; l--) {
1628                     if (lsps[l] <= tmp) break;
1629                     lsps[l + 1] = lsps[l];
1630                 }
1631                 lsps[l + 1] = tmp;
1632             }
1633             break;
1634         }
1635     }
1636 }
1637
1638 /**
1639  * Test if there's enough bits to read 1 superframe.
1640  *
1641  * @param orig_gb bit I/O context used for reading. This function
1642  *                does not modify the state of the bitreader; it
1643  *                only uses it to copy the current stream position
1644  * @param s WMA Voice decoding context private data
1645  * @return -1 if unsupported, 1 on not enough bits or 0 if OK.
1646  */
1647 static int check_bits_for_superframe(GetBitContext *orig_gb,
1648                                      WMAVoiceContext *s)
1649 {
1650     GetBitContext s_gb, *gb = &s_gb;
1651     int n, need_bits, bd_idx;
1652     const struct frame_type_desc *frame_desc;
1653
1654     /* initialize a copy */
1655     init_get_bits(gb, orig_gb->buffer, orig_gb->size_in_bits);
1656     skip_bits_long(gb, get_bits_count(orig_gb));
1657     assert(get_bits_left(gb) == get_bits_left(orig_gb));
1658
1659     /* superframe header */
1660     if (get_bits_left(gb) < 14)
1661         return 1;
1662     if (!get_bits1(gb))
1663         return -1;                        // WMAPro-in-WMAVoice superframe
1664     if (get_bits1(gb)) skip_bits(gb, 12); // number of  samples in superframe
1665     if (s->has_residual_lsps) {           // residual LSPs (for all frames)
1666         if (get_bits_left(gb) < s->sframe_lsp_bitsize)
1667             return 1;
1668         skip_bits_long(gb, s->sframe_lsp_bitsize);
1669     }
1670
1671     /* frames */
1672     for (n = 0; n < MAX_FRAMES; n++) {
1673         int aw_idx_is_ext = 0;
1674
1675         if (!s->has_residual_lsps) {     // independent LSPs (per-frame)
1676            if (get_bits_left(gb) < s->frame_lsp_bitsize) return 1;
1677            skip_bits_long(gb, s->frame_lsp_bitsize);
1678         }
1679         bd_idx = s->vbm_tree[get_vlc2(gb, frame_type_vlc.table, 6, 3)];
1680         if (bd_idx < 0)
1681             return -1;                   // invalid frame type VLC code
1682         frame_desc = &frame_descs[bd_idx];
1683         if (frame_desc->acb_type == ACB_TYPE_ASYMMETRIC) {
1684             if (get_bits_left(gb) < s->pitch_nbits)
1685                 return 1;
1686             skip_bits_long(gb, s->pitch_nbits);
1687         }
1688         if (frame_desc->fcb_type == FCB_TYPE_SILENCE) {
1689             skip_bits(gb, 8);
1690         } else if (frame_desc->fcb_type == FCB_TYPE_AW_PULSES) {
1691             int tmp = get_bits(gb, 6);
1692             if (tmp >= 0x36) {
1693                 skip_bits(gb, 2);
1694                 aw_idx_is_ext = 1;
1695             }
1696         }
1697
1698         /* blocks */
1699         if (frame_desc->acb_type == ACB_TYPE_HAMMING) {
1700             need_bits = s->block_pitch_nbits +
1701                 (frame_desc->n_blocks - 1) * s->block_delta_pitch_nbits;
1702         } else if (frame_desc->fcb_type == FCB_TYPE_AW_PULSES) {
1703             need_bits = 2 * !aw_idx_is_ext;
1704         } else
1705             need_bits = 0;
1706         need_bits += frame_desc->frame_size;
1707         if (get_bits_left(gb) < need_bits)
1708             return 1;
1709         skip_bits_long(gb, need_bits);
1710     }
1711
1712     return 0;
1713 }
1714
1715 /**
1716  * Synthesize output samples for a single superframe. If we have any data
1717  * cached in s->sframe_cache, that will be used instead of whatever is loaded
1718  * in s->gb.
1719  *
1720  * WMA Voice superframes contain 3 frames, each containing 160 audio samples,
1721  * to give a total of 480 samples per frame. See #synth_frame() for frame
1722  * parsing. In addition to 3 frames, superframes can also contain the LSPs
1723  * (if these are globally specified for all frames (residually); they can
1724  * also be specified individually per-frame. See the s->has_residual_lsps
1725  * option), and can specify the number of samples encoded in this superframe
1726  * (if less than 480), usually used to prevent blanks at track boundaries.
1727  *
1728  * @param ctx WMA Voice decoder context
1729  * @return 0 on success, <0 on error or 1 if there was not enough data to
1730  *         fully parse the superframe
1731  */
1732 static int synth_superframe(AVCodecContext *ctx, int *got_frame_ptr)
1733 {
1734     WMAVoiceContext *s = ctx->priv_data;
1735     GetBitContext *gb = &s->gb, s_gb;
1736     int n, res, n_samples = 480;
1737     double lsps[MAX_FRAMES][MAX_LSPS];
1738     const double *mean_lsf = s->lsps == 16 ?
1739         wmavoice_mean_lsf16[s->lsp_def_mode] : wmavoice_mean_lsf10[s->lsp_def_mode];
1740     float excitation[MAX_SIGNAL_HISTORY + MAX_SFRAMESIZE + 12];
1741     float synth[MAX_LSPS + MAX_SFRAMESIZE];
1742     float *samples;
1743
1744     memcpy(synth,      s->synth_history,
1745            s->lsps             * sizeof(*synth));
1746     memcpy(excitation, s->excitation_history,
1747            s->history_nsamples * sizeof(*excitation));
1748
1749     if (s->sframe_cache_size > 0) {
1750         gb = &s_gb;
1751         init_get_bits(gb, s->sframe_cache, s->sframe_cache_size);
1752         s->sframe_cache_size = 0;
1753     }
1754
1755     if ((res = check_bits_for_superframe(gb, s)) == 1) {
1756         *got_frame_ptr = 0;
1757         return 1;
1758     }
1759
1760     /* First bit is speech/music bit, it differentiates between WMAVoice
1761      * speech samples (the actual codec) and WMAVoice music samples, which
1762      * are really WMAPro-in-WMAVoice-superframes. I've never seen those in
1763      * the wild yet. */
1764     if (!get_bits1(gb)) {
1765         av_log_missing_feature(ctx, "WMAPro-in-WMAVoice support", 1);
1766         return -1;
1767     }
1768
1769     /* (optional) nr. of samples in superframe; always <= 480 and >= 0 */
1770     if (get_bits1(gb)) {
1771         if ((n_samples = get_bits(gb, 12)) > 480) {
1772             av_log(ctx, AV_LOG_ERROR,
1773                    "Superframe encodes >480 samples (%d), not allowed\n",
1774                    n_samples);
1775             return -1;
1776         }
1777     }
1778     /* Parse LSPs, if global for the superframe (can also be per-frame). */
1779     if (s->has_residual_lsps) {
1780         double prev_lsps[MAX_LSPS], a1[MAX_LSPS * 2], a2[MAX_LSPS * 2];
1781
1782         for (n = 0; n < s->lsps; n++)
1783             prev_lsps[n] = s->prev_lsps[n] - mean_lsf[n];
1784
1785         if (s->lsps == 10) {
1786             dequant_lsp10r(gb, lsps[2], prev_lsps, a1, a2, s->lsp_q_mode);
1787         } else /* s->lsps == 16 */
1788             dequant_lsp16r(gb, lsps[2], prev_lsps, a1, a2, s->lsp_q_mode);
1789
1790         for (n = 0; n < s->lsps; n++) {
1791             lsps[0][n]  = mean_lsf[n] + (a1[n]           - a2[n * 2]);
1792             lsps[1][n]  = mean_lsf[n] + (a1[s->lsps + n] - a2[n * 2 + 1]);
1793             lsps[2][n] += mean_lsf[n];
1794         }
1795         for (n = 0; n < 3; n++)
1796             stabilize_lsps(lsps[n], s->lsps);
1797     }
1798
1799     /* get output buffer */
1800     s->frame.nb_samples = 480;
1801     if ((res = ctx->get_buffer(ctx, &s->frame)) < 0) {
1802         av_log(ctx, AV_LOG_ERROR, "get_buffer() failed\n");
1803         return res;
1804     }
1805     s->frame.nb_samples = n_samples;
1806     samples = (float *)s->frame.data[0];
1807
1808     /* Parse frames, optionally preceded by per-frame (independent) LSPs. */
1809     for (n = 0; n < 3; n++) {
1810         if (!s->has_residual_lsps) {
1811             int m;
1812
1813             if (s->lsps == 10) {
1814                 dequant_lsp10i(gb, lsps[n]);
1815             } else /* s->lsps == 16 */
1816                 dequant_lsp16i(gb, lsps[n]);
1817
1818             for (m = 0; m < s->lsps; m++)
1819                 lsps[n][m] += mean_lsf[m];
1820             stabilize_lsps(lsps[n], s->lsps);
1821         }
1822
1823         if ((res = synth_frame(ctx, gb, n,
1824                                &samples[n * MAX_FRAMESIZE],
1825                                lsps[n], n == 0 ? s->prev_lsps : lsps[n - 1],
1826                                &excitation[s->history_nsamples + n * MAX_FRAMESIZE],
1827                                &synth[s->lsps + n * MAX_FRAMESIZE]))) {
1828             *got_frame_ptr = 0;
1829             return res;
1830         }
1831     }
1832
1833     /* Statistics? FIXME - we don't check for length, a slight overrun
1834      * will be caught by internal buffer padding, and anything else
1835      * will be skipped, not read. */
1836     if (get_bits1(gb)) {
1837         res = get_bits(gb, 4);
1838         skip_bits(gb, 10 * (res + 1));
1839     }
1840
1841     *got_frame_ptr = 1;
1842
1843     /* Update history */
1844     memcpy(s->prev_lsps,           lsps[2],
1845            s->lsps             * sizeof(*s->prev_lsps));
1846     memcpy(s->synth_history,      &synth[MAX_SFRAMESIZE],
1847            s->lsps             * sizeof(*synth));
1848     memcpy(s->excitation_history, &excitation[MAX_SFRAMESIZE],
1849            s->history_nsamples * sizeof(*excitation));
1850     if (s->do_apf)
1851         memmove(s->zero_exc_pf,       &s->zero_exc_pf[MAX_SFRAMESIZE],
1852                 s->history_nsamples * sizeof(*s->zero_exc_pf));
1853
1854     return 0;
1855 }
1856
1857 /**
1858  * Parse the packet header at the start of each packet (input data to this
1859  * decoder).
1860  *
1861  * @param s WMA Voice decoding context private data
1862  * @return 1 if not enough bits were available, or 0 on success.
1863  */
1864 static int parse_packet_header(WMAVoiceContext *s)
1865 {
1866     GetBitContext *gb = &s->gb;
1867     unsigned int res;
1868
1869     if (get_bits_left(gb) < 11)
1870         return 1;
1871     skip_bits(gb, 4);          // packet sequence number
1872     s->has_residual_lsps = get_bits1(gb);
1873     do {
1874         res = get_bits(gb, 6); // number of superframes per packet
1875                                // (minus first one if there is spillover)
1876         if (get_bits_left(gb) < 6 * (res == 0x3F) + s->spillover_bitsize)
1877             return 1;
1878     } while (res == 0x3F);
1879     s->spillover_nbits   = get_bits(gb, s->spillover_bitsize);
1880
1881     return 0;
1882 }
1883
1884 /**
1885  * Copy (unaligned) bits from gb/data/size to pb.
1886  *
1887  * @param pb target buffer to copy bits into
1888  * @param data source buffer to copy bits from
1889  * @param size size of the source data, in bytes
1890  * @param gb bit I/O context specifying the current position in the source.
1891  *           data. This function might use this to align the bit position to
1892  *           a whole-byte boundary before calling #avpriv_copy_bits() on aligned
1893  *           source data
1894  * @param nbits the amount of bits to copy from source to target
1895  *
1896  * @note after calling this function, the current position in the input bit
1897  *       I/O context is undefined.
1898  */
1899 static void copy_bits(PutBitContext *pb,
1900                       const uint8_t *data, int size,
1901                       GetBitContext *gb, int nbits)
1902 {
1903     int rmn_bytes, rmn_bits;
1904
1905     rmn_bits = rmn_bytes = get_bits_left(gb);
1906     if (rmn_bits < nbits)
1907         return;
1908     if (nbits > pb->size_in_bits - put_bits_count(pb))
1909         return;
1910     rmn_bits &= 7; rmn_bytes >>= 3;
1911     if ((rmn_bits = FFMIN(rmn_bits, nbits)) > 0)
1912         put_bits(pb, rmn_bits, get_bits(gb, rmn_bits));
1913     avpriv_copy_bits(pb, data + size - rmn_bytes,
1914                  FFMIN(nbits - rmn_bits, rmn_bytes << 3));
1915 }
1916
1917 /**
1918  * Packet decoding: a packet is anything that the (ASF) demuxer contains,
1919  * and we expect that the demuxer / application provides it to us as such
1920  * (else you'll probably get garbage as output). Every packet has a size of
1921  * ctx->block_align bytes, starts with a packet header (see
1922  * #parse_packet_header()), and then a series of superframes. Superframe
1923  * boundaries may exceed packets, i.e. superframes can split data over
1924  * multiple (two) packets.
1925  *
1926  * For more information about frames, see #synth_superframe().
1927  */
1928 static int wmavoice_decode_packet(AVCodecContext *ctx, void *data,
1929                                   int *got_frame_ptr, AVPacket *avpkt)
1930 {
1931     WMAVoiceContext *s = ctx->priv_data;
1932     GetBitContext *gb = &s->gb;
1933     int size, res, pos;
1934
1935     /* Packets are sometimes a multiple of ctx->block_align, with a packet
1936      * header at each ctx->block_align bytes. However, Libav's ASF demuxer
1937      * feeds us ASF packets, which may concatenate multiple "codec" packets
1938      * in a single "muxer" packet, so we artificially emulate that by
1939      * capping the packet size at ctx->block_align. */
1940     for (size = avpkt->size; size > ctx->block_align; size -= ctx->block_align);
1941     if (!size) {
1942         *got_frame_ptr = 0;
1943         return 0;
1944     }
1945     init_get_bits(&s->gb, avpkt->data, size << 3);
1946
1947     /* size == ctx->block_align is used to indicate whether we are dealing with
1948      * a new packet or a packet of which we already read the packet header
1949      * previously. */
1950     if (size == ctx->block_align) { // new packet header
1951         if ((res = parse_packet_header(s)) < 0)
1952             return res;
1953
1954         /* If the packet header specifies a s->spillover_nbits, then we want
1955          * to push out all data of the previous packet (+ spillover) before
1956          * continuing to parse new superframes in the current packet. */
1957         if (s->spillover_nbits > 0) {
1958             if (s->sframe_cache_size > 0) {
1959                 int cnt = get_bits_count(gb);
1960                 copy_bits(&s->pb, avpkt->data, size, gb, s->spillover_nbits);
1961                 flush_put_bits(&s->pb);
1962                 s->sframe_cache_size += s->spillover_nbits;
1963                 if ((res = synth_superframe(ctx, got_frame_ptr)) == 0 &&
1964                     *got_frame_ptr) {
1965                     cnt += s->spillover_nbits;
1966                     s->skip_bits_next = cnt & 7;
1967                     *(AVFrame *)data = s->frame;
1968                     return cnt >> 3;
1969                 } else
1970                     skip_bits_long (gb, s->spillover_nbits - cnt +
1971                                     get_bits_count(gb)); // resync
1972             } else
1973                 skip_bits_long(gb, s->spillover_nbits);  // resync
1974         }
1975     } else if (s->skip_bits_next)
1976         skip_bits(gb, s->skip_bits_next);
1977
1978     /* Try parsing superframes in current packet */
1979     s->sframe_cache_size = 0;
1980     s->skip_bits_next = 0;
1981     pos = get_bits_left(gb);
1982     if ((res = synth_superframe(ctx, got_frame_ptr)) < 0) {
1983         return res;
1984     } else if (*got_frame_ptr) {
1985         int cnt = get_bits_count(gb);
1986         s->skip_bits_next = cnt & 7;
1987         *(AVFrame *)data = s->frame;
1988         return cnt >> 3;
1989     } else if ((s->sframe_cache_size = pos) > 0) {
1990         /* rewind bit reader to start of last (incomplete) superframe... */
1991         init_get_bits(gb, avpkt->data, size << 3);
1992         skip_bits_long(gb, (size << 3) - pos);
1993         assert(get_bits_left(gb) == pos);
1994
1995         /* ...and cache it for spillover in next packet */
1996         init_put_bits(&s->pb, s->sframe_cache, SFRAME_CACHE_MAXSIZE);
1997         copy_bits(&s->pb, avpkt->data, size, gb, s->sframe_cache_size);
1998         // FIXME bad - just copy bytes as whole and add use the
1999         // skip_bits_next field
2000     }
2001
2002     return size;
2003 }
2004
2005 static av_cold int wmavoice_decode_end(AVCodecContext *ctx)
2006 {
2007     WMAVoiceContext *s = ctx->priv_data;
2008
2009     if (s->do_apf) {
2010         ff_rdft_end(&s->rdft);
2011         ff_rdft_end(&s->irdft);
2012         ff_dct_end(&s->dct);
2013         ff_dct_end(&s->dst);
2014     }
2015
2016     return 0;
2017 }
2018
2019 static av_cold void wmavoice_flush(AVCodecContext *ctx)
2020 {
2021     WMAVoiceContext *s = ctx->priv_data;
2022     int n;
2023
2024     s->postfilter_agc    = 0;
2025     s->sframe_cache_size = 0;
2026     s->skip_bits_next    = 0;
2027     for (n = 0; n < s->lsps; n++)
2028         s->prev_lsps[n] = M_PI * (n + 1.0) / (s->lsps + 1.0);
2029     memset(s->excitation_history, 0,
2030            sizeof(*s->excitation_history) * MAX_SIGNAL_HISTORY);
2031     memset(s->synth_history,      0,
2032            sizeof(*s->synth_history)      * MAX_LSPS);
2033     memset(s->gain_pred_err,      0,
2034            sizeof(s->gain_pred_err));
2035
2036     if (s->do_apf) {
2037         memset(&s->synth_filter_out_buf[MAX_LSPS_ALIGN16 - s->lsps], 0,
2038                sizeof(*s->synth_filter_out_buf) * s->lsps);
2039         memset(s->dcf_mem,              0,
2040                sizeof(*s->dcf_mem)              * 2);
2041         memset(s->zero_exc_pf,          0,
2042                sizeof(*s->zero_exc_pf)          * s->history_nsamples);
2043         memset(s->denoise_filter_cache, 0, sizeof(s->denoise_filter_cache));
2044     }
2045 }
2046
2047 AVCodec ff_wmavoice_decoder = {
2048     .name           = "wmavoice",
2049     .type           = AVMEDIA_TYPE_AUDIO,
2050     .id             = AV_CODEC_ID_WMAVOICE,
2051     .priv_data_size = sizeof(WMAVoiceContext),
2052     .init           = wmavoice_decode_init,
2053     .close          = wmavoice_decode_end,
2054     .decode         = wmavoice_decode_packet,
2055     .capabilities   = CODEC_CAP_SUBFRAMES | CODEC_CAP_DR1,
2056     .flush          = wmavoice_flush,
2057     .long_name      = NULL_IF_CONFIG_SMALL("Windows Media Audio Voice"),
2058 };