天天看點

webRTC AEC 濾波器的自适應

權值更新就是使用SSE指令集加速完成NLMS算法的頻域計算,具體實作細節見代碼注釋。

static void FilterAdaptationSSE2(AecCore* aec,
                                 float* fft,
                                 float ef[2][PART_LEN1]) {
  int i, j;
  const int num_partitions = aec->num_partitions;  // 分塊數
  for (i = 0; i < num_partitions; i++) {
	  // xpos指向分塊的第一個元素,每塊65個元素
    int xPos = (i + aec->xfBufBlockPos) * (PART_LEN1);
    int pos = i * PART_LEN1;
    // 邊界檢查,如果分塊位置超過塊數,回到緩存的第一塊的起始位址
    if (i + aec->xfBufBlockPos >= num_partitions) {
      xPos -= num_partitions * PART_LEN1;
    }

    // Process the whole array...
    for (j = 0; j < PART_LEN; j += 4) {
      // Load xfBuf(遠場資料的頻域表達) and ef(消回聲估計誤差的頻域表達).

		/*_mm_load_ps它的輸入是一個指向float的指針,傳回的就是一個__m128類型的資料,從函數的角度了解,
		就是把一個float數組的四個元素依次讀取,傳回一個組合的__m128類型的SSE資料類型,
		進而可以使用這個傳回的結果傳遞給其它的SSE指令進行運算,比如加法等;從彙編的角度了解,
		它對應的就是讀取記憶體中連續四個位址的float資料,将其放入SSE新的暫存器(XMM0~8)中,
		進而給其他的指令準備好資料進行計算,_ps表示并行。*/

		// 加載資料
      const __m128 xfBuf_re = _mm_loadu_ps(&aec->xfBuf[0][xPos + j]);
      const __m128 xfBuf_im = _mm_loadu_ps(&aec->xfBuf[1][xPos + j]);
      const __m128 ef_re = _mm_loadu_ps(&ef[0][j]);
      const __m128 ef_im = _mm_loadu_ps(&ef[1][j]);

      // Calculate the product of conjugate(xfBuf) by ef.
      //   re(conjugate(a) * b) = aRe * bRe + aIm * bIm 實部
      //   im(conjugate(a) * b)=  aRe * bIm - aIm * bRe 虛部
	  /*__m128 _mm_mul_ps(__m128 a , __m128 b )将a和b的四個單精度浮點值相乘。
	  傳回值r0 := a0 * b0;r1 := a1 * b1;r2 := a2 * b2;r3 := a3 * b3*/
      const __m128 a = _mm_mul_ps(xfBuf_re, ef_re);	//	遠場資料實部與誤差實部相乘
      const __m128 b = _mm_mul_ps(xfBuf_im, ef_im);	//  遠場資料虛部與誤差虛部相乘
      const __m128 c = _mm_mul_ps(xfBuf_re, ef_im); //  遠場資料實部與誤差虛部相乘
      const __m128 d = _mm_mul_ps(xfBuf_im, ef_re); //  遠場資料虛部與誤差實部相乘
      const __m128 e = _mm_add_ps(a, b);			//  a+b
      const __m128 f = _mm_sub_ps(c, d);			//  c-d

      // Interleave real and imaginary parts.
		// extern __m128 _mm_unpacklo_ps(__m128 _A, __m128 _B);
	    // 傳回一個__m128的寄存器,Selects and interleaves the lower two single-precision,  
		// floating-point values from a and b  
		// r0=_A0, r1=_B0, r2=_A1, r3=_B1  
      const __m128 g = _mm_unpacklo_ps(e, f);

	  //	extern __m128 _mm_unpackhi_ps(__m128 _A, __m128 _B);
	  //	傳回一個__m128的寄存器,Selects and interleaves the upper two single-precision,  
	  //	floating-point values from a and b  
	  //	r0=_A2, r1=_B2, r2=_A3, r3=_B3 
      const __m128 h = _mm_unpackhi_ps(e, f);

      // Store,extern void _mm_storeu_ps(float *_V, __m128 _A); 
	  //傳回為空,Stores four single-precision, floating-point values,  
	  //The address does not need to be 16-byte aligned  
	  //_V[0]=_A0, _V[1]=_A1, _V[2]=_A2, _V[3]=_A3  
      _mm_storeu_ps(&fft[2 * j + 0], g);
      _mm_storeu_ps(&fft[2 * j + 4], h);
    }

    // ... and fixup the first imaginary entry.
    fft[1] = MulRe(aec->xfBuf[0][xPos + PART_LEN],
                   -aec->xfBuf[1][xPos + PART_LEN],
                   ef[0][PART_LEN],
                   ef[1][PART_LEN]);

	// 對fft進行傅裡葉逆變換
    aec_rdft_inverse_128(fft);
    memset(fft + PART_LEN, 0, sizeof(float) * PART_LEN);

    // fft scaling
    {
      float scale = 2.0f / PART_LEN2;
      const __m128 scale_ps = _mm_load_ps1(&scale);
      for (j = 0; j < PART_LEN; j += 4) {
        const __m128 fft_ps = _mm_loadu_ps(&fft[j]);
        const __m128 fft_scale = _mm_mul_ps(fft_ps, scale_ps);
        _mm_storeu_ps(&fft[j], fft_scale);
      }
    }
    aec_rdft_forward_128(fft);

	// 權值更新,權值為頻域表達形式
    {
      float wt1 = aec->wfBuf[1][pos];
      aec->wfBuf[0][pos + PART_LEN] += fft[1];
      for (j = 0; j < PART_LEN; j += 4) {
        __m128 wtBuf_re = _mm_loadu_ps(&aec->wfBuf[0][pos + j]);
        __m128 wtBuf_im = _mm_loadu_ps(&aec->wfBuf[1][pos + j]);
        const __m128 fft0 = _mm_loadu_ps(&fft[2 * j + 0]);
        const __m128 fft4 = _mm_loadu_ps(&fft[2 * j + 4]);
        const __m128 fft_re =
            _mm_shuffle_ps(fft0, fft4, _MM_SHUFFLE(2, 0, 2, 0));
        const __m128 fft_im =
            _mm_shuffle_ps(fft0, fft4, _MM_SHUFFLE(3, 1, 3, 1));
        wtBuf_re = _mm_add_ps(wtBuf_re, fft_re);
        wtBuf_im = _mm_add_ps(wtBuf_im, fft_im);
        _mm_storeu_ps(&aec->wfBuf[0][pos + j], wtBuf_re);
        _mm_storeu_ps(&aec->wfBuf[1][pos + j], wtBuf_im);
      }
      aec->wfBuf[1][pos] = wt1;
    }
  }
}
           

繼續閱讀