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