Coverage Report

Created: 2023-03-26 06:13

/src/aac/libFDK/src/mdct.cpp
Line
Count
Source (jump to first uncovered line)
1
/* -----------------------------------------------------------------------------
2
Software License for The Fraunhofer FDK AAC Codec Library for Android
3
4
© Copyright  1995 - 2019 Fraunhofer-Gesellschaft zur Förderung der angewandten
5
Forschung e.V. All rights reserved.
6
7
 1.    INTRODUCTION
8
The Fraunhofer FDK AAC Codec Library for Android ("FDK AAC Codec") is software
9
that implements the MPEG Advanced Audio Coding ("AAC") encoding and decoding
10
scheme for digital audio. This FDK AAC Codec software is intended to be used on
11
a wide variety of Android devices.
12
13
AAC's HE-AAC and HE-AAC v2 versions are regarded as today's most efficient
14
general perceptual audio codecs. AAC-ELD is considered the best-performing
15
full-bandwidth communications codec by independent studies and is widely
16
deployed. AAC has been standardized by ISO and IEC as part of the MPEG
17
specifications.
18
19
Patent licenses for necessary patent claims for the FDK AAC Codec (including
20
those of Fraunhofer) may be obtained through Via Licensing
21
(www.vialicensing.com) or through the respective patent owners individually for
22
the purpose of encoding or decoding bit streams in products that are compliant
23
with the ISO/IEC MPEG audio standards. Please note that most manufacturers of
24
Android devices already license these patent claims through Via Licensing or
25
directly from the patent owners, and therefore FDK AAC Codec software may
26
already be covered under those patent licenses when it is used for those
27
licensed purposes only.
28
29
Commercially-licensed AAC software libraries, including floating-point versions
30
with enhanced sound quality, are also available from Fraunhofer. Users are
31
encouraged to check the Fraunhofer website for additional applications
32
information and documentation.
33
34
2.    COPYRIGHT LICENSE
35
36
Redistribution and use in source and binary forms, with or without modification,
37
are permitted without payment of copyright license fees provided that you
38
satisfy the following conditions:
39
40
You must retain the complete text of this software license in redistributions of
41
the FDK AAC Codec or your modifications thereto in source code form.
42
43
You must retain the complete text of this software license in the documentation
44
and/or other materials provided with redistributions of the FDK AAC Codec or
45
your modifications thereto in binary form. You must make available free of
46
charge copies of the complete source code of the FDK AAC Codec and your
47
modifications thereto to recipients of copies in binary form.
48
49
The name of Fraunhofer may not be used to endorse or promote products derived
50
from this library without prior written permission.
51
52
You may not charge copyright license fees for anyone to use, copy or distribute
53
the FDK AAC Codec software or your modifications thereto.
54
55
Your modified versions of the FDK AAC Codec must carry prominent notices stating
56
that you changed the software and the date of any change. For modified versions
57
of the FDK AAC Codec, the term "Fraunhofer FDK AAC Codec Library for Android"
58
must be replaced by the term "Third-Party Modified Version of the Fraunhofer FDK
59
AAC Codec Library for Android."
60
61
3.    NO PATENT LICENSE
62
63
NO EXPRESS OR IMPLIED LICENSES TO ANY PATENT CLAIMS, including without
64
limitation the patents of Fraunhofer, ARE GRANTED BY THIS SOFTWARE LICENSE.
65
Fraunhofer provides no warranty of patent non-infringement with respect to this
66
software.
67
68
You may use this FDK AAC Codec software or modifications thereto only for
69
purposes that are authorized by appropriate patent licenses.
70
71
4.    DISCLAIMER
72
73
This FDK AAC Codec software is provided by Fraunhofer on behalf of the copyright
74
holders and contributors "AS IS" and WITHOUT ANY EXPRESS OR IMPLIED WARRANTIES,
75
including but not limited to the implied warranties of merchantability and
76
fitness for a particular purpose. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
77
CONTRIBUTORS BE LIABLE for any direct, indirect, incidental, special, exemplary,
78
or consequential damages, including but not limited to procurement of substitute
79
goods or services; loss of use, data, or profits, or business interruption,
80
however caused and on any theory of liability, whether in contract, strict
81
liability, or tort (including negligence), arising in any way out of the use of
82
this software, even if advised of the possibility of such damage.
83
84
5.    CONTACT INFORMATION
85
86
Fraunhofer Institute for Integrated Circuits IIS
87
Attention: Audio and Multimedia Departments - FDK AAC LL
88
Am Wolfsmantel 33
89
91058 Erlangen, Germany
90
91
www.iis.fraunhofer.de/amm
92
amm-info@iis.fraunhofer.de
93
----------------------------------------------------------------------------- */
94
95
/******************* Library for basic calculation routines ********************
96
97
   Author(s):   Josef Hoepfl, Manuel Jander, Youliy Ninov, Daniel Hagel
98
99
   Description: MDCT/MDST routines
100
101
*******************************************************************************/
102
103
#include "mdct.h"
104
105
#include "FDK_tools_rom.h"
106
#include "dct.h"
107
#include "fixpoint_math.h"
108
109
261k
void mdct_init(H_MDCT hMdct, FIXP_DBL *overlap, INT overlapBufferSize) {
110
261k
  hMdct->overlap.freq = overlap;
111
  // FDKmemclear(overlap, overlapBufferSize*sizeof(FIXP_DBL));
112
261k
  hMdct->prev_fr = 0;
113
261k
  hMdct->prev_nr = 0;
114
261k
  hMdct->prev_tl = 0;
115
261k
  hMdct->ov_size = overlapBufferSize;
116
261k
  hMdct->prevAliasSymmetry = 0;
117
261k
  hMdct->prevPrevAliasSymmetry = 0;
118
261k
  hMdct->pFacZir = NULL;
119
261k
  hMdct->pAsymOvlp = NULL;
120
261k
}
121
122
/*
123
This program implements the forward MDCT transform on an input block of data.
124
The input block is in a form (A,B,C,D) where A,B,C and D are the respective
125
1/4th segments of the block. The program takes the input block and folds it in
126
the form:
127
(-D-Cr,A-Br). This block is twice shorter and here the 'r' suffix denotes
128
flipping of the sequence (reversing the order of the samples). While folding the
129
input block in the above mentioned shorter block the program windows the data.
130
Because the two operations (windowing and folding) are not implemented
131
sequentially, but together the program's structure is not easy to understand.
132
Once the output (already windowed) block (-D-Cr,A-Br) is ready it is passed to
133
the DCT IV for processing.
134
*/
135
INT mdct_block(H_MDCT hMdct, const INT_PCM *RESTRICT timeData,
136
               const INT noInSamples, FIXP_DBL *RESTRICT mdctData,
137
               const INT nSpec, const INT tl, const FIXP_WTP *pRightWindowPart,
138
0
               const INT fr, SHORT *pMdctData_e) {
139
0
  int i, n;
140
  /* tl: transform length
141
     fl: left window slope length
142
     nl: left window slope offset
143
     fr: right window slope length
144
     nr: right window slope offset
145
     See FDK_tools/doc/intern/mdct.tex for more detail. */
146
0
  int fl, nl, nr;
147
0
  const FIXP_WTP *wls, *wrs;
148
149
0
  wrs = pRightWindowPart;
150
151
  /* Detect FRprevious / FL mismatches and override parameters accordingly */
152
0
  if (hMdct->prev_fr ==
153
0
      0) { /* At start just initialize and pass parameters as they are */
154
0
    hMdct->prev_fr = fr;
155
0
    hMdct->prev_wrs = wrs;
156
0
    hMdct->prev_tl = tl;
157
0
  }
158
159
  /* Derive NR */
160
0
  nr = (tl - fr) >> 1;
161
162
  /* Skip input samples if tl is smaller than block size */
163
0
  timeData += (noInSamples - tl) >> 1;
164
165
  /* windowing */
166
0
  for (n = 0; n < nSpec; n++) {
167
    /*
168
     * MDCT scale:
169
     * + 1: fMultDiv2() in windowing.
170
     * + 1: Because of factor 1/2 in Princen-Bradley compliant windowed TDAC.
171
     */
172
0
    INT mdctData_e = 1 + 1;
173
174
    /* Derive left parameters */
175
0
    wls = hMdct->prev_wrs;
176
0
    fl = hMdct->prev_fr;
177
0
    nl = (tl - fl) >> 1;
178
179
    /* Here we implement a simplified version of what happens after the this
180
    piece of code (see the comments below). We implement the folding of A and B
181
    segments to (A-Br) but A is zero, because in this part of the MDCT sequence
182
    the window coefficients with which A must be multiplied are zero.    */
183
0
    for (i = 0; i < nl; i++) {
184
#if SAMPLE_BITS == DFRACT_BITS /* SPC_BITS and DFRACT_BITS should be equal. */
185
      mdctData[(tl / 2) + i] = -((FIXP_DBL)timeData[tl - i - 1] >> (1));
186
#else
187
0
      mdctData[(tl / 2) + i] = -(FIXP_DBL)timeData[tl - i - 1]
188
0
                               << (DFRACT_BITS - SAMPLE_BITS - 1); /* 0(A)-Br */
189
0
#endif
190
0
    }
191
192
    /* Implements the folding and windowing of the left part of the sequence,
193
    that is segments A and B. The A segment is multiplied by the respective left
194
    window coefficient and placed in a temporary variable.
195
196
    tmp0 = fMultDiv2((FIXP_PCM)timeData[i+nl], pLeftWindowPart[i].v.im);
197
198
    After this the B segment taken in reverse order is multiplied by the left
199
    window and subtracted from the previously derived temporary variable, so
200
    that finally we implement the A-Br operation. This output is written to the
201
    right part of the MDCT output : (-D-Cr,A-Br).
202
203
    mdctData[(tl/2)+i+nl] = fMultSubDiv2(tmp0, (FIXP_PCM)timeData[tl-nl-i-1],
204
    pLeftWindowPart[i].v.re);//A*window-Br*window
205
206
    The (A-Br) data is written to the output buffer (mdctData) without being
207
    flipped.     */
208
0
    for (i = 0; i < fl / 2; i++) {
209
0
      FIXP_DBL tmp0;
210
0
      tmp0 = fMultDiv2((FIXP_PCM)timeData[i + nl], wls[i].v.im); /* a*window */
211
0
      mdctData[(tl / 2) + i + nl] =
212
0
          fMultSubDiv2(tmp0, (FIXP_PCM)timeData[tl - nl - i - 1],
213
0
                       wls[i].v.re); /* A*window-Br*window */
214
0
    }
215
216
    /* Right window slope offset */
217
    /* Here we implement a simplified version of what happens after the this
218
    piece of code (see the comments below). We implement the folding of C and D
219
    segments to (-D-Cr) but D is zero, because in this part of the MDCT sequence
220
    the window coefficients with which D must be multiplied are zero.    */
221
0
    for (i = 0; i < nr; i++) {
222
#if SAMPLE_BITS == \
223
    DFRACT_BITS /* This should be SPC_BITS instead of DFRACT_BITS. */
224
      mdctData[(tl / 2) - 1 - i] = -((FIXP_DBL)timeData[tl + i] >> (1));
225
#else
226
0
      mdctData[(tl / 2) - 1 - i] =
227
0
          -(FIXP_DBL)timeData[tl + i]
228
0
          << (DFRACT_BITS - SAMPLE_BITS - 1); /* -C flipped at placing */
229
0
#endif
230
0
    }
231
232
    /* Implements the folding and windowing of the right part of the sequence,
233
    that is, segments C and D. The C segment is multiplied by the respective
234
    right window coefficient and placed in a temporary variable.
235
236
    tmp1 = fMultDiv2((FIXP_PCM)timeData[tl+nr+i], pRightWindowPart[i].v.re);
237
238
    After this the D segment taken in reverse order is multiplied by the right
239
    window and added from the previously derived temporary variable, so that we
240
    get (C+Dr) operation. This output is negated to get (-C-Dr) and written to
241
    the left part of the MDCT output while being reversed (flipped) at the same
242
    time, so that from (-C-Dr) we get (-D-Cr)=> (-D-Cr,A-Br).
243
244
    mdctData[(tl/2)-nr-i-1] = -fMultAddDiv2(tmp1,
245
    (FIXP_PCM)timeData[(tl*2)-nr-i-1], pRightWindowPart[i].v.im);*/
246
0
    for (i = 0; i < fr / 2; i++) {
247
0
      FIXP_DBL tmp1;
248
0
      tmp1 = fMultDiv2((FIXP_PCM)timeData[tl + nr + i],
249
0
                       wrs[i].v.re); /* C*window */
250
0
      mdctData[(tl / 2) - nr - i - 1] =
251
0
          -fMultAddDiv2(tmp1, (FIXP_PCM)timeData[(tl * 2) - nr - i - 1],
252
0
                        wrs[i].v.im); /* -(C*window+Dr*window) and flip before
253
                                         placing -> -Cr - D */
254
0
    }
255
256
    /* We pass the shortened folded data (-D-Cr,A-Br) to the MDCT function */
257
0
    dct_IV(mdctData, tl, &mdctData_e);
258
259
0
    pMdctData_e[n] = (SHORT)mdctData_e;
260
261
0
    timeData += tl;
262
0
    mdctData += tl;
263
264
0
    hMdct->prev_wrs = wrs;
265
0
    hMdct->prev_fr = fr;
266
0
    hMdct->prev_tl = tl;
267
0
  }
268
269
0
  return nSpec * tl;
270
0
}
271
272
1.00M
void imdct_gain(FIXP_DBL *pGain_m, int *pGain_e, int tl) {
273
1.00M
  FIXP_DBL gain_m = *pGain_m;
274
1.00M
  int gain_e = *pGain_e;
275
1.00M
  int log2_tl;
276
277
1.00M
  gain_e += -MDCT_OUTPUT_GAIN - MDCT_OUT_HEADROOM + 1;
278
1.00M
  if (tl == 0) {
279
    /* Dont regard the 2/N factor from the IDCT. It is compensated for somewhere
280
     * else. */
281
36.6k
    *pGain_e = gain_e;
282
36.6k
    return;
283
36.6k
  }
284
285
971k
  log2_tl = DFRACT_BITS - 1 - fNormz((FIXP_DBL)tl);
286
971k
  gain_e += -log2_tl;
287
288
  /* Detect non-radix 2 transform length and add amplitude compensation factor
289
     which cannot be included into the exponent above */
290
971k
  switch ((tl) >> (log2_tl - 2)) {
291
280k
    case 0x7: /* 10 ms, 1/tl = 1.0/(FDKpow(2.0, -log2_tl) *
292
                 0.53333333333333333333) */
293
280k
      if (gain_m == (FIXP_DBL)0) {
294
280k
        gain_m = FL2FXCONST_DBL(0.53333333333333333333f);
295
280k
      } else {
296
0
        gain_m = fMult(gain_m, FL2FXCONST_DBL(0.53333333333333333333f));
297
0
      }
298
280k
      break;
299
238k
    case 0x6: /* 3/4 of radix 2, 1/tl = 1.0/(FDKpow(2.0, -log2_tl) * 2.0/3.0) */
300
238k
      if (gain_m == (FIXP_DBL)0) {
301
206k
        gain_m = FL2FXCONST_DBL(2.0 / 3.0f);
302
206k
      } else {
303
32.2k
        gain_m = fMult(gain_m, FL2FXCONST_DBL(2.0 / 3.0f));
304
32.2k
      }
305
238k
      break;
306
199
    case 0x5: /* 0.8 of radix 2 (e.g. tl 160), 1/tl = 1.0/(FDKpow(2.0, -log2_tl)
307
               * 0.8/1.5) */
308
199
      if (gain_m == (FIXP_DBL)0) {
309
199
        gain_m = FL2FXCONST_DBL(0.53333333333333333333f);
310
199
      } else {
311
0
        gain_m = fMult(gain_m, FL2FXCONST_DBL(0.53333333333333333333f));
312
0
      }
313
199
      break;
314
452k
    case 0x4:
315
      /* radix 2, nothing to do. */
316
452k
      break;
317
0
    default:
318
      /* unsupported */
319
0
      FDK_ASSERT(0);
320
0
      break;
321
971k
  }
322
323
971k
  *pGain_m = gain_m;
324
971k
  *pGain_e = gain_e;
325
971k
}
326
327
105k
INT imdct_drain(H_MDCT hMdct, FIXP_DBL *output, INT nrSamplesRoom) {
328
105k
  int buffered_samples = 0;
329
330
105k
  if (nrSamplesRoom > 0) {
331
45.4k
    buffered_samples = hMdct->ov_offset;
332
333
45.4k
    FDK_ASSERT(buffered_samples <= nrSamplesRoom);
334
335
45.4k
    if (buffered_samples > 0) {
336
12.1k
      FDKmemcpy(output, hMdct->overlap.time,
337
12.1k
                buffered_samples * sizeof(FIXP_DBL));
338
12.1k
      hMdct->ov_offset = 0;
339
12.1k
    }
340
45.4k
  }
341
0
  return buffered_samples;
342
105k
}
343
344
72.1k
INT imdct_copy_ov_and_nr(H_MDCT hMdct, FIXP_DBL *pTimeData, INT nrSamples) {
345
72.1k
  FIXP_DBL *pOvl;
346
72.1k
  int nt, nf, i;
347
348
72.1k
  nt = fMin(hMdct->ov_offset, nrSamples);
349
72.1k
  nrSamples -= nt;
350
72.1k
  nf = fMin(hMdct->prev_nr, nrSamples);
351
72.1k
  FDKmemcpy(pTimeData, hMdct->overlap.time, nt * sizeof(FIXP_DBL));
352
72.1k
  pTimeData += nt;
353
354
72.1k
  pOvl = hMdct->overlap.freq + hMdct->ov_size - 1;
355
72.1k
  if (hMdct->prevPrevAliasSymmetry == 0) {
356
438k
    for (i = 0; i < nf; i++) {
357
365k
      FIXP_DBL x = -(*pOvl--);
358
365k
      *pTimeData = IMDCT_SCALE_DBL(x);
359
365k
      pTimeData++;
360
365k
    }
361
72.1k
  } else {
362
0
    for (i = 0; i < nf; i++) {
363
0
      FIXP_DBL x = (*pOvl--);
364
0
      *pTimeData = IMDCT_SCALE_DBL(x);
365
0
      pTimeData++;
366
0
    }
367
0
  }
368
369
72.1k
  return (nt + nf);
370
72.1k
}
371
372
void imdct_adapt_parameters(H_MDCT hMdct, int *pfl, int *pnl, int tl,
373
142k
                            const FIXP_WTP *wls, int noOutSamples) {
374
142k
  int fl = *pfl, nl = *pnl;
375
142k
  int window_diff, use_current = 0, use_previous = 0;
376
142k
  if (hMdct->prev_tl == 0) {
377
38.7k
    hMdct->prev_wrs = wls;
378
38.7k
    hMdct->prev_fr = fl;
379
38.7k
    hMdct->prev_nr = (noOutSamples - fl) >> 1;
380
38.7k
    hMdct->prev_tl = noOutSamples;
381
38.7k
    hMdct->ov_offset = 0;
382
38.7k
    use_current = 1;
383
38.7k
  }
384
385
142k
  window_diff = (hMdct->prev_fr - fl) >> 1;
386
387
  /* check if the previous window slope can be adjusted to match the current
388
   * window slope */
389
142k
  if (hMdct->prev_nr + window_diff > 0) {
390
58.3k
    use_current = 1;
391
58.3k
  }
392
  /* check if the current window slope can be adjusted to match the previous
393
   * window slope */
394
142k
  if (nl - window_diff > 0) {
395
66.4k
    use_previous = 1;
396
66.4k
  }
397
398
  /* if both is possible choose the larger of both window slope lengths */
399
142k
  if (use_current && use_previous) {
400
5.63k
    if (fl < hMdct->prev_fr) {
401
344
      use_current = 0;
402
344
    }
403
5.63k
  }
404
  /*
405
   * If the previous transform block is big enough, enlarge previous window
406
   * overlap, if not, then shrink current window overlap.
407
   */
408
142k
  if (use_current) {
409
81.4k
    hMdct->prev_nr += window_diff;
410
81.4k
    hMdct->prev_fr = fl;
411
81.4k
    hMdct->prev_wrs = wls;
412
81.4k
  } else {
413
61.1k
    nl -= window_diff;
414
61.1k
    fl = hMdct->prev_fr;
415
61.1k
  }
416
417
142k
  *pfl = fl;
418
142k
  *pnl = nl;
419
142k
}
420
421
/*
422
This program implements the inverse modulated lapped transform, a generalized
423
version of the inverse MDCT transform. Setting none of the MLT_*_ALIAS_FLAG
424
flags computes the IMDCT, setting all of them computes the IMDST. Other
425
combinations of these flags compute type III transforms used by the RSVD60
426
multichannel tool for transitions between MDCT/MDST. The following description
427
relates to the IMDCT only.
428
429
If we pass the data block (A,B,C,D,E,F) to the FORWARD MDCT it will produce two
430
outputs. The first one will be over the (A,B,C,D) part =>(-D-Cr,A-Br) and the
431
second one will be over the (C,D,E,F) part => (-F-Er,C-Dr), since there is a
432
overlap between consequtive passes of the algorithm. This overlap is over the
433
(C,D) segments. The two outputs will be given sequentially to the DCT IV
434
algorithm. At the INVERSE MDCT side we get two consecutive outputs from the IDCT
435
IV algorithm, namely the same blocks: (-D-Cr,A-Br) and (-F-Er,C-Dr). The first
436
of them lands in the Overlap buffer and the second is in the working one, which,
437
one algorithm pass later will substitute the one residing in the overlap
438
register. The IMDCT algorithm has to produce the C and D segments from the two
439
buffers. In order to do this we take the left part of the overlap
440
buffer(-D-Cr,A-Br), namely (-D-Cr) and add it appropriately to the right part of
441
the working buffer (-F-Er,C-Dr), namely (C-Dr), so that we get first the C
442
segment and later the D segment. We do this in the following way: From the right
443
part of the working buffer(C-Dr) we subtract the flipped left part of the
444
overlap buffer(-D-Cr):
445
446
Result = (C-Dr) - flipped(-D-Cr) = C -Dr + Dr + C = 2C
447
We divide by two and get the C segment. What we did is adding the right part of
448
the first frame to the left part of the second one.   While applying these
449
operation we multiply the respective segments with the appropriate window
450
functions.
451
452
In order to get the D segment we do the following:
453
From the negated second part of the working buffer(C-Dr) we subtract the flipped
454
first part of the overlap buffer (-D-Cr):
455
456
Result= - (C -Dr) - flipped(-D-Cr)= -C +Dr +Dr +C = 2Dr.
457
After dividing by two and flipping we get the D segment.What we did is adding
458
the right part of the first frame to the left part of the second one.   While
459
applying these operation we multiply the respective segments with the
460
appropriate window functions.
461
462
Once we have obtained the C and D segments the overlap buffer is emptied and the
463
current buffer is sent in it, so that the E and F segments are available for
464
decoding in the next algorithm pass.*/
465
INT imlt_block(H_MDCT hMdct, FIXP_DBL *output, FIXP_DBL *spectrum,
466
               const SHORT scalefactor[], const INT nSpec,
467
               const INT noOutSamples, const INT tl, const FIXP_WTP *wls,
468
               INT fl, const FIXP_WTP *wrs, const INT fr, FIXP_DBL gain,
469
477k
               int flags) {
470
477k
  FIXP_DBL *pOvl;
471
477k
  FIXP_DBL *pOut0 = output, *pOut1;
472
477k
  INT nl, nr;
473
477k
  int w, i, nrSamples = 0, specShiftScale, transform_gain_e = 0;
474
477k
  int currAliasSymmetry = (flags & MLT_FLAG_CURR_ALIAS_SYMMETRY);
475
476
  /* Derive NR and NL */
477
477k
  nr = (tl - fr) >> 1;
478
477k
  nl = (tl - fl) >> 1;
479
480
  /* Include 2/N IMDCT gain into gain factor and exponent. */
481
477k
  imdct_gain(&gain, &transform_gain_e, tl);
482
483
  /* Detect FRprevious / FL mismatches and override parameters accordingly */
484
477k
  if (hMdct->prev_fr != fl) {
485
128k
    imdct_adapt_parameters(hMdct, &fl, &nl, tl, wls, noOutSamples);
486
128k
  }
487
488
477k
  pOvl = hMdct->overlap.freq + hMdct->ov_size - 1;
489
490
477k
  if (noOutSamples > nrSamples) {
491
    /* Purge buffered output. */
492
74.2M
    for (i = 0; i < hMdct->ov_offset; i++) {
493
73.7M
      *pOut0 = hMdct->overlap.time[i];
494
73.7M
      pOut0++;
495
73.7M
    }
496
473k
    nrSamples = hMdct->ov_offset;
497
473k
    hMdct->ov_offset = 0;
498
473k
  }
499
500
2.17M
  for (w = 0; w < nSpec; w++) {
501
1.69M
    FIXP_DBL *pSpec, *pCurr;
502
1.69M
    const FIXP_WTP *pWindow;
503
504
    /* Detect FRprevious / FL mismatches and override parameters accordingly */
505
1.69M
    if (hMdct->prev_fr != fl) {
506
0
      imdct_adapt_parameters(hMdct, &fl, &nl, tl, wls, noOutSamples);
507
0
    }
508
509
1.69M
    specShiftScale = transform_gain_e;
510
511
    /* Setup window pointers */
512
1.69M
    pWindow = hMdct->prev_wrs;
513
514
    /* Current spectrum */
515
1.69M
    pSpec = spectrum + w * tl;
516
517
    /* DCT IV of current spectrum. */
518
1.69M
    if (currAliasSymmetry == 0) {
519
1.69M
      if (hMdct->prevAliasSymmetry == 0) {
520
1.69M
        dct_IV(pSpec, tl, &specShiftScale);
521
1.69M
      } else {
522
0
        FIXP_DBL _tmp[1024 + ALIGNMENT_DEFAULT / sizeof(FIXP_DBL)];
523
0
        FIXP_DBL *tmp = (FIXP_DBL *)ALIGN_PTR(_tmp);
524
0
        C_ALLOC_ALIGNED_REGISTER(tmp, sizeof(_tmp));
525
0
        dct_III(pSpec, tmp, tl, &specShiftScale);
526
0
        C_ALLOC_ALIGNED_UNREGISTER(tmp);
527
0
      }
528
1.69M
    } else {
529
0
      if (hMdct->prevAliasSymmetry == 0) {
530
0
        FIXP_DBL _tmp[1024 + ALIGNMENT_DEFAULT / sizeof(FIXP_DBL)];
531
0
        FIXP_DBL *tmp = (FIXP_DBL *)ALIGN_PTR(_tmp);
532
0
        C_ALLOC_ALIGNED_REGISTER(tmp, sizeof(_tmp));
533
0
        dst_III(pSpec, tmp, tl, &specShiftScale);
534
0
        C_ALLOC_ALIGNED_UNREGISTER(tmp);
535
0
      } else {
536
0
        dst_IV(pSpec, tl, &specShiftScale);
537
0
      }
538
0
    }
539
540
    /* Optional scaling of time domain - no yet windowed - of current spectrum
541
     */
542
    /* and de-scale current spectrum signal (time domain, no yet windowed) */
543
1.69M
    if (gain != (FIXP_DBL)0) {
544
168M
      for (i = 0; i < tl; i++) {
545
167M
        pSpec[i] = fMult(pSpec[i], gain);
546
167M
      }
547
939k
    }
548
549
1.69M
    {
550
1.69M
      int loc_scale =
551
1.69M
          fixmin_I(scalefactor[w] + specShiftScale, (INT)DFRACT_BITS - 1);
552
1.69M
      DWORD_ALIGNED(pSpec);
553
1.69M
      scaleValuesSaturate(pSpec, tl, loc_scale);
554
1.69M
    }
555
556
1.69M
    if (noOutSamples <= nrSamples) {
557
      /* Divert output first half to overlap buffer if we already got enough
558
       * output samples. */
559
527k
      pOut0 = hMdct->overlap.time + hMdct->ov_offset;
560
527k
      hMdct->ov_offset += hMdct->prev_nr + fl / 2;
561
1.17M
    } else {
562
      /* Account output samples */
563
1.17M
      nrSamples += hMdct->prev_nr + fl / 2;
564
1.17M
    }
565
566
    /* NR output samples 0 .. NR. -overlap[TL/2..TL/2-NR] */
567
1.69M
    if ((hMdct->pFacZir != 0) && (hMdct->prev_nr == fl / 2)) {
568
      /* In the case of ACELP -> TCX20 -> FD short add FAC ZIR on nr signal part
569
       */
570
69.0k
      for (i = 0; i < hMdct->prev_nr; i++) {
571
67.8k
        FIXP_DBL x = -(*pOvl--);
572
67.8k
        *pOut0 = fAddSaturate(x, IMDCT_SCALE_DBL(hMdct->pFacZir[i]));
573
67.8k
        pOut0++;
574
67.8k
      }
575
1.25k
      hMdct->pFacZir = NULL;
576
1.69M
    } else {
577
      /* Here we implement a simplified version of what happens after the this
578
      piece of code (see the comments below). We implement the folding of C and
579
      D segments from (-D-Cr) but D is zero, because in this part of the MDCT
580
      sequence the window coefficients with which D must be multiplied are zero.
581
      "pOut0" writes sequentially the C block from left to right.   */
582
1.69M
      if (hMdct->prevPrevAliasSymmetry == 0) {
583
34.6M
        for (i = 0; i < hMdct->prev_nr; i++) {
584
32.9M
          FIXP_DBL x = -(*pOvl--);
585
32.9M
          *pOut0 = IMDCT_SCALE_DBL(x);
586
32.9M
          pOut0++;
587
32.9M
        }
588
1.69M
      } else {
589
0
        for (i = 0; i < hMdct->prev_nr; i++) {
590
0
          FIXP_DBL x = *pOvl--;
591
0
          *pOut0 = IMDCT_SCALE_DBL(x);
592
0
          pOut0++;
593
0
        }
594
0
      }
595
1.69M
    }
596
597
1.69M
    if (noOutSamples <= nrSamples) {
598
      /* Divert output second half to overlap buffer if we already got enough
599
       * output samples. */
600
710k
      pOut1 = hMdct->overlap.time + hMdct->ov_offset + fl / 2 - 1;
601
710k
      hMdct->ov_offset += fl / 2 + nl;
602
989k
    } else {
603
989k
      pOut1 = pOut0 + (fl - 1);
604
989k
      nrSamples += fl / 2 + nl;
605
989k
    }
606
607
    /* output samples before window crossing point NR .. TL/2.
608
     * -overlap[TL/2-NR..TL/2-NR-FL/2] + current[NR..TL/2] */
609
    /* output samples after window crossing point TL/2 .. TL/2+FL/2.
610
     * -overlap[0..FL/2] - current[TL/2..FL/2] */
611
1.69M
    pCurr = pSpec + tl - fl / 2;
612
1.69M
    DWORD_ALIGNED(pCurr);
613
1.69M
    C_ALLOC_ALIGNED_REGISTER(pWindow, fl);
614
1.69M
    DWORD_ALIGNED(pWindow);
615
1.69M
    C_ALLOC_ALIGNED_UNREGISTER(pWindow);
616
617
1.69M
    if (hMdct->prevPrevAliasSymmetry == 0) {
618
1.69M
      if (hMdct->prevAliasSymmetry == 0) {
619
1.69M
        if (!hMdct->pAsymOvlp) {
620
181M
          for (i = 0; i < fl / 2; i++) {
621
180M
            FIXP_DBL x0, x1;
622
180M
            cplxMultDiv2(&x1, &x0, *pCurr++, -*pOvl--, pWindow[i]);
623
180M
            *pOut0 = IMDCT_SCALE_DBL_LSH1(x0);
624
180M
            *pOut1 = IMDCT_SCALE_DBL_LSH1(-x1);
625
180M
            pOut0++;
626
180M
            pOut1--;
627
180M
          }
628
1.69M
        } else {
629
0
          FIXP_DBL *pAsymOvl = hMdct->pAsymOvlp + fl / 2 - 1;
630
0
          for (i = 0; i < fl / 2; i++) {
631
0
            FIXP_DBL x0, x1;
632
0
            x1 = -fMultDiv2(*pCurr, pWindow[i].v.re) +
633
0
                 fMultDiv2(*pAsymOvl, pWindow[i].v.im);
634
0
            x0 = fMultDiv2(*pCurr, pWindow[i].v.im) -
635
0
                 fMultDiv2(*pOvl, pWindow[i].v.re);
636
0
            pCurr++;
637
0
            pOvl--;
638
0
            pAsymOvl--;
639
0
            *pOut0++ = IMDCT_SCALE_DBL_LSH1(x0);
640
0
            *pOut1-- = IMDCT_SCALE_DBL_LSH1(x1);
641
0
          }
642
0
          hMdct->pAsymOvlp = NULL;
643
0
        }
644
1.69M
      } else { /* prevAliasingSymmetry == 1 */
645
0
        for (i = 0; i < fl / 2; i++) {
646
0
          FIXP_DBL x0, x1;
647
0
          cplxMultDiv2(&x1, &x0, *pCurr++, -*pOvl--, pWindow[i]);
648
0
          *pOut0 = IMDCT_SCALE_DBL_LSH1(x0);
649
0
          *pOut1 = IMDCT_SCALE_DBL_LSH1(x1);
650
0
          pOut0++;
651
0
          pOut1--;
652
0
        }
653
0
      }
654
1.69M
    } else { /* prevPrevAliasingSymmetry == 1 */
655
0
      if (hMdct->prevAliasSymmetry == 0) {
656
0
        for (i = 0; i < fl / 2; i++) {
657
0
          FIXP_DBL x0, x1;
658
0
          cplxMultDiv2(&x1, &x0, *pCurr++, *pOvl--, pWindow[i]);
659
0
          *pOut0 = IMDCT_SCALE_DBL_LSH1(x0);
660
0
          *pOut1 = IMDCT_SCALE_DBL_LSH1(-x1);
661
0
          pOut0++;
662
0
          pOut1--;
663
0
        }
664
0
      } else { /* prevAliasingSymmetry == 1 */
665
0
        for (i = 0; i < fl / 2; i++) {
666
0
          FIXP_DBL x0, x1;
667
0
          cplxMultDiv2(&x1, &x0, *pCurr++, *pOvl--, pWindow[i]);
668
0
          *pOut0 = IMDCT_SCALE_DBL_LSH1(x0);
669
0
          *pOut1 = IMDCT_SCALE_DBL_LSH1(x1);
670
0
          pOut0++;
671
0
          pOut1--;
672
0
        }
673
0
      }
674
0
    }
675
676
1.69M
    if (hMdct->pFacZir != 0) {
677
      /* add FAC ZIR of previous ACELP -> mdct transition */
678
6.89k
      FIXP_DBL *pOut = pOut0 - fl / 2;
679
6.89k
      FDK_ASSERT(fl / 2 <= 128);
680
853k
      for (i = 0; i < fl / 2; i++) {
681
846k
        pOut[i] = fAddSaturate(pOut[i], IMDCT_SCALE_DBL(hMdct->pFacZir[i]));
682
846k
      }
683
6.89k
      hMdct->pFacZir = NULL;
684
6.89k
    }
685
0
    pOut0 += (fl / 2) + nl;
686
687
    /* NL output samples TL/2+FL/2..TL. - current[FL/2..0] */
688
1.69M
    pOut1 += (fl / 2) + 1;
689
1.69M
    pCurr = pSpec + tl - fl / 2 - 1;
690
    /* Here we implement a simplified version of what happens above the this
691
    piece of code (see the comments above). We implement the folding of C and D
692
    segments from (C-Dr) but C is zero, because in this part of the MDCT
693
    sequence the window coefficients with which C must be multiplied are zero.
694
    "pOut1" writes sequentially the D block from left to right.   */
695
1.69M
    if (hMdct->prevAliasSymmetry == 0) {
696
34.2M
      for (i = 0; i < nl; i++) {
697
32.5M
        FIXP_DBL x = -(*pCurr--);
698
32.5M
        *pOut1++ = IMDCT_SCALE_DBL(x);
699
32.5M
      }
700
1.69M
    } else {
701
0
      for (i = 0; i < nl; i++) {
702
0
        FIXP_DBL x = *pCurr--;
703
0
        *pOut1++ = IMDCT_SCALE_DBL(x);
704
0
      }
705
0
    }
706
707
    /* Set overlap source pointer for next window pOvl = pSpec + tl/2 - 1; */
708
1.69M
    pOvl = pSpec + tl / 2 - 1;
709
710
    /* Previous window values. */
711
1.69M
    hMdct->prev_nr = nr;
712
1.69M
    hMdct->prev_fr = fr;
713
1.69M
    hMdct->prev_tl = tl;
714
1.69M
    hMdct->prev_wrs = wrs;
715
716
    /* Previous aliasing symmetry */
717
1.69M
    hMdct->prevPrevAliasSymmetry = hMdct->prevAliasSymmetry;
718
1.69M
    hMdct->prevAliasSymmetry = currAliasSymmetry;
719
1.69M
  }
720
721
  /* Save overlap */
722
723
477k
  pOvl = hMdct->overlap.freq + hMdct->ov_size - tl / 2;
724
477k
  FDKmemcpy(pOvl, &spectrum[(nSpec - 1) * tl], (tl / 2) * sizeof(FIXP_DBL));
725
726
477k
  return nrSamples;
727
477k
}