1. 程式人生 > >webRTC AEC 非線性濾波處理演算法

webRTC AEC 非線性濾波處理演算法

先把程式碼貼上,有空時候回來註釋

static void NonLinearProcessing(AecCore* aec, short* output, short* outputH) {
  float efw[2][PART_LEN1], dfw[2][PART_LEN1], xfw[2][PART_LEN1];
  complex_t comfortNoiseHband[PART_LEN1];
  float fft[PART_LEN2];
  float scale, dtmp;
  float nlpGainHband;
  int i, j, pos;

  // Coherence and non-linear filter
  float cohde[PART_LEN1], cohxd[PART_LEN1];
  float hNlDeAvg, hNlXdAvg;
  float hNl[PART_LEN1];
  float hNlPref[kPrefBandSize];
  float hNlFb = 0, hNlFbLow = 0;
  const float prefBandQuant = 0.75f, prefBandQuantLow = 0.5f;
  const int prefBandSize = kPrefBandSize / aec->mult;
  const int minPrefBand = 4 / aec->mult;

  // Near and error power sums
  float sdSum = 0, seSum = 0;

  // Power estimate smoothing coefficients.
  const float* ptrGCoh = aec->extended_filter_enabled
                             ? kExtendedSmoothingCoefficients[aec->mult - 1]
                             : kNormalSmoothingCoefficients[aec->mult - 1];
  const float* min_overdrive = aec->extended_filter_enabled
                                   ? kExtendedMinOverDrive
                                   : kNormalMinOverDrive;

  // Filter energy
  float wfEnMax = 0, wfEn = 0;
  const int delayEstInterval = 10 * aec->mult;

  float* xfw_ptr = NULL;

  aec->delayEstCtr++;
  if (aec->delayEstCtr == delayEstInterval) {
    aec->delayEstCtr = 0;
  }

  // initialize comfort noise for H band
  memset(comfortNoiseHband, 0, sizeof(comfortNoiseHband));
  nlpGainHband = (float)0.0;
  dtmp = (float)0.0;

  // Measure energy in each filter partition to determine delay.
  // TODO: Spread by computing one partition per block?
  if (aec->delayEstCtr == 0) {
    wfEnMax = 0;
    aec->delayIdx = 0;
    for (i = 0; i < aec->num_partitions; i++) {
      pos = i * PART_LEN1;
      wfEn = 0;
      for (j = 0; j < PART_LEN1; j++) {
        wfEn += aec->wfBuf[0][pos + j] * aec->wfBuf[0][pos + j] +
                aec->wfBuf[1][pos + j] * aec->wfBuf[1][pos + j];
      }

      if (wfEn > wfEnMax) {
        wfEnMax = wfEn;
        aec->delayIdx = i;
      }
    }
  }

  // We should always have at least one element stored in |far_buf|.
  assert(WebRtc_available_read(aec->far_buf_windowed) > 0);
  // NLP
  WebRtc_ReadBuffer(aec->far_buf_windowed, (void**)&xfw_ptr, &xfw[0][0], 1);

  // TODO(bjornv): Investigate if we can reuse |far_buf_windowed| instead of
  // |xfwBuf|.
  // Buffer far.
  memcpy(aec->xfwBuf, xfw_ptr, sizeof(float) * 2 * PART_LEN1);

  // Use delayed far.
  memcpy(xfw, aec->xfwBuf + aec->delayIdx * PART_LEN1, sizeof(xfw));

  // Windowed near fft
  for (i = 0; i < PART_LEN; i++) {
    fft[i] = aec->dBuf[i] * sqrtHanning[i];
    fft[PART_LEN + i] = aec->dBuf[PART_LEN + i] * sqrtHanning[PART_LEN - i];
  }
  aec_rdft_forward_128(fft);

  dfw[1][0] = 0;
  dfw[1][PART_LEN] = 0;
  dfw[0][0] = fft[0];
  dfw[0][PART_LEN] = fft[1];
  for (i = 1; i < PART_LEN; i++) {
    dfw[0][i] = fft[2 * i];
    dfw[1][i] = fft[2 * i + 1];
  }

  // Windowed error fft
  for (i = 0; i < PART_LEN; i++) {
    fft[i] = aec->eBuf[i] * sqrtHanning[i];
    fft[PART_LEN + i] = aec->eBuf[PART_LEN + i] * sqrtHanning[PART_LEN - i];
  }
  aec_rdft_forward_128(fft);
  efw[1][0] = 0;
  efw[1][PART_LEN] = 0;
  efw[0][0] = fft[0];
  efw[0][PART_LEN] = fft[1];
  for (i = 1; i < PART_LEN; i++) {
    efw[0][i] = fft[2 * i];
    efw[1][i] = fft[2 * i + 1];
  }

  // Smoothed PSD
  for (i = 0; i < PART_LEN1; i++) {
    aec->sd[i] = ptrGCoh[0] * aec->sd[i] +
                 ptrGCoh[1] * (dfw[0][i] * dfw[0][i] + dfw[1][i] * dfw[1][i]);
    aec->se[i] = ptrGCoh[0] * aec->se[i] +
                 ptrGCoh[1] * (efw[0][i] * efw[0][i] + efw[1][i] * efw[1][i]);
    // We threshold here to protect against the ill-effects of a zero farend.
    // The threshold is not arbitrarily chosen, but balances protection and
    // adverse interaction with the algorithm's tuning.
    // TODO: investigate further why this is so sensitive.
    aec->sx[i] =
        ptrGCoh[0] * aec->sx[i] +
        ptrGCoh[1] *
            WEBRTC_SPL_MAX(xfw[0][i] * xfw[0][i] + xfw[1][i] * xfw[1][i], 15);

    aec->sde[i][0] =
        ptrGCoh[0] * aec->sde[i][0] +
        ptrGCoh[1] * (dfw[0][i] * efw[0][i] + dfw[1][i] * efw[1][i]);
    aec->sde[i][1] =
        ptrGCoh[0] * aec->sde[i][1] +
        ptrGCoh[1] * (dfw[0][i] * efw[1][i] - dfw[1][i] * efw[0][i]);

    aec->sxd[i][0] =
        ptrGCoh[0] * aec->sxd[i][0] +
        ptrGCoh[1] * (dfw[0][i] * xfw[0][i] + dfw[1][i] * xfw[1][i]);
    aec->sxd[i][1] =
        ptrGCoh[0] * aec->sxd[i][1] +
        ptrGCoh[1] * (dfw[0][i] * xfw[1][i] - dfw[1][i] * xfw[0][i]);

    sdSum += aec->sd[i];
    seSum += aec->se[i];
  }

  // Divergent filter safeguard.
  if (aec->divergeState == 0) {
    if (seSum > sdSum) {
      aec->divergeState = 1;
    }
  } else {
    if (seSum * 1.05f < sdSum) {
      aec->divergeState = 0;
    }
  }

  if (aec->divergeState == 1) {
    memcpy(efw, dfw, sizeof(efw));
  }

  if (!aec->extended_filter_enabled) {
    // Reset if error is significantly larger than nearend (13 dB).
    if (seSum > (19.95f * sdSum)) {
      memset(aec->wfBuf, 0, sizeof(aec->wfBuf));
    }
  }

  // Subband coherence
  for (i = 0; i < PART_LEN1; i++) {
    cohde[i] =
        (aec->sde[i][0] * aec->sde[i][0] + aec->sde[i][1] * aec->sde[i][1]) /
        (aec->sd[i] * aec->se[i] + 1e-10f);
    cohxd[i] =
        (aec->sxd[i][0] * aec->sxd[i][0] + aec->sxd[i][1] * aec->sxd[i][1]) /
        (aec->sx[i] * aec->sd[i] + 1e-10f);
  }

  hNlXdAvg = 0;
  for (i = minPrefBand; i < prefBandSize + minPrefBand; i++) {
    hNlXdAvg += cohxd[i];
  }
  hNlXdAvg /= prefBandSize;
  hNlXdAvg = 1 - hNlXdAvg;

  hNlDeAvg = 0;
  for (i = minPrefBand; i < prefBandSize + minPrefBand; i++) {
    hNlDeAvg += cohde[i];
  }
  hNlDeAvg /= prefBandSize;

  if (hNlXdAvg < 0.75f && hNlXdAvg < aec->hNlXdAvgMin) {
    aec->hNlXdAvgMin = hNlXdAvg;
  }

  if (hNlDeAvg > 0.98f && hNlXdAvg > 0.9f) {
    aec->stNearState = 1;
  } else if (hNlDeAvg < 0.95f || hNlXdAvg < 0.8f) {
    aec->stNearState = 0;
  }

  if (aec->hNlXdAvgMin == 1) {
    aec->echoState = 0;
    aec->overDrive = min_overdrive[aec->nlp_mode];

    if (aec->stNearState == 1) {
      memcpy(hNl, cohde, sizeof(hNl));
      hNlFb = hNlDeAvg;
      hNlFbLow = hNlDeAvg;
    } else {
      for (i = 0; i < PART_LEN1; i++) {
        hNl[i] = 1 - cohxd[i];
      }
      hNlFb = hNlXdAvg;
      hNlFbLow = hNlXdAvg;
    }
  } else {

    if (aec->stNearState == 1) {
      aec->echoState = 0;
      memcpy(hNl, cohde, sizeof(hNl));
      hNlFb = hNlDeAvg;
      hNlFbLow = hNlDeAvg;
    } else {
      aec->echoState = 1;
      for (i = 0; i < PART_LEN1; i++) {
        hNl[i] = WEBRTC_SPL_MIN(cohde[i], 1 - cohxd[i]);
      }

      // Select an order statistic from the preferred bands.
      // TODO: Using quicksort now, but a selection algorithm may be preferred.
      memcpy(hNlPref, &hNl[minPrefBand], sizeof(float) * prefBandSize);
      qsort(hNlPref, prefBandSize, sizeof(float), CmpFloat);
      hNlFb = hNlPref[(int)floor(prefBandQuant * (prefBandSize - 1))];
      hNlFbLow = hNlPref[(int)floor(prefBandQuantLow * (prefBandSize - 1))];
    }
  }

  // Track the local filter minimum to determine suppression overdrive.
  if (hNlFbLow < 0.6f && hNlFbLow < aec->hNlFbLocalMin) {
    aec->hNlFbLocalMin = hNlFbLow;
    aec->hNlFbMin = hNlFbLow;
    aec->hNlNewMin = 1;
    aec->hNlMinCtr = 0;
  }
  aec->hNlFbLocalMin =
      WEBRTC_SPL_MIN(aec->hNlFbLocalMin + 0.0008f / aec->mult, 1);
  aec->hNlXdAvgMin = WEBRTC_SPL_MIN(aec->hNlXdAvgMin + 0.0006f / aec->mult, 1);

  if (aec->hNlNewMin == 1) {
    aec->hNlMinCtr++;
  }
  if (aec->hNlMinCtr == 2) {
    aec->hNlNewMin = 0;
    aec->hNlMinCtr = 0;
    aec->overDrive =
        WEBRTC_SPL_MAX(kTargetSupp[aec->nlp_mode] /
                           ((float)log(aec->hNlFbMin + 1e-10f) + 1e-10f),
                       min_overdrive[aec->nlp_mode]);
  }

  // Smooth the overdrive.
  if (aec->overDrive < aec->overDriveSm) {
    aec->overDriveSm = 0.99f * aec->overDriveSm + 0.01f * aec->overDrive;
  } else {
    aec->overDriveSm = 0.9f * aec->overDriveSm + 0.1f * aec->overDrive;
  }

  WebRtcAec_OverdriveAndSuppress(aec, hNl, hNlFb, efw);

  // Add comfort noise.
  ComfortNoise(aec, efw, comfortNoiseHband, aec->noisePow, hNl);

  // TODO(bjornv): Investigate how to take the windowing below into account if
  // needed.
  if (aec->metricsMode == 1) {
    // Note that we have a scaling by two in the time domain |eBuf|.
    // In addition the time domain signal is windowed before transformation,
    // losing half the energy on the average. We take care of the first
    // scaling only in UpdateMetrics().
    UpdateLevel(&aec->nlpoutlevel, efw);
  }
  // Inverse error fft.
  fft[0] = efw[0][0];
  fft[1] = efw[0][PART_LEN];
  for (i = 1; i < PART_LEN; i++) {
    fft[2 * i] = efw[0][i];
    // Sign change required by Ooura fft.
    fft[2 * i + 1] = -efw[1][i];
  }
  aec_rdft_inverse_128(fft);

  // Overlap and add to obtain output.
  scale = 2.0f / PART_LEN2;
  for (i = 0; i < PART_LEN; i++) {
    fft[i] *= scale;  // fft scaling
    fft[i] = fft[i] * sqrtHanning[i] + aec->outBuf[i];

    // Saturation protection
    output[i] = (short)WEBRTC_SPL_SAT(
        WEBRTC_SPL_WORD16_MAX, fft[i], WEBRTC_SPL_WORD16_MIN);

    fft[PART_LEN + i] *= scale;  // fft scaling
    aec->outBuf[i] = fft[PART_LEN + i] * sqrtHanning[PART_LEN - i];
  }

  // For H band
  if (aec->sampFreq == 32000) {

    // H band gain
    // average nlp over low band: average over second half of freq spectrum
    // (4->8khz)
    GetHighbandGain(hNl, &nlpGainHband);

    // Inverse comfort_noise
    if (flagHbandCn == 1) {
      fft[0] = comfortNoiseHband[0][0];
      fft[1] = comfortNoiseHband[PART_LEN][0];
      for (i = 1; i < PART_LEN; i++) {
        fft[2 * i] = comfortNoiseHband[i][0];
        fft[2 * i + 1] = comfortNoiseHband[i][1];
      }
      aec_rdft_inverse_128(fft);
      scale = 2.0f / PART_LEN2;
    }

    // compute gain factor
    for (i = 0; i < PART_LEN; i++) {
      dtmp = (float)aec->dBufH[i];
      dtmp = (float)dtmp * nlpGainHband;  // for variable gain

      // add some comfort noise where Hband is attenuated
      if (flagHbandCn == 1) {
        fft[i] *= scale;  // fft scaling
        dtmp += cnScaleHband * fft[i];
      }

      // Saturation protection
      outputH[i] = (short)WEBRTC_SPL_SAT(
          WEBRTC_SPL_WORD16_MAX, dtmp, WEBRTC_SPL_WORD16_MIN);
    }
  }

  // Copy the current block to the old position.
  memcpy(aec->dBuf, aec->dBuf + PART_LEN, sizeof(float) * PART_LEN);
  memcpy(aec->eBuf, aec->eBuf + PART_LEN, sizeof(float) * PART_LEN);

  // Copy the current block to the old position for H band
  if (aec->sampFreq == 32000) {
    memcpy(aec->dBufH, aec->dBufH + PART_LEN, sizeof(float) * PART_LEN);
  }

  memmove(aec->xfwBuf + PART_LEN1,
          aec->xfwBuf,
          sizeof(aec->xfwBuf) - sizeof(complex_t) * PART_LEN1);
}

static void GetHighbandGain(const float* lambda, float* nlpGainHband) {
  int i;

  nlpGainHband[0] = (float)0.0;
  for (i = freqAvgIc; i < PART_LEN1 - 1; i++) {
    nlpGainHband[0] += lambda[i];
  }
  nlpGainHband[0] /= (float)(PART_LEN1 - 1 - freqAvgIc);
}

相關推薦

webRTC AEC 非線性濾波處理演算法

先把程式碼貼上,有空時候回來註釋 static void NonLinearProcessing(AecCore* aec, short* output, short* outputH) { float efw[2][PART_LEN1], dfw[2][PART_LE

webrtc語音增強處理演算法綜述

    作為實時音視訊通訊框架的webrtc,裡面有著豐富的語音處理演算法,其中主要涉及到AEC(聲學回聲抑制),NS(噪音抑制),AGC(自增益控制),VAD(語活檢測)和CNG(舒適噪聲)等。語音處理資料分為近端和遠端,近端資料是採集到的音訊資料,遠端資料為接受到並播放的

OpenCV4Android開發實錄(5):影象邊緣處理非線性濾波(中值、雙邊)

在OpenCV4Android開發實錄(4):影象去噪與線性濾波(均值、方框、高斯)文章中,我們較為詳細地介紹了OpenCV中幾種常用的線性濾波方法原理和相關API的使用,本文將在此基礎上繼續講解OpenCV中中值和雙邊兩種非線性濾波,以及在影象濾波過

字串處理演算法(六)求2個字串最長公共部分

基礎演算法連結快速通道,不斷更新中: 整型陣列處理演算法部分: 整型陣列處理演算法(一)按照正態分佈來排列整型陣列元素 整型陣列處理演算法(二)檔案中有一組整數,要求排序後輸出到另一個檔案中 整型陣列處理演算法(三)把一個數組裡的所有元素,插入到另一個數組的指定位置 整型陣列

線性濾波非線性濾波區別

數字影象處理線性濾波: 輸出影象fo(x,y)= T[ fi(x,y) ],T是線性運算元,即:輸出影象上每個畫素點的值都是由輸入影象各畫素點值加權求和的結果。 非線性濾波的運算元中包含了取絕對值、置零等非線性運算。 線性濾波器的原始資料與濾波結果是一種算術運算,即用加減乘除等運算實現,

通俗理解卡爾曼濾波及其演算法實現(帶例項解析)

1.簡介(Brief Introduction) 在學習卡爾曼濾波器之前,首先看看為什麼叫“卡爾曼”。跟其他著名的理論(例如傅立葉變換,泰勒級數等等)一樣,卡爾曼也是一個人的名字,而跟他們不同的是,他是個現代人! 卡爾曼全名Rudolf Emil Kalman,匈牙利數學家,1930年出生於

影象處理演算法工程師——必備技能總結

                                 影象處理演算法工程師   職位要

海量資料處理演算法—Bit-Map

分享一下我老師大神的人工智慧教程!零基礎,通俗易懂!http://blog.csdn.net/jiangjunshow 也歡迎大家轉載本篇文章。分享知識,造福人民,實現我們中華民族偉大復興!        

Python自然語言處理演算法基礎

本章主要介紹文字分析的演算法設計過程中會用到的一些技巧,我只把書中對我來說有意思的例子拿出來了。 一  遞迴 遞迴就是迴圈的一種,為了實現某種目的反覆呼叫自身。下面這個例子的有意思的地方不僅限於迭代,還用了yield,可以參考廖雪峰老師關於Yield的解釋https://www.i

影象資料處理演算法總結

本博文主要介紹了影象處理的一些基礎知識 一. 影象儲存原理 影象儲存原理主要分為5種 1. RGB顏色空間,使用最為廣泛的顏色空間。 2. CMY(K)顏色空間,主要用於印刷行業。 3. HSV/HSL(I)顏色空間,人類視覺,和畫家配色領域。 4. CIE-XYZ顏色空間,

線上處理演算法 | HDU 1003 Max Sum 最大子列和

Max Sum Time Limit: 2000/1000 MS (Java/Others)    Memory Limit: 65536/32768 K (Java/Others) Total Submission(s): 299181    Accepted Submission(s)

【OpenCV入門教程之九】 非線性濾波專場 中值濾波 雙邊濾波

                正如我們上一篇文章中講到的,線性濾波可以實現很多種不同的影象變換。然而非線性濾波,如中值濾波器和雙邊濾波器,有時可以達到更好的實現效果。鄰域運算元的其他一些例子還有對二值影象進行操作的形態學運算元,用於計算距離變換和尋找連通量的半全域性運算元。 先上一張截圖:一、理論與概念講解

海量資料處理演算法—Bloom Filter

                1. Bloom-Filter演算法簡介        Bloom-Filter,即布隆過濾器,1970年由Bloom中提出。它可以用於檢索一個元素是否在一個集合中。Bloom Filter(BF)是一種空間效率很高的隨機資料結構,它利用位陣列很簡潔地表示一個集合,並能判斷一個

驗證碼處理演算法(一)

       在面對那種有許多幹擾線或者干擾點的驗證碼,或者各種各樣的驗證碼的時候,往往一個閥值是無法精確的處理圖形驗證碼的,這裡,我們主要使用一個範圍縮圈進行畫素點的比對,因為之前有使用計算畫素的平均值K作為閾值,但是會導致部分物件畫素或者背景畫素丟失,故這

影象處理演算法之美顏

       和濾鏡一樣,美顏也是影象類app必不可少的功能之一,也有的app叫人像美容,主要包括美膚及美白等幾大功能。甚至有很多專門美顏的app,比如美顏相機什麼的,可見美顏功能需求量之大。很多女孩

影象處理演算法4——Sobel 邊緣檢測運算元

Sobel 運算元是一個離散微分運算元 (discrete differentiation operator)。 它結合了高斯平滑和微分求導,用來計算影象灰度函式的近似梯度。     影象邊緣,相素值會發生顯著的變化了。表示這一改變的一個方法是使用 導數 。 梯度值的大

PTA 資料結構題目(1):最大子列和問題(分而治之、線上處理演算法

題目來源: 問題描述: 問題分析: 對於一般的問題,原始解 都能通過一種 蠻力演算法,即窮舉法的思想得到。這題也不例外。 如果我們,把輸入的陣列,所有的子列都歷遍,並從中找出最大,即可得出我們的演算法。也就是版本一。 學習要點: 1、如何

最大連續子數列和(線上處理演算法

問題描述 最大連續子數列和一道很經典的演算法問題,給定一個數列,其中可能有正數也可能有負數,我們的任務是找出其中連續的一個子數列(不允許空序列),使它們的和儘可能大。我們一起用多種方式,逐步優化解決這個問題。 暴力方法 求出所有可能連續子列的和,時間複

大資料處理演算法三:分而治之/hash對映 + hash統計 + 堆/快速/歸併排序

百度面試題1、海量日誌資料,提取出某日訪問百度次數最多的那個IP。 IP 是32位的,最多有個2^32個IP。同樣可以採用對映的方法,比如模1000,把整個大檔案對映為1000個小檔案,再找出每個小文中出現頻率最大的 IP(可以採用hash_map進行頻率統計,然後再找出頻

驗證碼識別演算法腐蝕和二值化處理演算法

腐蝕演算法實現程式碼 private static int[][] erosion(int[][] source,int threshold){ int width = source[0].length; int height = source.leng