權值更新就是使用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;
}
}
}