mp3enc_psychoacoustic_fp.c

audio-video-codecs.rar语音编解码器

/*//////////////////////////////////////////////////////////////////////////////
//
//                  INTEL CORPORATION PROPRIETARY INFORMATION
//     This software is supplied under the terms of a license agreement or
//     nondisclosure agreement with Intel Corporation and may not be copied
//     or disclosed except in accordance with the terms of that agreement.
//          Copyright(c) 2002-2007 Intel Corporation. All Rights Reserved.
//
*/

#include "mp3enc_own_fp.h"

/****************************************************************************/

static Ipp32f sprdngf(Ipp32f b1,
                      Ipp32f b2,
                      Ipp32s layer)
{
  Ipp32f tmpx,tmpy,tmpz;

  if (layer == 3) {
    tmpx = (b2 >= b1) ? 3*(b2-b1) : 1.5f*(b2-b1);
  } else {
    tmpx = 1.05f * (b2-b1);
  }

  tmpz = 8 * MIN((tmpx-0.5f)*(tmpx-0.5f) - 2*(tmpx-0.5f), 0);

  tmpy = 15.811389f + 7.5f*(tmpx + 0.474f)-
         17.5f*(Ipp32f)sqrt(1 + (tmpx+0.474f)*(tmpx+0.474f));

  return (tmpy < -100 ? 0 : (Ipp32f)pow(10.0,(tmpz + tmpy)/10.0));
}

/****************************************************************************/

static void mp3enc_psychoacousticLongWindow(MP3Enc *state,
                                            Ipp32f *input_data,
                                            Ipp32f* rsqr_long, Ipp32s ch)
{
  VM_ALIGN16_DECL(Ipp32f) sw[1024];
  Ipp32f *fft_line = state->fft_line_long[ch & 1];

  VM_ALIGN16_DECL(Ipp32f) c_w[512];
  VM_ALIGN16_DECL(Ipp32f) minthres[32];
  Ipp32f x2[MAX_PPT_LONG];
  Ipp32f e_b[MAX_PPT_LONG];
  Ipp32f c_b[MAX_PPT_LONG];
  Ipp32f ecb[MAX_PPT_LONG];
  Ipp32f ct[MAX_PPT_LONG];
  Ipp32f d_cb[MAX_PPT_LONG];
  Ipp32f tmp_ecb[MAX_PPT_LONG];
  Ipp32f tmp_ct[MAX_PPT_LONG];
  Ipp32f cb[MAX_PPT_LONG];
  Ipp32f ecb_h_limit[MAX_PPT_LONG];
  Ipp32f ecb_l_limit[MAX_PPT_LONG];
  Ipp32f tmp[NUM_UNPRED_LINES_LONG];
  Ipp32f tmp0[NUM_UNPRED_LINES_LONG];
  Ipp32f tmp1[NUM_UNPRED_LINES_LONG];
  Ipp32f tmp2[NUM_UNPRED_LINES_LONG];
  Ipp32f num[NUM_UNPRED_LINES_LONG];
  Ipp32f denum[NUM_UNPRED_LINES_LONG];
  Ipp32f r_pred[NUM_UNPRED_LINES_LONG];
  Ipp32f re_pred[NUM_UNPRED_LINES_LONG];
  Ipp32f im_pred[NUM_UNPRED_LINES_LONG];
  Ipp32f *r, *r_prev, *r_prev_prev;
  Ipp32f *re, *re_prev, *re_prev_prev;
  Ipp32f *im, *im_prev, *im_prev_prev;
  Ipp32f *tmp_dst[2];
  Ipp32f *nb,*nb_l;
  Ipp32s i, j, b, sb, index;
  Ipp32s current_f_r_index, prev_f_r_index, prev_prev_f_r_index;
  Ipp32s pa_num_ptt_long = state->pa_num_ptt_long;
  Ipp32f sum1, sum2;
  Ipp32s freq_ind, layer, sblimit_real;
  Ipp32s absthr_len[3] = {130, 126, 132};
  Ipp32s upsample = state->com.upsample;

  freq_ind = state->com.header.samplingFreq;
  layer = state->com.header.layer;

  if (ch < 2 || state->com.mc_channel)
    ippsMul_32f(input_data, state->pa_hann_window_long,
               sw, state->pa_iblen_long * 2);

  if (ch < 2 || state->com.mc_channel) {
    ippsFFTFwd_RToCCS_32f(sw, fft_line, state->pa_pFFTSpecLong, state->pa_pBuffer);
    if (upsample)
      ippsZero_32f(fft_line + (1024 >> upsample), 1024 - (1024 >> upsample));
  }

/* **** */
  fft_line[0] = fft_line[1] = fft_line[2] = fft_line[3] = 0;
  fft_line[1020] = fft_line[1021] = fft_line[1022] = fft_line[1023] = 0;
/* **** */

  current_f_r_index = state->pa_current_f_r_index[ch];
  prev_f_r_index = current_f_r_index - 1;
  if (prev_f_r_index < 0) prev_f_r_index = 2;
  prev_prev_f_r_index = prev_f_r_index - 1;
  if (prev_prev_f_r_index < 0) prev_prev_f_r_index = 2;

  r            = (Ipp32f*)state->pa_r[ch][current_f_r_index];
  r_prev       = (Ipp32f*)state->pa_r[ch][prev_f_r_index];
  r_prev_prev  = (Ipp32f*)state->pa_r[ch][prev_prev_f_r_index];

  re           = (Ipp32f*)state->pa_re[ch][current_f_r_index];
  re_prev      = (Ipp32f*)state->pa_re[ch][prev_f_r_index];
  re_prev_prev = (Ipp32f*)state->pa_re[ch][prev_prev_f_r_index];

  im           = (Ipp32f*)state->pa_im[ch][current_f_r_index];
  im_prev      = (Ipp32f*)state->pa_im[ch][prev_f_r_index];
  im_prev_prev = (Ipp32f*)state->pa_im[ch][prev_prev_f_r_index];

  ippsMagnitude_32fc((Ipp32fc *)fft_line, (Ipp32f*)r, 512);
  if (ch >= 2 && layer == 3)
    ippsSum_32f(r, 512, &(state->ms_pwr[ch-2]), ippAlgHintNone);

  ippsThreshold_LT_32f_I((Ipp32f*)r, 512, 1);

  tmp_dst[0] = re;
  tmp_dst[1] = im;

  /* Calculate the unpredictebility measure c(w)                */
  /* Some transformations:                                      */
  /* re((2*r_prev-r_prev_prev)*exp(-j(2*f_prev-f_prev_prev))) = */
  /* ((2*r_prev-r_prev_prev)/(r_prev*r_prev*r_prev_prev))*      */
  /* (2*im_prev_prev*re_prev*im_prev +                          */
  /* re_prev_prev*(re_prev*re_prev-im_prev*im_prev))            */
  /*                                                            */
  /* im((2*r_prev-r_prev_prev)*exp(-j(2*f_prev-f_prev_prev))) = */
  /* ((2*r_prev-r_prev_prev)/(r_prev*r_prev*r_prev_prev))*      */
  /* (2*re_prev_prev*re_prev*im_prev -                          */
  /* im_prev_prev*(re_prev*re_prev-im_prev*im_prev))            */
  /*                                                            */
  /* where re_prev_prev = prev_prev_r*cos(prev_prev_f),         */
  /*       im_prev_prev = prev_prev_r*sin(prev_prev_f),         */
  /*       re_prev = prev_r*cos(prev_prev_f),                   */
  /*       im_prev = prev_r*sin(prev_prev_f)                    */

  ippsDeinterleave_32f(fft_line, 2,  NUM_UNPRED_LINES_LONG, tmp_dst);

  /* tmp0 = (2*r_prev-r_prev_prev)/(r_prev*r_prev*r_prev_prev) */
  ippsMulC_32f(r_prev, (Ipp32f)2, r_pred,  NUM_UNPRED_LINES_LONG);
  ippsSub_32f_I(r_prev_prev, r_pred,  NUM_UNPRED_LINES_LONG);
  ippsMul_32f(r_prev, r_prev, denum,  NUM_UNPRED_LINES_LONG);
  ippsMul_32f_I(r_prev_prev, denum,  NUM_UNPRED_LINES_LONG);
  ippsDiv_32f(denum, r_pred, tmp0,  NUM_UNPRED_LINES_LONG);

  /* tmp1 = 2*re_prev*im_prev */
  ippsMulC_32f(re_prev, (Ipp32f)2, tmp1,  NUM_UNPRED_LINES_LONG);
  ippsMul_32f_I(im_prev, tmp1,  NUM_UNPRED_LINES_LONG);

  /* tmp2 = re_prev*re_prev-im_prev*im_prev */
  ippsMul_32f(re_prev, re_prev, tmp,  NUM_UNPRED_LINES_LONG);
  ippsMul_32f(im_prev, im_prev, tmp2,  NUM_UNPRED_LINES_LONG);
  ippsSub_32f_I(tmp, tmp2,  NUM_UNPRED_LINES_LONG);

  ippsMul_32f(im_prev_prev, tmp1, re_pred,  NUM_UNPRED_LINES_LONG);
  ippsMul_32f(re_prev_prev, tmp2, tmp,  NUM_UNPRED_LINES_LONG);
  ippsAdd_32f_I(tmp, re_pred,  NUM_UNPRED_LINES_LONG);
  ippsMul_32f_I(tmp0, re_pred,  NUM_UNPRED_LINES_LONG);

  ippsMul_32f(re_prev_prev, tmp1, im_pred,  NUM_UNPRED_LINES_LONG);
  ippsMul_32f(im_prev_prev, tmp2, tmp,  NUM_UNPRED_LINES_LONG);
  ippsSub_32f_I(tmp, im_pred,  NUM_UNPRED_LINES_LONG);
  ippsMul_32f_I(tmp0, im_pred,  NUM_UNPRED_LINES_LONG);

  ippsSub_32f(re, re_pred, tmp0,  NUM_UNPRED_LINES_LONG);
  ippsMul_32f_I(tmp0, tmp0,  NUM_UNPRED_LINES_LONG);

  ippsSub_32f(im, im_pred, tmp1,  NUM_UNPRED_LINES_LONG);
  ippsMul_32f_I(tmp1, tmp1,  NUM_UNPRED_LINES_LONG);

  ippsAdd_32f(tmp0, tmp1, num,  NUM_UNPRED_LINES_LONG);
  ippsSqrt_32f_I(num,  NUM_UNPRED_LINES_LONG);

  ippsAbs_32f(r_pred, denum,  NUM_UNPRED_LINES_LONG);
  ippsAdd_32f_I(r, denum,  NUM_UNPRED_LINES_LONG);

  ippsDiv_32f(denum, num, c_w,  NUM_UNPRED_LINES_LONG);

  ippsSet_32f(layer == 3 ? 0.4f : 0.3f, &c_w[ NUM_UNPRED_LINES_LONG],
               512 -  NUM_UNPRED_LINES_LONG);

  ippsSqr_32f((Ipp32f*)r, rsqr_long, 512);

  index = 0;

  for (b = 0; b < pa_num_ptt_long; b++) {
    Ipp32f *tmp_rsqr = &rsqr_long[index];
    Ipp32f *tmp_c_w = (Ipp32f*)&c_w[index];
    Ipp32s len;

    if (layer == 3) {
      len = mp3enc_ptbl_l[freq_ind][b].numlines;
    } else {
      len = mp3enc_ptbl_numlines_l12[freq_ind][b];
    }

    ippsSum_32f(tmp_rsqr, len, &e_b[b], ippAlgHintNone);
    ippsDotProd_32f(tmp_rsqr, tmp_c_w, len, &c_b[b]);
    index += len;
  }

  nb   = state->pa_nb_long[ch][state->pa_nb_curr_index[ch]];
  nb_l = state->pa_nb_long[ch][state->pa_nb_prev_index[ch]];

  for (b = 0; b < pa_num_ptt_long; b++) {
    Ipp32f *tmp_ptr = (Ipp32f*)state->pa_sprdngf_long + b * pa_num_ptt_long;

    ippsDotProd_32f(e_b, tmp_ptr, pa_num_ptt_long, &ecb[b]);
    ippsDotProd_32f(c_b, tmp_ptr, pa_num_ptt_long, &ct[b]);
  }

  if (layer == 3) {
    state->pa_next_frame_PE[ch] = 0;

    ippsMulC_32f(ecb, (Ipp32f)0.00964116, tmp_ecb, pa_num_ptt_long);
    ippsMulC_32f(ct, (Ipp32f)0.52280978, tmp_ct, pa_num_ptt_long);

    ippsAdd_32f(tmp_ct, tmp_ecb, cb, pa_num_ptt_long);

    ippsMulC_32f(ecb, (Ipp32f)0.25118864, ecb_h_limit, pa_num_ptt_long);
    ippsMulC_32f(ecb, (Ipp32f)0.0158489, ecb_l_limit, pa_num_ptt_long);
    ippsMaxEvery_32f_I(ecb_l_limit, cb, pa_num_ptt_long);
    ippsMinEvery_32f_I(ecb_h_limit, cb, pa_num_ptt_long);
    ippsMul_32f((Ipp32f*)state->pa_rnorm_long, cb,
                (Ipp32f*)nb, pa_num_ptt_long);

    /* instead of nb[b] = MAX( pow(10.0,pBlock->qsthr_long[b]/10.0)/32767.0, */
    /*                         MIN(nb[b],2.0*nb_l[b]));                      */
    /* we use only nb[b] = MIN(nb[b],2.0*nb_l[b]) yet                        */

    ippsMulC_32f_I(2, (Ipp32f*)nb_l, pa_num_ptt_long);
    ippsMinEvery_32f_I((Ipp32f*)nb_l, (Ipp32f*)nb, pa_num_ptt_long);

    for (b = 0; b < pa_num_ptt_long; b++) {
        if (nb[b] < e_b[b])
            state->pa_next_frame_PE[ch] -= (Ipp32f)(mp3enc_ptbl_l[freq_ind][b].numlines * log((nb[b]) / (e_b[b])));
    }

    if (state->com.ns_mode) {
      for ( sb = 0; sb < SBBND_L; sb++ ) {
        sum1 = mp3enc_p2sb_l[freq_ind][sb].w1 * e_b[mp3enc_p2sb_l[freq_ind][sb].bu]
          + mp3enc_p2sb_l[freq_ind][sb].w2 * e_b[mp3enc_p2sb_l[freq_ind][sb].bo];
        sum2 = mp3enc_p2sb_l[freq_ind][sb].w1 * nb[mp3enc_p2sb_l[freq_ind][sb].bu]
          + mp3enc_p2sb_l[freq_ind][sb].w2 * nb[mp3enc_p2sb_l[freq_ind][sb].bo];

        for ( b = mp3enc_p2sb_l[freq_ind][sb].bu+1; b < mp3enc_p2sb_l[freq_ind][sb].bo; b++ ) {
          sum1 += e_b[b];
          sum2 += nb[b];
        }

        if (sum1 > 0)
          state->pa_ratio_l_next[ch][sb] = sum2/sum1;
        else
          state->pa_ratio_l_next[ch][sb] = 0;
      }
    }
  } else {
    Ipp32f nmt = 5.5f;

    ippsDiv_32f(ecb, ct, d_cb, pa_num_ptt_long);
    ippsSqr_32f(d_cb, x2, pa_num_ptt_long);
    ippsMulC_32f_I(/*5.8734f*/4.3097f, x2, pa_num_ptt_long);
    ippsMulC_32f_I(/*-5.6788f*/-4.274f, d_cb, pa_num_ptt_long);
    ippsAdd_32f_I(x2, d_cb, pa_num_ptt_long);
    ippsAddC_32f_I(/*1.352f*/1.1002f, d_cb, pa_num_ptt_long);

    ippsThreshold_LT_32f_I(d_cb, pa_num_ptt_long, 0.0f);
    ippsThreshold_GT_32f_I(d_cb, pa_num_ptt_long, 1.0f);

    ippsMul_32f(d_cb, mp3enc_ptbl_TMN_l2[freq_ind], x2, pa_num_ptt_long);
    ippsSubCRev_32f_I(1.0f, d_cb, pa_num_ptt_long);
    ippsMulC_32f_I(nmt, d_cb, pa_num_ptt_long);
    ippsAdd_32f_I(x2, d_cb, pa_num_ptt_long);

    ippsMaxEvery_32f_I(mp3enc_ptbl_minval_l12[freq_ind], d_cb, pa_num_ptt_long);
    ippsMulC_32f_I(-0.230258509f, d_cb, pa_num_ptt_long);
    ippsExp_32f_I(d_cb, pa_num_ptt_long);
    ippsMul_32f_I(state->pa_rnorm_long, d_cb, pa_num_ptt_long);
    ippsMul_32f_I(ecb, d_cb, pa_num_ptt_long);

    index = 0;
    for (b = 0; b < pa_num_ptt_long; b++) {
      Ipp32s len = mp3enc_ptbl_numlines_l12[freq_ind][b];

      ippsSet_32f(d_cb[b], sw + index, len);
      index += len;
    }

    index = 0;
    for(j = 0; j < absthr_len[freq_ind]; j++) {
      Ipp32f max;
      Ipp32s len = mp3enc_absthr[freq_ind][j].numlines;

      for (i = 0; i < len; i++) {
        max = sw[index + i];
        if(max < mp3enc_absthr[freq_ind][j].absthr)
          max = mp3enc_absthr[freq_ind][j].absthr;
        sw[index + j] = max;
      }
      index += len;
    }

    ippsSet_32f(mp3enc_absthr[freq_ind][absthr_len[freq_ind]-1].absthr,
      sw + index, 512 - index);

    sblimit_real = state->com.sblimit_real;

    for(j = 0; j < sblimit_real; j++) {
      ippsMin_32f(sw + (j << 4), 16, &(minthres[j]));
//      ippsSum_32f(rsqr_long + (j << 4), 16, &(state->pa_snr[ch][j]), ippAlgHintNone);
      ippsSum_32f(rsqr_long + (j << 4), 16, &(state->snr[j]), ippAlgHintNone);
    }

//    ippsDiv_32f_I(minthres, state->pa_snr[ch], sblimit_real);
    ippsDiv_32f_I(minthres, state->snr, sblimit_real);
  }
}

/****************************************************************************/

static void mp3enc_psychoacousticShortWindow(MP3Enc *state,
                                             Ipp32f *input_data,
                                             Ipp32f *s_en,
                                             Ipp32s ch)
{
  VM_ALIGN16_DECL(Ipp32f) sw_short[3][256];
  Ipp32f *ptr_input_data;
  Ipp32f eb[MAX_PPT_SHORT];
  Ipp32f thr[MAX_PPT_SHORT];
  Ipp32f *r;
  Ipp32s win_counter;
  Ipp32s j, b, sb;
  Ipp32s pa_num_ptt_short = state->pa_num_ptt_short;
  Ipp32f ecb, nb, sum1, sum2;
  Ipp32s freq_ind;
  Ipp32s upsample = state->com.upsample;

  freq_ind = state->com.header.samplingFreq;

  ptr_input_data = (Ipp32f*)input_data;

  if (ch < 2)
  for (win_counter = 0; win_counter < 3; win_counter++) {
    ptr_input_data += 192;
    ippsMul_32f(ptr_input_data, state->pa_hann_window_short,
                (Ipp32f*)sw_short[win_counter], 2 * state->pa_iblen_short);
  }

  for (win_counter = 0; win_counter < 3; win_counter++) {
    Ipp32f *fft_line = state->fft_line_short[ch & 1][win_counter];

    if (ch < 2) {
      ippsFFTFwd_RToCCS_32f(sw_short[win_counter],
                            fft_line, state->pa_pFFTSpecShort,
                            state->pa_pBuffer);
      if (upsample)
        ippsZero_32f(fft_line + (256 >> upsample), 256 - (256 >> upsample));
    }

    r = state->pa_r_short[ch][state->pa_current_f_r_index[ch]][win_counter];

    ippsMagnitude_32fc((Ipp32fc *)fft_line, r, state->pa_iblen_short);
    ippsThreshold_LT_32f_I(r, state->pa_iblen_short, 1);
    ippsDotProd_32f(&r[32], &r[32], 96, &s_en[win_counter]);

    if (state->com.ns_mode) {
      for ( j=0, b = 0; b < pa_num_ptt_short; b++ ) {
        Ipp32s tp = mp3enc_ptbl_s[freq_ind][b].numlines;
        ippsDotProd_32f(&r[j], &r[j], tp, &eb[b]);
        j += tp;
      }

      for ( b = 0; b < pa_num_ptt_short; b++ ){
        Ipp32f *tmp_ptr = (Ipp32f*)state->pa_sprdngf_short + b * pa_num_ptt_short;
        ippsDotProd_32f(eb, tmp_ptr, pa_num_ptt_short, &ecb);

        nb = ecb * mp3enc_ptbl_s[freq_ind][b].SNR;
        thr[b] = IPP_MAX(mp3enc_ptbl_s[freq_ind][b].qthr,nb);
      }

      for ( sb = 0; sb < SBBND_S; sb++ ) {
        sum1 = mp3enc_p2sb_s[freq_ind][sb].w1 * eb[mp3enc_p2sb_s[freq_ind][sb].bu]
          + mp3enc_p2sb_s[freq_ind][sb].w2 * eb[mp3enc_p2sb_s[freq_ind][sb].bo];