Speex降噪代码分析2——FFT_speex 降噪

 

微信扫一扫,关注公众号“音频算法与工程实践”

FFT基础知识：

1、x(n)是实数，X(k)是复数，实部为偶对称，虚部为奇对称。

2、M点实数(M是偶数，因为代码中使用的N不等于FFT长度，因此此处不使用常见的N表示FFT长度)做FFT，得到结果X，X[0]和X[M/2]都是实数，X[1,...,M/2-1]和X[M/2+1,...,M-1]关于x轴对称，因此结果只保留（M/2+1）点即可。在matlab或者https://octave-online.net/在线执行下面两行代码：X[0]和X[3]是实数，X[1]、X[2]和X[4]、X[5]关于x轴对称。

x=1:6;

fft(x)'

ans =

   21.0000 - 0.0000i

   -3.0000 - 5.1962i

   -3.0000 - 1.7321i

   -3.0000 - 0.0000i

   -3.0000 + 1.7321i

   -3.0000 + 5.1962i

Speex降噪代码的speex_preprocess_state_init()完成FFT初始化，FFT长度M=2N=320。

   st->frame =(spx_word16_t*)speex_alloc(2*N*sizeof(spx_word16_t));

   st->window =(spx_word16_t*)speex_alloc(2*N*sizeof(spx_word16_t));

   st->ft = (spx_word16_t*)speex_alloc(2*N*sizeof(spx_word16_t));

   st->fft_lookup = spx_fft_init(2*N);

preprocess_analysis()函数完成输入信号的交叠(overlap)、加窗、归一化、FFT、计算频域能量ps。spx_fft封装了多种fft代码，通过宏可以随意切换。注意与fft公式的差别是spx_fft结果等价于fft/M这样做的结果是避免FFT结果溢出，可以理解为ifft公式中的除M操作移到fft，因此spx_ifft结果无需再除M。从计算ps的代码来看。FFT结果存放数组元素如下：r[0],r[1],i[1],r[2],i[2],… ,r[N-1],i[N-1],r[N]。总共（N-1）个复数点，2个实数点X[0]和X[N]的虚部0i丢弃，只保存实数分量。st->ft合计所需数组长度2N。https://github.com/xiph/speexdsp/blob/master/libspeexdsp/preprocess.c，代码在614行。

   /* 'Build' input frame */

   for (i=0;i<N3;i++)

      st->frame[i]=st->inbuf[i];

   for (i=0;i<st->frame_size;i++)

      st->frame[N3+i]=x[i];

   /* Update inbuf */

   for (i=0;i<N3;i++)

      st->inbuf[i]=x[N4+i];

   /* Windowing */

   for (i=0;i<2*N;i++)

      st->frame[i] = MULT16_16_Q15(st->frame[i], st->window[i]);

#ifdef FIXED_POINT

   {

      spx_word16_t max_val=0;

      for (i=0;i<2*N;i++)

         max_val = MAX16(max_val, ABS16(st->frame[i]));

      st->frame_shift = 14-spx_ilog2(EXTEND32(max_val));

      for (i=0;i<2*N;i++)

         st->frame[i] = SHL16(st->frame[i], st->frame_shift);

   }

#endif

   /* Perform FFT */

   spx_fft(st->fft_lookup, st->frame, st->ft);

   /* Power spectrum */

   ps[0]=MULT16_16(st->ft[0],st->ft[0]);

   for (i=1;i<N;i++)

      ps[i]=MULT16_16(st->ft[2*i-1],st->ft[2*i-1]) + MULT16_16(st->ft[2*i],st->ft[2*i]);

   for (i=0;i<N;i++)

      st->ps[i] = PSHR32(st->ps[i], 2*st->frame_shift);

从kiss_fft.c的开头注释（Copyright (c) 2003-2004, MarkBorgerding）来看，speex使用的是kiss_fft的早期版本，支持浮点和定点，定点只支持16位#define kiss_fft_scalar spx_int16_t。kiss_fft后续又有更新，新版定点代码支持32位#define kiss_fft_scalar int32_t。正好满足了如今嵌入式芯片ARM和DSP朝32位发展的趋势，比如高通蓝牙芯片QCC5100系列使用了32位QualcommKalimBA DSP，部分电视机ARM芯片集成的音频DSP是32位。因此新算法开发建议下载使用最新的kiss_fft，kiss_fft代码细节将会另辟章节分析。

下面通过测试代码验证前面讲解的spx_fft/spx_ifft两个函数。

//arch.h

#defineFLOATING_POINT

//#defineFIXED_POINT

// preprocess.c

#define FFTLEN 8

void fft_test()

{

       int i;

       void *fft_lookup = spx_fft_init(FFTLEN);

#ifdefFLOATING_POINT

       spx_word16_tin[FFTLEN]={1,2,3,4,5,6,7,8};

       spx_word16_t out[FFTLEN];

    spx_fft(fft_lookup, in, out);

    for (i=0;i<FFTLEN;i++)

    {

        printf("%1.2f ", out[i]);

    }

       printf("\n");

    spx_ifft(fft_lookup, out, in);

    for (i=0;i<FFTLEN;i++)

    {

        printf("%1.0f ", in[i]);

    }

       printf("\n");

#else

       spx_word16_tin[FFTLEN]={10010,10020,10030,10040,10050,10060,10070,10080,10090,10100};

       spx_word16_t out[FFTLEN];

    spx_fft(fft_lookup, in, out);

    for (i=0;i<FFTLEN;i++)

    {

        printf("%d ", out[i]);

    }

       printf("\n");

    spx_ifft(fft_lookup, out, in);

    for (i=0;i<FFTLEN;i++)

    {

        printf("%d ", in[i]);

    }

       printf("\n");

#endif

       spx_fft_destroy(fft_lookup);

}

浮点结果如下：

4.5000 -0.5000 1.2071 -0.5000 0.5000 -0.5000 0.2071 -0.5000

1 2 3 4 5 6 7 8

定点结果如下：

10044 -5 12 -5 5 -5 2 -5

10009 10019 10029 10039 10049 10059 10069 10079

对应matlab代码如下：

更多内容，请关注微信公众号“音频算法与工程实践”。