sconceal.cpp
上传用户:zhongxx05
上传日期:2007-06-06
资源大小:33641k
文件大小:15k
源码类别:

Symbian

开发平台:

C/C++

  1. /* ***** BEGIN LICENSE BLOCK ***** 
  2.  * Version: RCSL 1.0/RPSL 1.0 
  3.  *  
  4.  * Portions Copyright (c) 1995-2002 RealNetworks, Inc. All Rights Reserved. 
  5.  *      
  6.  * The contents of this file, and the files included with this file, are 
  7.  * subject to the current version of the RealNetworks Public Source License 
  8.  * Version 1.0 (the "RPSL") available at 
  9.  * http://www.helixcommunity.org/content/rpsl unless you have licensed 
  10.  * the file under the RealNetworks Community Source License Version 1.0 
  11.  * (the "RCSL") available at http://www.helixcommunity.org/content/rcsl, 
  12.  * in which case the RCSL will apply. You may also obtain the license terms 
  13.  * directly from RealNetworks.  You may not use this file except in 
  14.  * compliance with the RPSL or, if you have a valid RCSL with RealNetworks 
  15.  * applicable to this file, the RCSL.  Please see the applicable RPSL or 
  16.  * RCSL for the rights, obligations and limitations governing use of the 
  17.  * contents of the file.  
  18.  *  
  19.  * This file is part of the Helix DNA Technology. RealNetworks is the 
  20.  * developer of the Original Code and owns the copyrights in the portions 
  21.  * it created. 
  22.  *  
  23.  * This file, and the files included with this file, is distributed and made 
  24.  * available on an 'AS IS' basis, WITHOUT WARRANTY OF ANY KIND, EITHER 
  25.  * EXPRESS OR IMPLIED, AND REALNETWORKS HEREBY DISCLAIMS ALL SUCH WARRANTIES, 
  26.  * INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY, FITNESS 
  27.  * FOR A PARTICULAR PURPOSE, QUIET ENJOYMENT OR NON-INFRINGEMENT. 
  28.  * 
  29.  * Technology Compatibility Kit Test Suite(s) Location: 
  30.  *    http://www.helixcommunity.org/content/tck 
  31.  * 
  32.  * Contributor(s): 
  33.  *  
  34.  * ***** END LICENSE BLOCK ***** */ 
  35. #include "sconceal.h"
  36. #include "hlxclib/string.h" // memcpy()
  37. #include "hlxclib/stdlib.h" // rand()
  38. #include "hlxclib/float.h"  // FLT_MIN and friends
  39. #include "hlxclib/math.h"   // sqrt, pow...
  40. #ifdef _VXWORKS
  41. #include "privatetrigP.h"
  42. #endif
  43. #if !defined(_WINDOWS) && !defined(_OPENWAVE)
  44. # define _copysign copysign
  45. #endif
  46. #ifndef max
  47. #define max(a,b) ((a)>=(b)?(a):(b))
  48. #endif
  49. #ifndef min
  50. #define min(a,b) ((a)<=(b)?(a):(b))
  51. #endif
  52. /*
  53.   nLinesMax denotes the maximum number of spectral lines this
  54.   object will ever be fed. No testing done!
  55.   Currently, we use dynamic memory allocation in here which raises
  56.   the question how we should fail if we run out of memory.
  57.   It might make sense to allocate memory statically.
  58.  */
  59. #ifdef _CARBON
  60. #pragma old_argmatch on
  61. #endif
  62. CConcealment::CConcealment(int nLinesMax)
  63. : m_offset(0), m_nLinesMax(nLinesMax), lastConcealed(0)
  64. {
  65.   int i ;
  66.   // allocate spectral history. Should really check if succeeded.
  67.   m_history[0].spectrum = new float[m_nLinesMax * sizeofBuffer];
  68.   // clear spectral history
  69.   zeroFloat(m_history[0].spectrum, m_nLinesMax * sizeofBuffer) ;
  70.   for (i = 0; i < sizeofBuffer ; i++)
  71.   {
  72.     m_history[i].spectrum        = m_history[0].spectrum + m_nLinesMax * i;
  73.     m_history[i].spectrumPresent = 1 ;
  74.     m_history[i].blockType       = 0 ;
  75.     m_history[i].nLinesActive    = 0 ;
  76.   }
  77. }
  78. CConcealment::~CConcealment()
  79. {
  80.   if (m_history[0].spectrum)
  81.     delete[] m_history[0].spectrum ;
  82. }
  83. /*
  84.   throw out old spectrum, making room for new. Implemented via a
  85.   circular buffer of pointers to spectra.
  86.  */
  87. void CConcealment::rotateHistory()
  88. {
  89.   int i ;
  90.   // kick time forward
  91.   if (++m_offset == sizeofBuffer) m_offset = 0 ;
  92.   for (i = 0 ; i < sizeofBuffer ; i++)
  93.   {
  94.     m_historyPtr[i] = &m_history[(i + m_offset) % sizeofBuffer] ;
  95.   }
  96. }
  97. /*
  98.   user-calleable function.
  99.   insert spectrum (MDCT / MLT coefficients) into our delay line.
  100.   Currently, blocktype is of Layer-3 flavor (2 == short block);
  101.   this should be generalized.
  102.   (Currently, spectra with block type of 2 are handed out unchanged,
  103.   but concealed afterwards because their layout is different. For G2,
  104.   this behaviour is unnecessary)
  105.   nLines is the highest active spectral line in this block
  106.   (i.e. the audio bandwidth)
  107.  */
  108. void CConcealment::insert(const float* s, int blocktype, int nLines)
  109. {
  110.   rotateHistory() ;
  111.   /*
  112.    * insert current information into the future
  113.    */
  114.   
  115.   if (nLines == -1)        // if we were passed -1
  116.     nLines = m_nLinesMax ; // it means "all"
  117.   nLinesActive(sizeofFuture)    = nLines ;
  118.   blockType(sizeofFuture)       = blocktype ;
  119.   spectrumPresent(sizeofFuture) = (s != 0) ;
  120.   // copy all lines containing energy
  121.   if (s) memcpy(spectrum(sizeofFuture), s, sizeof(float)*nLines); /* Flawfinder: ignore */
  122.   // zero out the rest
  123.   if (!s) nLines = 0 ;
  124.   zeroFloat(&(spectrum(sizeofFuture)[nLines]), m_nLinesMax - nLines);
  125. }
  126. /*
  127.   no spectrum present, schedule this spectrum to be concealed.
  128.   call this if you have a defect, or missing, spectrum.
  129.  */
  130. void CConcealment::insert()
  131. {
  132.   insert(0,0,0) ;
  133. }
  134. /*
  135.   retrieve a valid spectrum, concealed if necessary
  136.   returns index+1 of highest spectral line containing energy
  137.   (i.e. audio bandwidth) and Layer-3 type blocktype.
  138.  */
  139. int CConcealment::retrieve(float* s, int &blocktype)
  140. {
  141.   int linesActive ;
  142.   // if spectrum is present, return it as-is
  143.   if (spectrumPresent(0))
  144.   {
  145.     memcpy(s, spectrum(0), sizeof(float)*m_nLinesMax); /* Flawfinder: ignore */
  146.     blocktype = blockType(0) ;
  147.     linesActive = nLinesActive(0) ;
  148.   }
  149.   // conceal this spectrum if it was a short block or damaged
  150.   if (!spectrumPresent(0) ||
  151.        blockType(0) == 2)
  152.   {
  153.     conceal(CONCEAL_PREDICTION);
  154.     lastConcealed = 1 ;
  155.   }
  156.   else
  157.   {
  158.     lastConcealed = 0 ;
  159.   }
  160.   // if spectrum was damaged, return concealed version
  161.   if (!spectrumPresent(0))
  162.   {
  163.     memcpy(s, spectrum(0), sizeof(float)*m_nLinesMax); /* Flawfinder: ignore */
  164.     blocktype = blockType(0) ;
  165.     linesActive = nLinesActive(0) ;
  166.   }
  167.   // if this spectrum was concealed because it was a short spectrum,
  168.   // mark it as non-present (so it won't be included in the prediction)
  169.   if (blockType(0) == 2)
  170.   {
  171.     spectrumPresent(0) = 0 ;
  172.     blockType(0) = 0 ;
  173.   }
  174.   // this would be a good time to update predictor coefficients
  175.   // (if we knew how to)
  176.   // return number of lines containing energy
  177.   return linesActive ;
  178. }
  179. /*
  180.   a dispatcher between different methods of concealment. This can
  181.   probably go away once I've settled on one method.
  182.  */
  183. void CConcealment::conceal(int method)
  184. {
  185.   switch(method)
  186.   {
  187.   case CONCEAL_PREDICTION:
  188.     concealByPrediction();
  189.     break ;
  190.   }
  191. }
  192. /*
  193.   magic numbers for concealment-by-prediction
  194.  */
  195. #define alpha 1.0f // weighting factor old vs. newer spectra
  196. #define kk    1.0f // max allowed deviation of prediction from precious values
  197. void CConcealment::concealByPrediction()
  198. {
  199.   prediction() ;
  200.   /* adapt energy */
  201.   adaptEnergy() ;
  202. }
  203. void CConcealment::prediction()
  204. {
  205.   int i ;
  206.   int line ;
  207.   int linesActive ;
  208.   int finalLinesActive = 0;
  209.   float weight = (float)pow(1.0f/alpha,sizeofHistory-2) ;
  210.   if (!lastConcealed)
  211.   {
  212.     zeroFloat(g0,  CONCEALLIMIT) ;
  213.     zeroFloat(g0p, CONCEALLIMIT) ;
  214.     zeroFloat(g1,  CONCEALLIMIT) ;
  215.     zeroFloat(g1p, CONCEALLIMIT) ;
  216.     zeroFloat(g2,  CONCEALLIMIT) ;
  217.     zeroFloat(ma,  CONCEALLIMIT) ;
  218.     const float* spec0, *spec1, *spec2 ;
  219.     spec1 = spectrum(-sizeofHistory) ;
  220.     spec2 = spectrum(-sizeofHistory+1) ;
  221.     linesActive = max(nLinesActive(-sizeofHistory),nLinesActive(-sizeofHistory+1));
  222.     for (line = 0; line < linesActive ; line++)
  223.     {
  224.       ma[line] = (fabs(spec1[line]) >= fabs(spec2[line])) ?
  225.            (float)fabs(spec1[line]) : (float)fabs(spec2[line]) ;
  226.     }
  227.     for (i = -sizeofHistory+2; i < 0; i++)
  228.     {
  229.       spec0 = spec1;
  230.       spec1 = spec2;
  231.       spec2 = spectrum(i) ;
  232.       if (spectrumPresent(i))
  233.       {
  234.         linesActive = max(nLinesActive(i),nLinesActive(i-1));
  235.         linesActive = max(linesActive,nLinesActive(i-2));
  236.         if (linesActive > finalLinesActive)
  237.           finalLinesActive = linesActive ;
  238.         for (line = 0; line < linesActive ; line++)
  239.         {
  240.           float s0 = spec0[line] ;
  241.           float s1 = spec1[line] ;
  242.           float s2 = spec2[line] ;
  243.           if ((float)fabs(s2) > ma[line]) ma[line] = (float)fabs(s2) ;
  244.           g0[line]  += s1 * s1 * weight;
  245.           g0p[line] += s0 * s0 * weight;
  246.           g1[line]  += s0 * s1 * weight;
  247.           g1p[line] += s1 * s2 * weight;
  248.           g2[line]  += s0 * s2 * weight;
  249.         } // for (line)
  250.       } // if spectrumPresent
  251.       weight *= alpha ;
  252.     } // for (i over history)
  253.     // now calculate the prediction coefficients
  254.     calculateCoefficients(finalLinesActive,g0,g0p,g1,g1p,g2,a0,a1);
  255.   }
  256.   else
  257.   {
  258.     finalLinesActive = nLinesActive(-1) ;
  259.   }
  260.   predict(finalLinesActive, spectrum(-2),spectrum(-1),a0,a1,ma,spectrum(0)) ;
  261.   if (!lastConcealed)
  262.   {
  263.     // take the future into account!
  264.     weight = 1.0f; // make it important!
  265.     if (spectrumPresent(1))
  266.     {
  267.       const float* spec0, *spec1, *spec2 ;
  268.       spec0 = spectrum(-1);
  269.       spec1 = spectrum(0);
  270.       spec2 = spectrum(1);
  271.       linesActive = max(nLinesActive(1),nLinesActive(0));
  272.       linesActive = max(linesActive,nLinesActive(-1));
  273.       if (linesActive > finalLinesActive)
  274.         finalLinesActive = linesActive ;
  275.       for (line = 0; line < linesActive ; line++)
  276.       {
  277.         float s0 = spec0[line] ;
  278.         float s1 = spec1[line] ;
  279.         float s2 = spec2[line] ;
  280.         g0[line]  += s1 * s1 * weight;
  281.         g0p[line] += s0 * s0 * weight;
  282.         g1[line]  += s0 * s1 * weight;
  283.         g1p[line] += s1 * s2 * weight;
  284.         g2[line]  += s0 * s2 * weight;
  285.         if (fabs(s2) > ma[line]) ma[line] = (float)fabs(s2) ;
  286.       } // for (line)
  287.     } // if spectrumPresent
  288.     // calculate the prediction coefficients again
  289.     calculateCoefficients(finalLinesActive,g0,g0p,g1,g1p,g2,a0,a1);
  290.     // the actual prediction, redone
  291.     predict(finalLinesActive,spectrum(-2),spectrum(-1),a0,a1,ma,spectrum(0)) ;
  292.   }
  293.   // clear spectrum above highest predicted line
  294.   zeroFloat(spectrum(0) + finalLinesActive, m_nLinesMax - finalLinesActive) ;
  295.   nLinesActive(0) = finalLinesActive ;
  296. }
  297. /*
  298.   tile the spectrum into equally-sized bands, adapting the energy
  299.   within to be equal to either
  300.   - the energy in the last frame
  301.   - geometric mean of last and future frame energies if available
  302.  */
  303. void CConcealment::adaptEnergy()
  304. {
  305.   enum
  306.   {
  307.     nBands = 100
  308.   };
  309.   int bWidth = m_nLinesMax / nBands;
  310.   int band ;
  311.   for (band = 0 ; band < nBands ; band++)
  312.   {
  313.     float nrgm1 ;
  314.     float nrg0  ;
  315.     float nrgp1 ;
  316.     int j, jstart = band * bWidth, jstop = min(m_nLinesMax, jstart + bWidth) ;
  317.     float nrg;
  318.     /* energy per band of last spectrum */
  319.     nrg = 0.0f ;
  320.     for (j = jstart ; j < jstop ; j++)
  321.     {
  322.       nrg += spectrum(-1)[j]*spectrum(-1)[j] ;
  323.     }
  324.     nrgm1 = nrg ;
  325.     /* energy per band of current (concealed/predicted) spectrum */
  326.     nrg = 0.0f ;
  327.     for (j = jstart ; j < jstop ; j++)
  328.     {
  329.       nrg += spectrum(0)[j]*spectrum(0)[j] ;
  330.     }
  331.     nrg0 = nrg ;
  332.     /* energy per band of next spectrum, if present */
  333.     int future = spectrumPresent(1) ? 1 : (spectrumPresent(2) ? 2 : 0) ;
  334.     if (future)
  335.     {
  336.       nrg = 0.0f ;
  337.       for (j = jstart ; j < jstop ; j++)
  338.       {
  339.         nrg += spectrum(future)[j]*spectrum(future)[j] ;
  340.       }
  341.       nrgp1 = nrg ;
  342.     }
  343.     else
  344.     {
  345.       nrgp1 = nrgm1 ;
  346.     }
  347.     /* adapt spectral energy. If we have two spectra, we adapt to the geometric
  348.        mean of signal energies. If there is just one (the last), we slowly fade
  349.        out
  350.     */
  351.     if (nrg0 > FLT_MIN)
  352.       /* if energy is zero already, don't adjust energy */
  353.     {
  354.       double templ  = pow(nrgm1, ((float)future)/(future+1)) * pow(nrgp1, ((float)1.0)/(future+1)) ;
  355.       double factor = sqrt(templ / nrg0) ;
  356.       if (factor > 10.0) factor = 10.0 ;
  357.       for (j = jstart ; j < jstop ; j++)
  358.       {
  359.         spectrum(0)[j] *= (float)factor ;
  360.       }
  361.     }
  362.   }
  363. }
  364. /*
  365.   Given prediction coefficients, do the actual prediction.
  366.  */
  367. void CConcealment::predict(/* input */
  368.                            int          nLines,
  369.                            const float *spec0,
  370.                            const float *spec1,
  371.                            const float *a0,
  372.                            const float *a1,
  373.                            const float *ma,
  374.                            /* output */
  375.                            float       *result)
  376. {
  377.   int line ;
  378.   for (line = 0 ; line < nLines ; line++)
  379.   {
  380.     float predict = spec0[line]*a1[line] + spec1[line]*a0[line] ;
  381.     /* sanity check on prediction: Do not let predicted value stray from interval
  382.       spanned by past mdct values */
  383.     if (fabs(predict) > kk*ma[line]) predict = (float)_copysign(kk*ma[line], predict) ;
  384.     result[line] = predict ;
  385.   }
  386. }
  387. /*
  388.   derive prediction coefficients from intermediate calculated values
  389.  */
  390. void CConcealment::calculateCoefficients(/* input */
  391.                                          int          nLines,
  392.                                          const float *g0,
  393.                                          const float *g0p,
  394.                                          const float *g1,
  395.                                          const float *g1p,
  396.                                          const float *g2,
  397.                                          /* output */
  398.                                          float       *a0,
  399.                                          float       *a1)
  400. {
  401.   int line ;
  402.   for (line = 0; line < nLines ; line++)
  403.   {
  404.     float den = (g0[line] *g0p[line] - g1[line]*g1[line] ) ;
  405.     float _a0 = (g1p[line]*g0p[line] - g2[line]*g1[line] ) ;
  406.     float _a1 = (g2[line] *g0[line]  - g1p[line]*g1[line]) ;
  407.     if (fabs(den) > 1.0f / FLT_MAX)
  408.     {
  409.       den = 1.0f / den ;
  410.       _a0 *= den ;
  411.       _a1 *= den ;
  412.     }
  413.     else
  414.     {
  415.       _a0 = 1.0f; _a1 = 0.0f; // repetition if no estimate possible
  416.     }
  417.     // limit the magnitude of the pole(s) to less than 2
  418.     poleLimitMagnitude(2.0f,&_a0,&_a1);
  419.     a0[line] = _a0 ; a1[line] = _a1 ; // update predictors
  420.   }
  421. }
  422. /*
  423.   limit the magnitudes of the two poles of the error concealment
  424.   prediction filter to make it stable.
  425.   Only works for two-pole filters.
  426.  */
  427. void CConcealment::poleLimitMagnitude(float mag, float *a0, float *a1)
  428. {
  429.   if (fabs(*a1) > mag)
  430.   {
  431.     *a0 /= (float)sqrt(fabs(*a1/mag));
  432.     *a1  = -mag ;
  433.   }
  434. }
  435. /*
  436.   clear (set to zero) a float vector
  437.  */
  438. void CConcealment::zeroFloat(float *f, int n)
  439. {
  440.   int i ;
  441.   for (i = 0 ; i < n ; i++) f[i] = 0.0f ;
  442. }
  443. #ifdef _CARBON
  444. #pragma old_argmatch reset
  445. #endif