RNNoise超详细解析

rnnosise笔记

论文阅读

​ 原始论文为 Jean-Marc Valin 等发表的A Hybrid DSP/Deep Learning Approach to Real-Time Full-Band Speech Enhancement一文，该文章提出一种将传统信号处理与深度学习结合的语音增强算法。项目地址为：GitHub - xiph/rnnoise: Recurrent neural network for audio noise reduction，下面简单总结一下论文。

1. INTRODUCTION

   ​ 略过

2. SIGNAL MODEL

   ​ 作者提出了一种用于噪声抑制的快速算法，在降噪时需要精细调整的部分使用深度学习方法，而其他部分使用传统的DSP方法。使用的窗长为20ms，窗之间的overlap为10ms，语音帧的分析（analysis）与合成（synthesis）用到的窗函数为Vorbis Window，定义为$\omega (n)=\sin [\frac{\pi }{2}{\sin }^2(\frac{\pi n}{N})]$，其中$N$为窗长。噪声抑制的主要部分是将RNN计算出的增益作用于分辨率较低的噪声频谱包络。后面还用pitch filter进行进一步地优化。

   A. band struction

   ​ 其他论文中用神经网络直接估计frequency bins需要的网络复杂度较高，从而计算量较大。为了避免此问题，作者假定频谱包络足够平坦，进而可以使用比较粗糙的分辨率。此外，并没有直接计算频谱幅度，而是对理想临界带增益（ideal critical band gains）进行估计。频带划分选择和Opus codec使用的Bark scale相同（实际上为了方便，文章作者直接使用了Opus的pitch计算代码），在低频区，每个频带最少有4个bins，并且使用的是三角频带（滤波）而非矩形频带，每个三角的峰值和其相邻三角的边界点重合。最终band的数量为22，也即网络输出为22个$[0,1]$的值。

   ​ 用${\omega }_b(k)$表示第$b$个band在频率$k$出的幅度，有${\sum }_b{\omega }_b(k)=1$，对于频域信号$X(k)$，某一个band的能量为$E(b)={\sum }_k{\omega }_b(k)|X(k){|}^2$，每个band的增益定义为$g_b=\sqrt{\frac{E_s(b)}{E_x(b)}}$，其中$E_s(b)$的纯净语音的band能量，$E_x(b)$的带噪语音的band能量。降噪即是对理想增益${\hat{g}}_b$进行插值$r(k)=\sum {\omega }_b(k){\hat{g}}_b$，然后将每个$r(k)$作用于第$k$个频点上。

   B. Pitch filtering

   ​ 上面提到的频带分辨率较低，无法获得更多更加细节的信息，这会对噪声抑制（noise suppression between pitch harmonics）产生不利影响，我们可以利用一个梳状滤波器来抑制谐波之间的噪声。每个频带的滤波器系数己为${\alpha }_b$，用$P(k)$表示基音延迟的信号$x(n-pindex)$进行DFT变换后的信号（原始信号为$x(n)$，进行pitch tracking后得到的index为pindex，则相应的基音延时信号位哦$x(n-pindex)$）。滤波操作为计算$X(k)+{\alpha }_{b}P(k)$。

   ​ 基音相关度为$p_b=\frac{{\sum }_k{\omega }_b(k)\Re [X(k)P^{\ast }(k)]}{\sqrt{{\sum }_k{\omega }_b(k)|X(k){|}^2\cdot {\sum }_k{\omega }_b(k)|P(k){|}^2}}$（可以理解为$x(n)$和$x(n-pindex)$在频带$b$的相似度），此公式等效于在时域上求相关。公式中的花体R表示取实部。

   ​ 想要获得最优的滤波器系数是比较困难的，使用最小化均方误差的算法也未必有较好的效果，作者采用了一种启发式的算法。因为噪声会导致基音相关度（pitch correlation）的降低，且平均来讲$p_b>g_b$，因而考虑以下几种情况：对于$p_b>g_b$的频带，说明$x(n)$本身噪声较大，需要利用pitch延时信号来加强语音，因而取${\alpha }_b=1$；当没有噪声时，尽量不使信号发生畸变，因而取${\alpha }_b=0$；当$p_b=0$也即基音相关度为0时，也取${\alpha }_b=0$。可以使用一个式子来包含这些情况：${\alpha }_b=\min (\sqrt{\frac{p_b^2(1-g_b^2)}{(1-p_b^2)g_b^2}},1)$。上面提到的滤波器为FIR滤波器，也可以使用IIR滤波器，这里不做介绍。

   C. Feature extraction

   ​ 前22个特征由22个频带能量做对数变换再做DCT变换得到，也即22个BFCC（Bark-frequency cepstral coefficients）。此外取BFCC的前六个一阶差分和二阶差分作为第23至34个特征，再将前六个基因相关度作为第35至第40个特征。最后两个特征为基音周期（pitch tracking得到的结果）和基音平稳度。因此输入特征一共有42个。

3. DEEP LEARNING ARCHITECTURE

   ​ 神经网络的数据如下图所示，主要包括三个循环层，作用分别为VAD、噪声谱估计、噪声消除，实验表明，实验表明使用GRU网络的效果优于LSTM网络。

​ 论文剩下部分略过

整体计算流程分析

​ 上图为rnnoise的整体计算流程。需要解释叙述的包括Pitch analysis、Pitch filtering、Feature extraction、Band gain interpolation几个部分。Pitch analysis即为基音（音高）追踪（pitch tracking），追踪得到的基音周期一方面用于Feature extraction作为一个特征，另一方面利用其得到$x(n-pindex)$进而进行后面的基音滤波（pitch filtering）操作。Feature extraction即为论文中求取42维特征的部分，将其作为输入送入RNN，输出为22维数据，其target便是前文提到的$g_b$。每一帧的增益为22维，因而为了将其作用于点数为481（帧长＋1）的$X(k)$上，需要对其进行插值，即Band gain interpolation。作用于$X$后再经过包括IFFT等一系列变换得到增强后的音频。下面对一些步骤进行详细分析。

Pitch analysis

​ 上图为一段语音（带噪）的频谱图，如果我们分析对其分析会发现人声有很明显的共振特征（图中高亮条状即为共振峰），这是由于声带振动所产生的。在很短的时间内，声带振动的频率是比较平稳的，也即有一个稳定的基频，pitch analysis的目的就是寻找该基频。因为基频的存在，语音的时域信号在很短的时间窗口内可以认为是有周期性的。如下图所示，对于一条语音，将红线部分放大便得到下面的时域图，可以看出信号在短时间内有很强的周期性。我们将某一个峰值对应的index记为index 0，将其下一个峰值记为index 1，音高追踪（pitch analysis/ pitch tracking）的目的便是得到两个index的差值。

周期检测函数（PDF, periodicity detection function）

​ 事实上，上述index差值的求解常常用到周期检测函数，常用的周期检测函数有：自相关函数、归一化方差函数、平均幅度方差函数等，如下：

1. 自相关函数（ACF, Autocorrelation function）：假设原本的frame为$s(t),t=0...n-1$，移动后的frame为$s(t-\tau )$也即向后移动$\tau$个采样点，则$acf(\tau )=\sum _{t=\tau }^{n-1}s(t)s(t-\tau )$ , 即对两帧的重合部分相乘相加。重合的长度为$n-\tau$，显然$\tau$越大，重合长度越短，则其实质为截断（truncated）ACF. 当然ACF还有多种不同形式，如归一化（即除以重叠长度）、半帧移等。

2. 归一化方差（Normalized squared difference）

   $nsdf(\tau )=\frac{2\sum s(t)s(t-\tau )}{\sum s^2(t)+\sum s^2(t-\tau )}$, 可以验证 $-1\le nsdf(\tau )\le 1$

3. 平均幅度方差函数（AMDF，Average Magnitude Difference Function）：$amdf(\tau )=\sum _{t=\tau }^{n-1}|s(t)-s(t-\tau )|$

rnnoise实际使用的是类似归一化方差的周期检测函数。

SIFT（simple inverse filter tracking）

​ 通道效应（channel effect）可能会导致pitch track出现错误，即声带信号经过口腔鼻腔等会发生变形，希望找到声带原始信号。我们可将当前信号$s(n)$表示为前面$m$个信号的线性组合: $s(n)=a_{1}s(n-1)+a_{2}s(n-2)+...+s_{m}s(n-m)+e(n)$，利用平方最小法来寻找最佳的系数$a_1,a_2,a_3,...,a_m$，使得$\sum e^2(n)$最小，则$e(n)$便是原始激发信号（excitation signal）。即：线性预测误差信号就是原始激发信号，我们利用原始激发信号，便可以更准确地进行音高追踪。

整体流程

​ 上图为音高追踪的整体流程，首先进行低通滤波和降采样操作，接着对数据求自相关，利用平方最小法进行线性预测分析（LPC analysis）得到线性预测系数（LPC cofficient），将之作用于低通滤波后的信号便得到了LPC残差（residual），此残差便是计算得到的声带原始激励信号。对残差求自相关，然后利用周期检测函数求出index值。

Feature extraction

​ 特征提取是rnnoise流程中很关键的的一步，送入RNN的数据便是通过这一步获得的，整个特征提取数据流动如下图所示。帧移为10ms，使用采样率为48kHz的音频时，每一帧的新数据长度变为480，即下图中的in，将其与上一次输入的480个点数据拼接得到长度为960的数据，此为一帧。对其做长度为960的傅里叶变换并取前481（代码中此值记作FREQ_SIZE）点数据得到$X$。$Ex$则是利用$X$求出的22个频带的能量。

数据p

​ 上图中有一段长度为1728的数据名为pitch_buf，每次有新数据in进入时，都对其进行更新。对于某一帧数据进行基音追踪得到index，则p=pitch_buf[1728 - WINDOW_SIZE - index : 1728 - WINDOW_SIZE - index]，其中WINDOW_SIZE=960)。即p为长度为960的数据，其最后一个数据为pitch_buf的第1728 - index个数据。同理对p进行傅里叶变换并取FREQ_SIZE个点得到P。

频带能量计算

​ 由$X$得到$Ex$以及由$P$得到$Ep$的做法相同，即将三角滤波组作用于频域数据，所使用的三角滤波器如图所示：

如图所示，上面的图展示了完整的22个三角滤波器，下面是前面若干个滤波器放大后的结果。对于每个频带，滤波时，将滤波器的所有点（图中用小圆圈表示）与频域数据$X$或$P$的对应部分的模长平方相乘相加即可，也即论文中的$E(b)={\sum }_k{\omega }_b(k)|X(k){|}^2$，同时由上图可以看出，每个滤波器的中间点和上一个滤波器的右端及下一个滤波器的左端对应同一个频点。此外还需要注意的是，第一个滤波器和最后一个滤波器只有三角形的一半，因而计算这两个频带的能量时需要乘以2。

BFCC相关数据

​ 对$X$进行三角滤波得到22个频带的能量$Ex$，接着对$Ex$取对数并做DCT变换得到$Ly$，同时利用一个结构（示意图中的环形buffer）对最新8帧的$Ly$进行储存，每当新一帧的数据进入后，都对其进行更新。

​ 将当前帧的$Ly$作为特征的前22个数据，同时对$Ly$进行一阶差分和二阶差分运算并分别取前6个值作为第23-34个特征，再使用最新8帧的$Ly$运算得到平稳性度量（spectral non-stationarity metric），其计算公式为：$snsm[i]=\frac{1}{8}\sum _{n=0}^7\sum _{m=0}^7\Vert Ly[i-n]-Ly[i-m]{\Vert }^2$，$Ly[i]$表示当前帧的能量倒谱系数。此度量是第42个特征。

基音相关数据

​ 对数据$p$进行傅里叶变换得到$P$，并利用$P$和$X$计算得到基音相关度$Exp$，计算公式为： $Exp[b]=\frac{{\sum }_k{\omega }_b(k)\Re [X(k)P^{\ast }(k)]}{\sqrt{{\sum }_k{\omega }_b(k)|X(k){|}^2\cdot {\sum }_k{\omega }_b(k)|P(k){|}^2}}$，$Exp[b]$即论文中的$p_b$

对$Exp$的理解：$p$可以理解为$x$去掉最新的index个新数据再添加index个旧数据得到的（也即$p[n]=x[n-index]$）。而前面所说的$x$是具有短时周期性的，相当于$p$即为$x$平移一个周期得到，因此理论上纯净语音的$p$和$x$的频谱（共振峰）应当近乎相同，因而二者的相关系数较大。很容易联想到，噪声较大时，二者的相关系数会有所降低（因为不同时间的噪声认为是不相关的）。

​ 取$Exp$的前六个数据（对应低频部分）作为第35-40个特征值，然后将前面基音分析（音高追踪）得到的index作为第41个特征。至此42个特征全部得到。

增益计算与插值

​ 将前面计算出的42维特征通过前文所述的RNN网络，输出22维增益，此输出的target即为前面所叙述的$g_b=\sqrt{\frac{E_s(b)}{E_x(b)}}$。对该值的理解：由于噪声和语音不相关，因而带噪语音的能量必然大于纯净语音的能量即$E_s<E_x$。将频域划分为22个频带计算能量，对于每个频带，$E_s$约接近$E_x$说明该频带所含的噪声越少，即语音越纯净，$g_b$越大。对于含噪语音，给每个频带乘以增益，其物理意义即为当该频带噪声较大时，将其乘以一个较小的增益，反之则乘以一个较大的增益，这样便可以增强语音、抑制噪声。

​ 注意$X$的长度维481，而$g$的长度为22，因而需要对其进行插值，使其长度变为481从而可以和$X$逐点相乘，插值方式如下：第$b$个频带的增益为$g_b$，假设此频带对应三角滤波器点数为4，即$[0,0.5,1,0,5]$（与事实上低频区域的三角滤波器一致），则此频带插值后的增益为$[0,0.5g_b,g_b,0.5g_b]$，对于每个频带做类似操作，便可将其插值到481的长度。

Pitch filtering

​ 由于生成特征的时候使用的时Bark’s scale进行频带处理，在计算频带能量时有相乘求和的操作，这样就会带来一定的平滑（smoothing）效果，使得共振峰凸显效果变弱，因而希望使用一个pitch滤波器对共振峰进行加强。滤波操作实际上便是$X(k)+{\alpha }_{b}P(k)$。根据前文所述，$X(k)$和$P(k)$有着相似的语音共振峰，而两者的噪声谱的相关程度就相对弱一些，因而此求和操作可以加强共振峰，如下图所示，红色为进行pitch滤波后的频谱，可以看出在共振峰处有加强效果。${\alpha }_b$的选取方法在前面论文介绍的pitch filtering部分有详细介绍。当然在求出${\alpha }_b$后，也需要像$g$一样进行插值后再与$P$相乘。滤波完成后进行逆傅里叶变换便可得到降噪后的数据。

一些值得注意的问题

1. 高频数据消失

​ 如上图所示为某含噪语音频谱及其经过rnnoise降噪后的频谱，在实现降噪的同时，可以看出上面高频区域的数据完全消失，可能的原因如下：观察三角滤波器的图我们不难发现，其最右边的滤波器的频点也只是到了400左右，远远未达到480，因此这后面的高频区域没有被三角滤波器所覆盖。最终的增益插值后，在该部分也全部为0，因而乘以增益后，该部分的频谱即为绝对的0，即为图片所呈现的样子。

原项目代码分析

​ 主要的降噪代码位于denoise.c中，下面对代码进行部分分析，主要以介绍函数功能为主，部分运算会被删去。

#define FRAME_SIZE_SHIFT 2
#define FRAME_SIZE (120<<FRAME_SIZE_SHIFT)      // FRAME_SIZE = 480
#define WINDOW_SIZE (2*FRAME_SIZE)              // WINDOW_SIZE = 960
#define FREQ_SIZE (FRAME_SIZE + 1)              // FREQ_SIZ = 481

#define PITCH_MIN_PERIOD 60
#define PITCH_MAX_PERIOD 768
#define PITCH_FRAME_SIZE 960
#define PITCH_BUF_SIZE (PITCH_MAX_PERIOD+PITCH_FRAME_SIZE)   // PITCH_BUF_SIZE = 1728

#define SQUARE(x) ((x)*(x))

#define NB_BANDS 22

#define CEPS_MEM 8
#define NB_DELTA_CEPS 6

#define NB_FEATURES (NB_BANDS+3*NB_DELTA_CEPS+2)       //NB_FEATURES = 42


#ifndef TRAINING
#define TRAINING 0
#endif

// 22个点，opus band相关，和频率有一定对应关系，用来进行三角滤波
static const opus_int16 eband5ms[] = {
/*0  200 400 600 800  1k 1.2 1.4 1.6  2k 2.4 2.8 3.2  4k 4.8 5.6 6.8  8k 9.6 12k 15.6 20k*/
        0,  1,  2,  3,  4,  5,  6,  7,  8, 10, 12, 14, 16, 20, 24, 28, 34, 40, 48, 60, 78, 100
};


typedef struct {
    int init;
    kiss_fft_state *kfft;
    float half_window[FRAME_SIZE];
    float dct_table[NB_BANDS*NB_BANDS];
} CommonState;

struct DenoiseState {
    float analysis_mem[FRAME_SIZE];
    float cepstral_mem[CEPS_MEM][NB_BANDS];     //[8][22]
    int memid;
    float synthesis_mem[FRAME_SIZE];
    float pitch_buf[PITCH_BUF_SIZE];        //[1728]
    float pitch_enh_buf[PITCH_BUF_SIZE];    //[1728]
    float last_gain;
    int last_period;
    float mem_hp_x[2];
    float lastg[NB_BANDS];
    RNNState rnn;
};

/* 
 * 计算每个band的能量
 * kiss_fft_cpx就是一个复数的结构体，包含r和i，数据类型是float，这里的参数实际上就是一个复数数组
 * bandE是储存每个band的能量的，有22个值
 */
void compute_band_energy(float *bandE, const kiss_fft_cpx *X) {
    int i;
    float sum[NB_BANDS] = {0};
    for (i=0;i<NB_BANDS-1;i++)        //对22个band进行操作
    {
        int j;
        int band_size;
        //band的点数为eband5ms每两个值之间的差乘以4
        band_size = (eband5ms[i+1]-eband5ms[i])<<FRAME_SIZE_SHIFT;    
        for (j=0;j<band_size;j++) {                                   // 每个band内部的计算
            float tmp;
            // frac实际上就是论文中的wb，将frac作用于X实际上就是将opusband作用于频谱
            float frac = (float)j/band_size;     
            tmp = SQUARE(X[(eband5ms[i]<<FRAME_SIZE_SHIFT) + j].r);   
            tmp += SQUARE(X[(eband5ms[i]<<FRAME_SIZE_SHIFT) + j].i);
            // 每两个三角窗都有一半区域是重叠的，因而对于某一段频点，其既被算入该频带，也被算入下一频带
            sum[i] += (1-frac)*tmp;
            sum[i+1] += frac*tmp;        
        }
    }
    sum[0] *= 2;             // 第一个band和最后一个band的窗只有一半因而能量乘以2
    sum[NB_BANDS-1] *= 2;
    for (i=0;i<NB_BANDS;i++)
    {
        bandE[i] = sum[i];
    }
}
/*
 * 计算band之间的相关度
 */
void compute_band_corr(float *bandE, const kiss_fft_cpx *X, const kiss_fft_cpx *P) {
    int i;
    float sum[NB_BANDS] = {0};
    for (i=0;i<NB_BANDS-1;i++)
    {
        int j;
        int band_size;
        band_size = (eband5ms[i+1]-eband5ms[i])<<FRAME_SIZE_SHIFT;
        for (j=0;j<band_size;j++) {
            float tmp;
            float frac = (float)j/band_size;        // 求解wb
            tmp = X[(eband5ms[i]<<FRAME_SIZE_SHIFT) + j].r * P[(eband5ms[i]<<FRAME_SIZE_SHIFT) + j].r;
            tmp += X[(eband5ms[i]<<FRAME_SIZE_SHIFT) + j].i * P[(eband5ms[i]<<FRAME_SIZE_SHIFT) + j].i;
            sum[i] += (1-frac)*tmp;
            sum[i+1] += frac*tmp;
        }
    }
    sum[0] *= 2;
    sum[NB_BANDS-1] *= 2;
    for (i=0;i<NB_BANDS;i++)
    {
        bandE[i] = sum[i];
    }
}
/*
 * 插值，即22点数据插值得到481点数据
 */
void interp_band_gain(float *g, const float *bandE) {
    // 省略
}


CommonState common;

static void check_init() { // 检查初始化状态，即要对fft运算分配内存空间，然后生成需要使用的dct table
    // 省略
}

// 对输入数据 in 做离散余弦变换
static void dct(float *out, const float *in) {       
    //省略
}

//对输入数据进行快速傅里叶变换，并将结果储存在out中，注意out的长度为481
static void forward_transform(kiss_fft_cpx *out, const float *in) {    // out = fft(in)
    int i;
    kiss_fft_cpx x[WINDOW_SIZE];
    kiss_fft_cpx y[WINDOW_SIZE];
    check_init();          // 生成fft的存储数据
    for (i=0;i<WINDOW_SIZE;i++) {        //将输入数据赋值给x
        x[i].r = in[i];
        x[i].i = 0;
    }
    opus_fft(common.kfft, x, y, 0);   // 对输入数据x进行fft，计算结果给y
    for (i=0;i<FREQ_SIZE;i++) {
        out[i] = y[i];                    // 储存计算结果
    }
}

// 对输入信号in做ifft
static void inverse_transform(float *out, const kiss_fft_cpx *in) {
   //省略
}

static void apply_window(float *x) {    //给输入数据x进行加窗
    //省略
}

// 获取denoise的大小（字节数）
int rnnoise_get_size() {     
    return sizeof(DenoiseState);
}

//获取帧长 480
int rnnoise_get_frame_size() {    
    return FRAME_SIZE;
}

// 对DenoiseState整个结构体进行初始化
int rnnoise_init(DenoiseState *st, RNNModel *model) {  
    //省略
}

// 创建一个Denoise实例并对其进行初始化
DenoiseState *rnnoise_create(RNNModel *model) { 
    //省略
}

//销毁数据
void rnnoise_destroy(DenoiseState *st) {    
    //省略
}


// 帧分析 ：in是长度为480的输入数据，对x做傅里叶变换得到X并求出X的频带能量Ex
static void frame_analysis(DenoiseState *st, kiss_fft_cpx *X, float *Ex, const float *in) {
    int i;
    float x[WINDOW_SIZE];
    // 从 analysis_mem拷贝FRAME_SIZE个数据赋给x的前半段，analysis_mem是上一次的输入
    RNN_COPY(x, st->analysis_mem, FRAME_SIZE);
    for (i=0;i<FRAME_SIZE;i++) x[FRAME_SIZE + i] = in[i];   // 将输入数据赋给x的后半段
    RNN_COPY(st->analysis_mem, in, FRAME_SIZE);    // 然后再将in拷贝给 analysis_mem
    apply_window(x);               // 对x进行加窗
    forward_transform(X, x);      // 对x进行fft并储存在X中

    compute_band_energy(Ex, X);      // 计算频带能量
}

/*
 * 计算42维特征
 */
static int compute_frame_features(DenoiseState *st, kiss_fft_cpx *X, kiss_fft_cpx *P,
                                  float *Ex, float *Ep, float *Exp, float *features, const float *in) {

    int i;
    float E = 0;
    float *ceps_0, *ceps_1, *ceps_2;
    float spec_variability = 0;
    float Ly[NB_BANDS];
    float p[WINDOW_SIZE];
    float pitch_buf[PITCH_BUF_SIZE>>1];
    int pitch_index;      // 这里设置的是整型，但下面的参数是其地址，则可改变原来的值
    float gain;
    float *(pre[1]);
    float tmp[NB_BANDS];
    float follow, logMax;
    frame_analysis(st, X, Ex, in);    // X是做完fft的数据， Ex储存频带能量
    

    // 也是从源src拷贝给dst n个字节数，不同的是，若src和dst内存有重叠，也能顺利拷贝
    // pitch_buffer长度是1728，这里的意思是将其后面(1728 - 480)个数据放到最前面
    RNN_MOVE(st->pitch_buf, &st->pitch_buf[FRAME_SIZE], PITCH_BUF_SIZE-FRAME_SIZE);
    // 对其后面的480个数据进行更新
    RNN_COPY(&st->pitch_buf[PITCH_BUF_SIZE-FRAME_SIZE], in, FRAME_SIZE);
    
    pre[0] = &st->pitch_buf[0];
    /*
     * 降采样，对pitch_buf平滑降采样，求自相关，利用自相关求lpc系数，然后进行lpc滤波，即得到lpc残差
     */
    pitch_downsample(pre, pitch_buf, PITCH_BUF_SIZE, 1);

    // 寻找基音周期   存入pitch_index   
    pitch_search(pitch_buf+(PITCH_MAX_PERIOD>>1), pitch_buf, PITCH_FRAME_SIZE,
                 PITCH_MAX_PERIOD-3*PITCH_MIN_PERIOD, &pitch_index);
    pitch_index = PITCH_MAX_PERIOD-pitch_index;
    // 去除高阶谐波影响
    gain = remove_doubling(pitch_buf, PITCH_MAX_PERIOD, PITCH_MIN_PERIOD,
                           PITCH_FRAME_SIZE, &pitch_index, st->last_period, st->last_gain);    
    

    st->last_period = pitch_index;
    st->last_gain = gain;
    // 根据index得到p[i]
    for (i=0;i<WINDOW_SIZE;i++)
        p[i] = st->pitch_buf[PITCH_BUF_SIZE-WINDOW_SIZE-pitch_index+i];
    apply_window(p);                          // pitch数据应用window
    forward_transform(P, p);                  // 对pitch数据进行傅里叶变换
    compute_band_energy(Ep, P);               // 计算pitch部分band能量
    compute_band_corr(Exp, X, P);             // 计算信号频域与pitch频域的相关band能量系数
   
    for (i=0;i<NB_BANDS;i++) Exp[i] = Exp[i]/sqrt(.001+Ex[i]*Ep[i]);     // Exp进行标准化
    dct(tmp, Exp);                          // 然后再做一次dct，实际上就是信号与pitch相关BFCC了

    for (i=0;i<NB_DELTA_CEPS;i++) features[NB_BANDS+2*NB_DELTA_CEPS+i] = tmp[i];
    features[NB_BANDS+2*NB_DELTA_CEPS] -= 1.3;
    features[NB_BANDS+2*NB_DELTA_CEPS+1] -= 0.9;
    features[NB_BANDS+3*NB_DELTA_CEPS] = .01*(pitch_index-300);
    //而feature的1-NB_BANDS（22）是由log10(Ex)再做一次DCT后填充的，代码如下
    logMax = -2;
    follow = -2;
    for (i=0;i<NB_BANDS;i++) {
        Ly[i] = log10(1e-2+Ex[i]);
        Ly[i] = MAX16(logMax-7, MAX16(follow-1.5, Ly[i]));
        logMax = MAX16(logMax, Ly[i]);
        follow = MAX16(follow-1.5, Ly[i]);
        E += Ex[i];
    }
    if (!TRAINING && E < 0.04) {
        /* If there's no audio, avoid messing up the state. */
        RNN_CLEAR(features, NB_FEATURES);
        return 1;
    }

    for (int j=0;j<22;j++) {  //求y的平方和
        //printf("%f----%d\n", Ly[j],j+1);
    }
    dct(features, Ly);
    //----
    /*
     * cepstral_mem是一个8*22的数组，每一次feature里的值填充到ceps_0,然后这个数组会往下再做一次。
     * ceps_0是float指针，它指向的是ceptral_mem第一个NB_BANDS数组，然后每次与相邻的band数组相见，做出一个delta差值。
     */
    features[0] -= 12;
    features[1] -= 4;
    ceps_0 = st->cepstral_mem[st->memid];
    ceps_1 = (st->memid < 1) ? st->cepstral_mem[CEPS_MEM+st->memid-1] : st->cepstral_mem[st->memid-1];
    ceps_2 = (st->memid < 2) ? st->cepstral_mem[CEPS_MEM+st->memid-2] : st->cepstral_mem[st->memid-2];
    for (i=0;i<NB_BANDS;i++) ceps_0[i] = features[i];
    st->memid++;



    for (i=0;i<NB_DELTA_CEPS;i++) {    //
        features[i] = ceps_0[i] + ceps_1[i] + ceps_2[i];
        features[NB_BANDS+i] = ceps_0[i] - ceps_2[i];
        features[NB_BANDS+NB_DELTA_CEPS+i] =  ceps_0[i] - 2*ceps_1[i] + ceps_2[i];
    }
    /* Spectral variability features. */
    if (st->memid == CEPS_MEM) st->memid = 0;
    //----
    // 最后一个特性值
    for (i=0;i<CEPS_MEM;i++)
    {
        int j;
        float mindist = 1e15f;
        for (j=0;j<CEPS_MEM;j++)
        {
            int k;
            float dist=0;
            for (k=0;k<NB_BANDS;k++)
            {
                float tmp;
                tmp = st->cepstral_mem[i][k] - st->cepstral_mem[j][k];
                dist += tmp*tmp;
            }
            if (j!=i)
                mindist = MIN32(mindist, dist);
        }
        spec_variability += mindist;
    }
    //应该是最后一个特征谱稳度
    features[NB_BANDS+3*NB_DELTA_CEPS+1] = spec_variability/CEPS_MEM-2.1;     
                // 返回值为判断是否有语音

    for(int u = 0; u < 42; u++){
        //printf("%f ", P[u].r);
    }
    return TRAINING && E < 0.1;
}

// 将处理好的频域数据处理得到时域数据
static void frame_synthesis(DenoiseState *st, float *out,  kiss_fft_cpx *y) {
    float x[WINDOW_SIZE];
    int i;
	
    inverse_transform(x, y);
    apply_window(x);
    for (i=0;i<FRAME_SIZE;i++) out[i] = x[i] + st->synthesis_mem[i];
    RNN_COPY(st->synthesis_mem, &x[FRAME_SIZE], FRAME_SIZE);
}
/*
 * biquad滤波，双二阶滤波器
 * y为输出信号
 * x为输入信号
 * a b 为滤波器系数
 * mem为调节参数，从训练模型中获取
 * N为信号长度
 */
static void biquad(float *y, float mem[2], const float *x, const float *b, const float *a, int N) {
	//略过
}

/*
 * 论文中的pitch filter部分，下面的r[i]实际上就是论文中的alpha
 */
void pitch_filter(kiss_fft_cpx *X, const kiss_fft_cpx *P, const float *Ex, const float *Ep,
                  const float *Exp, const float *g) {
  // 略过
}

//对每一帧及进行处理
float rnnoise_process_frame(DenoiseState *st, float *out, const float *in) {
    for (int j = 0; j < 2; j++) {
        //printf("%f---%d\n", st->mem_hp_x[j],j+1);
    }
    int i;
    kiss_fft_cpx X[FREQ_SIZE];
    kiss_fft_cpx P[WINDOW_SIZE];
    float x[FRAME_SIZE];
    float Ex[NB_BANDS], Ep[NB_BANDS];
    float Exp[NB_BANDS];
    float features[NB_FEATURES]；
    float g[NB_BANDS];
    float gf[FREQ_SIZE]={1};
    float vad_prob = 0;
    int silence;
    static const float a_hp[2] = {-1.99599, 0.99600};
    static const float b_hp[2] = {-2, 1};
    //先进行双二阶滤波
    biquad(x, st->mem_hp_x, in, b_hp, a_hp, FRAME_SIZE);
    // 求出特征值并计算返回vad（利用能量阈值判断）
    silence = compute_frame_features(st, X, P, Ex, Ep, Exp, features, x);
	//若为非静默，则进行处理
    if (!silence) {
        //用rnn网络得到增益
        compute_rnn(&st->rnn, g, &vad_prob, features);
        //pitch滤波
        pitch_filter(X, P, Ex, Ep, Exp, g);
        for (i=0;i<NB_BANDS;i++) {
            float alpha = .6f;
            g[i] = MAX16(g[i], alpha*st->lastg[i]);
            st->lastg[i] = g[i];
        }
        //对增益进行插值
        interp_band_gain(gf, g);
#if 1
        //频谱乘以增益
        for (i=0;i<FREQ_SIZE;i++) {
            X[i].r *= gf[i];
            X[i].i *= gf[i];
        }
#endif
    }
	//合成时域数据
    frame_synthesis(st, out, X);

    return vad_prob;

}