Subversion Repositories Kolibri OS

Rev

Blame | Last modification | View Log | RSS feed

  1. /*
  2.  * AAC Spectral Band Replication decoding functions
  3.  * Copyright (c) 2008-2009 Robert Swain ( rob opendot cl )
  4.  * Copyright (c) 2009-2010 Alex Converse <alex.converse@gmail.com>
  5.  *
  6.  * This file is part of FFmpeg.
  7.  *
  8.  * FFmpeg is free software; you can redistribute it and/or
  9.  * modify it under the terms of the GNU Lesser General Public
  10.  * License as published by the Free Software Foundation; either
  11.  * version 2.1 of the License, or (at your option) any later version.
  12.  *
  13.  * FFmpeg is distributed in the hope that it will be useful,
  14.  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  15.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  16.  * Lesser General Public License for more details.
  17.  *
  18.  * You should have received a copy of the GNU Lesser General Public
  19.  * License along with FFmpeg; if not, write to the Free Software
  20.  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  21.  */
  22.  
  23. /**
  24.  * @file
  25.  * AAC Spectral Band Replication decoding functions
  26.  * @author Robert Swain ( rob opendot cl )
  27.  */
  28. #define USE_FIXED 0
  29.  
  30. #include "aac.h"
  31. #include "sbr.h"
  32. #include "aacsbr.h"
  33. #include "aacsbrdata.h"
  34. #include "aacsbr_tablegen.h"
  35. #include "fft.h"
  36. #include "aacps.h"
  37. #include "sbrdsp.h"
  38. #include "libavutil/internal.h"
  39. #include "libavutil/libm.h"
  40. #include "libavutil/avassert.h"
  41.  
  42. #include <stdint.h>
  43. #include <float.h>
  44. #include <math.h>
  45.  
  46. #if ARCH_MIPS
  47. #include "mips/aacsbr_mips.h"
  48. #endif /* ARCH_MIPS */
  49.  
  50. static VLC vlc_sbr[10];
  51. static void aacsbr_func_ptr_init(AACSBRContext *c);
  52.  
  53. static void make_bands(int16_t* bands, int start, int stop, int num_bands)
  54. {
  55.     int k, previous, present;
  56.     float base, prod;
  57.  
  58.     base = powf((float)stop / start, 1.0f / num_bands);
  59.     prod = start;
  60.     previous = start;
  61.  
  62.     for (k = 0; k < num_bands-1; k++) {
  63.         prod *= base;
  64.         present  = lrintf(prod);
  65.         bands[k] = present - previous;
  66.         previous = present;
  67.     }
  68.     bands[num_bands-1] = stop - previous;
  69. }
  70.  
  71. /// Dequantization and stereo decoding (14496-3 sp04 p203)
  72. static void sbr_dequant(SpectralBandReplication *sbr, int id_aac)
  73. {
  74.     int k, e;
  75.     int ch;
  76.  
  77.     if (id_aac == TYPE_CPE && sbr->bs_coupling) {
  78.         float alpha      = sbr->data[0].bs_amp_res ?  1.0f :  0.5f;
  79.         float pan_offset = sbr->data[0].bs_amp_res ? 12.0f : 24.0f;
  80.         for (e = 1; e <= sbr->data[0].bs_num_env; e++) {
  81.             for (k = 0; k < sbr->n[sbr->data[0].bs_freq_res[e]]; k++) {
  82.                 float temp1 = exp2f(sbr->data[0].env_facs[e][k] * alpha + 7.0f);
  83.                 float temp2 = exp2f((pan_offset - sbr->data[1].env_facs[e][k]) * alpha);
  84.                 float fac;
  85.                 if (temp1 > 1E20) {
  86.                     av_log(NULL, AV_LOG_ERROR, "envelope scalefactor overflow in dequant\n");
  87.                     temp1 = 1;
  88.                 }
  89.                 fac   = temp1 / (1.0f + temp2);
  90.                 sbr->data[0].env_facs[e][k] = fac;
  91.                 sbr->data[1].env_facs[e][k] = fac * temp2;
  92.             }
  93.         }
  94.         for (e = 1; e <= sbr->data[0].bs_num_noise; e++) {
  95.             for (k = 0; k < sbr->n_q; k++) {
  96.                 float temp1 = exp2f(NOISE_FLOOR_OFFSET - sbr->data[0].noise_facs[e][k] + 1);
  97.                 float temp2 = exp2f(12 - sbr->data[1].noise_facs[e][k]);
  98.                 float fac;
  99.                 if (temp1 > 1E20) {
  100.                     av_log(NULL, AV_LOG_ERROR, "envelope scalefactor overflow in dequant\n");
  101.                     temp1 = 1;
  102.                 }
  103.                 fac = temp1 / (1.0f + temp2);
  104.                 sbr->data[0].noise_facs[e][k] = fac;
  105.                 sbr->data[1].noise_facs[e][k] = fac * temp2;
  106.             }
  107.         }
  108.     } else { // SCE or one non-coupled CPE
  109.         for (ch = 0; ch < (id_aac == TYPE_CPE) + 1; ch++) {
  110.             float alpha = sbr->data[ch].bs_amp_res ? 1.0f : 0.5f;
  111.             for (e = 1; e <= sbr->data[ch].bs_num_env; e++)
  112.                 for (k = 0; k < sbr->n[sbr->data[ch].bs_freq_res[e]]; k++){
  113.                     sbr->data[ch].env_facs[e][k] =
  114.                         exp2f(alpha * sbr->data[ch].env_facs[e][k] + 6.0f);
  115.                     if (sbr->data[ch].env_facs[e][k] > 1E20) {
  116.                         av_log(NULL, AV_LOG_ERROR, "envelope scalefactor overflow in dequant\n");
  117.                         sbr->data[ch].env_facs[e][k] = 1;
  118.                     }
  119.                 }
  120.  
  121.             for (e = 1; e <= sbr->data[ch].bs_num_noise; e++)
  122.                 for (k = 0; k < sbr->n_q; k++)
  123.                     sbr->data[ch].noise_facs[e][k] =
  124.                         exp2f(NOISE_FLOOR_OFFSET - sbr->data[ch].noise_facs[e][k]);
  125.         }
  126.     }
  127. }
  128.  
  129. /** High Frequency Generation (14496-3 sp04 p214+) and Inverse Filtering
  130.  * (14496-3 sp04 p214)
  131.  * Warning: This routine does not seem numerically stable.
  132.  */
  133. static void sbr_hf_inverse_filter(SBRDSPContext *dsp,
  134.                                   float (*alpha0)[2], float (*alpha1)[2],
  135.                                   const float X_low[32][40][2], int k0)
  136. {
  137.     int k;
  138.     for (k = 0; k < k0; k++) {
  139.         LOCAL_ALIGNED_16(float, phi, [3], [2][2]);
  140.         float dk;
  141.  
  142.         dsp->autocorrelate(X_low[k], phi);
  143.  
  144.         dk =  phi[2][1][0] * phi[1][0][0] -
  145.              (phi[1][1][0] * phi[1][1][0] + phi[1][1][1] * phi[1][1][1]) / 1.000001f;
  146.  
  147.         if (!dk) {
  148.             alpha1[k][0] = 0;
  149.             alpha1[k][1] = 0;
  150.         } else {
  151.             float temp_real, temp_im;
  152.             temp_real = phi[0][0][0] * phi[1][1][0] -
  153.                         phi[0][0][1] * phi[1][1][1] -
  154.                         phi[0][1][0] * phi[1][0][0];
  155.             temp_im   = phi[0][0][0] * phi[1][1][1] +
  156.                         phi[0][0][1] * phi[1][1][0] -
  157.                         phi[0][1][1] * phi[1][0][0];
  158.  
  159.             alpha1[k][0] = temp_real / dk;
  160.             alpha1[k][1] = temp_im   / dk;
  161.         }
  162.  
  163.         if (!phi[1][0][0]) {
  164.             alpha0[k][0] = 0;
  165.             alpha0[k][1] = 0;
  166.         } else {
  167.             float temp_real, temp_im;
  168.             temp_real = phi[0][0][0] + alpha1[k][0] * phi[1][1][0] +
  169.                                        alpha1[k][1] * phi[1][1][1];
  170.             temp_im   = phi[0][0][1] + alpha1[k][1] * phi[1][1][0] -
  171.                                        alpha1[k][0] * phi[1][1][1];
  172.  
  173.             alpha0[k][0] = -temp_real / phi[1][0][0];
  174.             alpha0[k][1] = -temp_im   / phi[1][0][0];
  175.         }
  176.  
  177.         if (alpha1[k][0] * alpha1[k][0] + alpha1[k][1] * alpha1[k][1] >= 16.0f ||
  178.            alpha0[k][0] * alpha0[k][0] + alpha0[k][1] * alpha0[k][1] >= 16.0f) {
  179.             alpha1[k][0] = 0;
  180.             alpha1[k][1] = 0;
  181.             alpha0[k][0] = 0;
  182.             alpha0[k][1] = 0;
  183.         }
  184.     }
  185. }
  186.  
  187. /// Chirp Factors (14496-3 sp04 p214)
  188. static void sbr_chirp(SpectralBandReplication *sbr, SBRData *ch_data)
  189. {
  190.     int i;
  191.     float new_bw;
  192.     static const float bw_tab[] = { 0.0f, 0.75f, 0.9f, 0.98f };
  193.  
  194.     for (i = 0; i < sbr->n_q; i++) {
  195.         if (ch_data->bs_invf_mode[0][i] + ch_data->bs_invf_mode[1][i] == 1) {
  196.             new_bw = 0.6f;
  197.         } else
  198.             new_bw = bw_tab[ch_data->bs_invf_mode[0][i]];
  199.  
  200.         if (new_bw < ch_data->bw_array[i]) {
  201.             new_bw = 0.75f    * new_bw + 0.25f    * ch_data->bw_array[i];
  202.         } else
  203.             new_bw = 0.90625f * new_bw + 0.09375f * ch_data->bw_array[i];
  204.         ch_data->bw_array[i] = new_bw < 0.015625f ? 0.0f : new_bw;
  205.     }
  206. }
  207.  
  208. /**
  209.  * Calculation of levels of additional HF signal components (14496-3 sp04 p219)
  210.  * and Calculation of gain (14496-3 sp04 p219)
  211.  */
  212. static void sbr_gain_calc(AACContext *ac, SpectralBandReplication *sbr,
  213.                           SBRData *ch_data, const int e_a[2])
  214. {
  215.     int e, k, m;
  216.     // max gain limits : -3dB, 0dB, 3dB, inf dB (limiter off)
  217.     static const float limgain[4] = { 0.70795, 1.0, 1.41254, 10000000000 };
  218.  
  219.     for (e = 0; e < ch_data->bs_num_env; e++) {
  220.         int delta = !((e == e_a[1]) || (e == e_a[0]));
  221.         for (k = 0; k < sbr->n_lim; k++) {
  222.             float gain_boost, gain_max;
  223.             float sum[2] = { 0.0f, 0.0f };
  224.             for (m = sbr->f_tablelim[k] - sbr->kx[1]; m < sbr->f_tablelim[k + 1] - sbr->kx[1]; m++) {
  225.                 const float temp = sbr->e_origmapped[e][m] / (1.0f + sbr->q_mapped[e][m]);
  226.                 sbr->q_m[e][m] = sqrtf(temp * sbr->q_mapped[e][m]);
  227.                 sbr->s_m[e][m] = sqrtf(temp * ch_data->s_indexmapped[e + 1][m]);
  228.                 if (!sbr->s_mapped[e][m]) {
  229.                     sbr->gain[e][m] = sqrtf(sbr->e_origmapped[e][m] /
  230.                                             ((1.0f + sbr->e_curr[e][m]) *
  231.                                              (1.0f + sbr->q_mapped[e][m] * delta)));
  232.                 } else {
  233.                     sbr->gain[e][m] = sqrtf(sbr->e_origmapped[e][m] * sbr->q_mapped[e][m] /
  234.                                             ((1.0f + sbr->e_curr[e][m]) *
  235.                                              (1.0f + sbr->q_mapped[e][m])));
  236.                 }
  237.             }
  238.             for (m = sbr->f_tablelim[k] - sbr->kx[1]; m < sbr->f_tablelim[k + 1] - sbr->kx[1]; m++) {
  239.                 sum[0] += sbr->e_origmapped[e][m];
  240.                 sum[1] += sbr->e_curr[e][m];
  241.             }
  242.             gain_max = limgain[sbr->bs_limiter_gains] * sqrtf((FLT_EPSILON + sum[0]) / (FLT_EPSILON + sum[1]));
  243.             gain_max = FFMIN(100000.f, gain_max);
  244.             for (m = sbr->f_tablelim[k] - sbr->kx[1]; m < sbr->f_tablelim[k + 1] - sbr->kx[1]; m++) {
  245.                 float q_m_max   = sbr->q_m[e][m] * gain_max / sbr->gain[e][m];
  246.                 sbr->q_m[e][m]  = FFMIN(sbr->q_m[e][m], q_m_max);
  247.                 sbr->gain[e][m] = FFMIN(sbr->gain[e][m], gain_max);
  248.             }
  249.             sum[0] = sum[1] = 0.0f;
  250.             for (m = sbr->f_tablelim[k] - sbr->kx[1]; m < sbr->f_tablelim[k + 1] - sbr->kx[1]; m++) {
  251.                 sum[0] += sbr->e_origmapped[e][m];
  252.                 sum[1] += sbr->e_curr[e][m] * sbr->gain[e][m] * sbr->gain[e][m]
  253.                           + sbr->s_m[e][m] * sbr->s_m[e][m]
  254.                           + (delta && !sbr->s_m[e][m]) * sbr->q_m[e][m] * sbr->q_m[e][m];
  255.             }
  256.             gain_boost = sqrtf((FLT_EPSILON + sum[0]) / (FLT_EPSILON + sum[1]));
  257.             gain_boost = FFMIN(1.584893192f, gain_boost);
  258.             for (m = sbr->f_tablelim[k] - sbr->kx[1]; m < sbr->f_tablelim[k + 1] - sbr->kx[1]; m++) {
  259.                 sbr->gain[e][m] *= gain_boost;
  260.                 sbr->q_m[e][m]  *= gain_boost;
  261.                 sbr->s_m[e][m]  *= gain_boost;
  262.             }
  263.         }
  264.     }
  265. }
  266.  
  267. /// Assembling HF Signals (14496-3 sp04 p220)
  268. static void sbr_hf_assemble(float Y1[38][64][2],
  269.                             const float X_high[64][40][2],
  270.                             SpectralBandReplication *sbr, SBRData *ch_data,
  271.                             const int e_a[2])
  272. {
  273.     int e, i, j, m;
  274.     const int h_SL = 4 * !sbr->bs_smoothing_mode;
  275.     const int kx = sbr->kx[1];
  276.     const int m_max = sbr->m[1];
  277.     static const float h_smooth[5] = {
  278.         0.33333333333333,
  279.         0.30150283239582,
  280.         0.21816949906249,
  281.         0.11516383427084,
  282.         0.03183050093751,
  283.     };
  284.     float (*g_temp)[48] = ch_data->g_temp, (*q_temp)[48] = ch_data->q_temp;
  285.     int indexnoise = ch_data->f_indexnoise;
  286.     int indexsine  = ch_data->f_indexsine;
  287.  
  288.     if (sbr->reset) {
  289.         for (i = 0; i < h_SL; i++) {
  290.             memcpy(g_temp[i + 2*ch_data->t_env[0]], sbr->gain[0], m_max * sizeof(sbr->gain[0][0]));
  291.             memcpy(q_temp[i + 2*ch_data->t_env[0]], sbr->q_m[0],  m_max * sizeof(sbr->q_m[0][0]));
  292.         }
  293.     } else if (h_SL) {
  294.         for (i = 0; i < 4; i++) {
  295.             memcpy(g_temp[i + 2 * ch_data->t_env[0]],
  296.                    g_temp[i + 2 * ch_data->t_env_num_env_old],
  297.                    sizeof(g_temp[0]));
  298.             memcpy(q_temp[i + 2 * ch_data->t_env[0]],
  299.                    q_temp[i + 2 * ch_data->t_env_num_env_old],
  300.                    sizeof(q_temp[0]));
  301.         }
  302.     }
  303.  
  304.     for (e = 0; e < ch_data->bs_num_env; e++) {
  305.         for (i = 2 * ch_data->t_env[e]; i < 2 * ch_data->t_env[e + 1]; i++) {
  306.             memcpy(g_temp[h_SL + i], sbr->gain[e], m_max * sizeof(sbr->gain[0][0]));
  307.             memcpy(q_temp[h_SL + i], sbr->q_m[e],  m_max * sizeof(sbr->q_m[0][0]));
  308.         }
  309.     }
  310.  
  311.     for (e = 0; e < ch_data->bs_num_env; e++) {
  312.         for (i = 2 * ch_data->t_env[e]; i < 2 * ch_data->t_env[e + 1]; i++) {
  313.             LOCAL_ALIGNED_16(float, g_filt_tab, [48]);
  314.             LOCAL_ALIGNED_16(float, q_filt_tab, [48]);
  315.             float *g_filt, *q_filt;
  316.  
  317.             if (h_SL && e != e_a[0] && e != e_a[1]) {
  318.                 g_filt = g_filt_tab;
  319.                 q_filt = q_filt_tab;
  320.                 for (m = 0; m < m_max; m++) {
  321.                     const int idx1 = i + h_SL;
  322.                     g_filt[m] = 0.0f;
  323.                     q_filt[m] = 0.0f;
  324.                     for (j = 0; j <= h_SL; j++) {
  325.                         g_filt[m] += g_temp[idx1 - j][m] * h_smooth[j];
  326.                         q_filt[m] += q_temp[idx1 - j][m] * h_smooth[j];
  327.                     }
  328.                 }
  329.             } else {
  330.                 g_filt = g_temp[i + h_SL];
  331.                 q_filt = q_temp[i];
  332.             }
  333.  
  334.             sbr->dsp.hf_g_filt(Y1[i] + kx, X_high + kx, g_filt, m_max,
  335.                                i + ENVELOPE_ADJUSTMENT_OFFSET);
  336.  
  337.             if (e != e_a[0] && e != e_a[1]) {
  338.                 sbr->dsp.hf_apply_noise[indexsine](Y1[i] + kx, sbr->s_m[e],
  339.                                                    q_filt, indexnoise,
  340.                                                    kx, m_max);
  341.             } else {
  342.                 int idx = indexsine&1;
  343.                 int A = (1-((indexsine+(kx & 1))&2));
  344.                 int B = (A^(-idx)) + idx;
  345.                 float *out = &Y1[i][kx][idx];
  346.                 float *in  = sbr->s_m[e];
  347.                 for (m = 0; m+1 < m_max; m+=2) {
  348.                     out[2*m  ] += in[m  ] * A;
  349.                     out[2*m+2] += in[m+1] * B;
  350.                 }
  351.                 if(m_max&1)
  352.                     out[2*m  ] += in[m  ] * A;
  353.             }
  354.             indexnoise = (indexnoise + m_max) & 0x1ff;
  355.             indexsine = (indexsine + 1) & 3;
  356.         }
  357.     }
  358.     ch_data->f_indexnoise = indexnoise;
  359.     ch_data->f_indexsine  = indexsine;
  360. }
  361.  
  362. #include "aacsbr_template.c"
  363.