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