语音处理：音频信号采样点白化方法初探

附：C代码在线调试工具

公式

滑动平均滤波器(Moving Average Filter, MA)：


公式说明如下：

  • Xi表示信号幅值，带尖的Xi表示白化后的值
  • Envi指一个滑窗内能量大小和
  • a通常取3
  • log转换是为了降低计算复杂度，变sqrt为log和幂运算

效果说明

  • 效果类似滑动平均法，整体会将能量缩小些
  • 将使能量曲线变平坦，拉平异常信号尖峰
  • 注意事项，使用该公式Xi绝对值必须大于1才有效。论证见效果实现中的右图。

大家可以思考下，该如何改进，可以让该公式或者方法对(-1, 1)范围内的值也能起到白化效果？个人思考见效果改进一节，供参考。

功能实现C代码

MA滑窗的C代码实现，还可以优化，去掉里面一个for循环，有兴趣的读者可以试试。利用以下代码实现对应公式，将处理后的数据导入Python进行折线绘制，观察效果。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

void SignalWhitening(float *mdct, int start, int stop) 
{
   
    int i, j, n;
    float div, envi;
    int a = 3; // 滑窗左边宽3， [i-3, i+3]
    div = 0;
    for (i = start; i <= stop - a; i++) {
   
        envi = 1e-3f;
        for (j = i - a; j < i + a + 1; j++) {
    // 滑窗
            envi += mdct[j] * mdct[j];
        }
        envi /= (float)(a * 2 + 1);

        // n = max(0.f, (int)(log(envi) * INV_LOG_2)); /* INV_LOG_2 = 1 / (float)log(2.0f)) */
        // n >>= 1; /* sqrt() */
        // div = (float)(pow(2.0f, (float)-n)); // -0.5log2(Envi)
        envi = sqrt(envi);
        div = 1 / envi;

        mdct[i] *= div;                                       /* same as shift */
    }
    for (; i <= stop; i++) {
   
        envi = 1e-3f;
        for (j = i - a; j < stop; j++) {
   
            envi += mdct[j] * mdct[j];
        }
        envi /= (float)(stop - (i - a));

        // n = max(0.f,(int)(log(envi) * INV_LOG_2)); /* INV_LOG_2 = 1 / (float)log(2.0f)) */
        // n >>= 1; /* sqrt() */
        // div = (float)(pow(2.0f, (float)(-n)));
        envi = sqrt(envi);
        div = 1 / envi;

        mdct[i] *= div;                                       /* same as shift */
    }
    return; 
}

int main(void)
{
   
    float f1[10] = {
   1.2, 3.7, 4, 1, 3, 5.7, 10, 4.5, 5.3, 3.1};
    float f2[10] = {
   0.2, 0.7, 0.3, 0.1, 0.5, 0.0008, 0.7, 0.7, 0.4, 0.6};
    SignalWhitening(f1, 0, 10);
    int i;
    for (i = 0; i < 10; i++) {
   
        printf("%f, ", f1[i]);
    }
    printf("\n");
    SignalWhitening(f2, 0, 10);
    for (i = 0; i < 10; i++) {
   
        printf("%f, ", f2[i]);
    }

    return 0;
}

效果实现Python代码

容易从下图方框中看到，异常的尖峰被平滑了，能量调整得更均匀，曲线平整许多，但曲线形状和趋势不会改变。效果如下：

注：上面右图标注无误，黄线确为 f2 after，蓝线为 f2 before。可以看出，当Xi的值都在(-1,1)范围时，会起反效果，变得更尖，故要注意应用的场景。

实现代码如下：

import numpy as np
import matplotlib.pyplot as plt

if __name__ == '__main__':
## 处理前的数据 
    x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] ## x坐标表示各时刻
    f1 = [1.2, 3.7, 4, 1, 3, 5.7, 10, 4.5, 5.3, 3.1] ## 信号幅值
    f2 = [0.2, 0.7, 0.3, 0.1, 0.5, 0.0008, 0.7, 0.7, 0.4, 0.6]
## 白化处理后的数据
    f1_1 = [0.556992,1.547797,1.352833,0.218198,0.616024,1.114535,2.094136,1.488723,1.756067,1.411408]
    f2_1 = [0.666138,1.633952,0.426430,0.131527,0.653326,0.001592,1.327920,1.035118,0.488252,0.647540]

## legend设置
    l1 = plt.plot(x,f1,'r-',label='f1 before')
    l2 = plt.plot(x,f1_1,'g-',label='f1 after')
## 画折线
    plt.plot(x, f1,'ro-', x, f1_1,'g+-')
    plt.legend()
    
## 新起一张图表
    plt.show()  
    l3 = plt.plot(x,f2,'b-',label='f2 before')
    l4 = plt.plot(x,f2_1,'y-',label='f2 after')
    plt.plot(x, f2,'b^-', x, f2_1, 'y+-')
    plt.legend()

效果改进

当Xi在(-1,1)范围时，仅需要将Xi / sqrt(Envi)变为Xi * sqrt(Envi)即可，因为第一个结果被归一化，绝对值在1范围附近，第二个已经在1范围内了，只要调节下系数即可。

代码修改，仅需将C代码中22行和37行代码

        div =  1 / envi;

改为

        if (abs(envi) < 1) {
   
            div =  envi;
        } else {
   
            div =  1 / envi;
        }

即可实现对任意情况的Xi调节。针对(-1, 1)范围时的值，调整效果如下：

原来有问题的右图，变成如今正确的右图。

遗留问题

  • 公式里envi开根号前应该是要除以2a + 1的，代码里已调整
  • Xi / sqrt(envi)，sqrt(envi)本质就是平均化后的一个Xi，感觉结果被归一化了，单位量纲为1，要注意

参考资料