576样本的PCM信号加窗之疑惑

 

从WAV文件的语音数据中读出1152*2个语音数据, 对应到一个音框(frame). 这个1152*2个语音数据在一个frame中是有如图的分析

 


 

数字输出python_filter

分组, 分声道处理的伪代码如下:

/* polyphase filtering 多相滤波, 对时域信号进行滤波*/ 
  
for(gr=0; gr<2; gr++)       // 组 group 
  
{ 
  
    for(channel=config.wave.channels; channel--; )  // 声道 channel 
  
    { 
  
        /* 每次处理576个样本 18 == ( 576 / 32 ) */ 
  
        for(i=0; i<18; i++ )  
  
        { 
  
            /* 
  
             * 对18循环是什么意思啊? 18 == (( 1152/2) / 32) ) 
  
             * 每次处理32个样本,对于buff[gr][channel]就需要处理18次, 因为是576个样本 
  
             * 在最内层循环, 也就是这个18次循环里头. 传递给滤波器的是576笔数据 
  
             */ 
  
             /*  
  
             * 对PCM信号buffer_window[channel]加窗,  
  
             * 得到加窗后的信号为win_que[channel][0]  
  
             */ 
  
             window_subband(&buffer_window[channel], &win_que[channel][0], channel); 
  
             /* 
  
              * 对加窗信号win_que进行滤波, 得到子带信号sb_sample 
  
              * sb_sample[2][3][18][SBLIMIT] 
  
              */ 
  
              filter_subband(&win_que[channel][0], &sb_sample[channel][gr+1][i][0]); 
  
        } 
   
    } 
  
} 
  
 
  
      
  
short int buffer_window [channel][ samp_per_frame] 
  
是用来存储PCM信号的, channel一般为2, 表示声道; samp_per_frame为576 
  
 
需要细分析的就是 
   
        for(i=0; i<18; i++ )  
  
        { 
  
            /* 
  
             * 对18循环是什么意思啊? 18 == (( 1152/2) / 32) ) 
  
             * 每次处理32个样本,对于buff[gr][channel]就需要处理18次, 因为是576个样本 
  
             * 在最内层循环, 也就是这个18次循环里头. 传递给滤波器的是576笔数据 
  
             */ 
  
             /*  
  
             * 对PCM信号buffer_window[channel]加窗,  
  
             * 得到加窗后的信号为win_que[channel][0]  
  
             */ 
  
             window_subband(&buffer_window[channel], &win_que[channel][0], channel); 
  
             /* 
  
              * 对加窗信号win_que进行滤波, 得到子带信号sb_sample 
  
              * sb_sample[2][3][18][SBLIMIT] 
  
              */ 
  
              filter_subband(&win_que[channel][0], &sb_sample[channel][gr+1][i][0]); 
  
        } 
  
如图2.10 
   
对于最内层循环, 即真正的处理PCM信号的那一层. 每一次循环都会将&buffer_window[channel]取出32个样本. 在对x加窗的时候我一直不知道这32个样本如果存储到x中去, 一直认为是 
  
直接平移到x中, FIFO的方式平移 
  
buffer_window[channel][0]-- buffer_window[channel][31]平移到x[0]—x[31] 
  
buffer_window[channel][32]-- buffer_window[channel][63]平移到x[32]—x[63] 
  
 . . . . 
  
如果是这样做的话, 那么PCM还有64个样本无法处理. 
  
FK, 我错了! 
  
 
  
      
  
加窗为什么要以这样的方式? 
  
 
  
      
  
Shine的处理方式是 
   
static int      off[ 2 ] = {0, 0};          // ??? 
   
/** 
  
 * PURPOSE:  Overlapping window on PCM samples 
  
 * SEMANTICS: 
  
 * 32 16-bit pcm samples are scaled to fractional 2's complement and 
  
 * concatenated to the end of the window buffer #x#. The updated window 
  
 * buffer #x# is then windowed by the analysis window #enwindow# to produce the 
  
 * windowed sample #z# 
  
 *  
  
 * 给PCM信号buffer加窗, 这个PCM信号具有576个样本. 输出加窗后的信号为z. 
  
 * 这个加窗函数会被循环调用18次, 因为每个样本空间具有576个样本值, 而加窗函数每次只处理32个样本, 
  
 * 所以必须连续处理18(576/32)次才能处理完成. 处理的顺序是  
 * 0-31, 511-480, 479-448, ... 63-32, 31-0, 511-480.(刚好18次) 
  
 * 发现对于buffer—x是如此的操作 
  
* x[31] = buffer[0] . . .  x[0] = buffer[31]; 
  
* x[511] = buffer[32] . . .  x[480] = buffer[63]; 
  
* . . . 
  
* x[63] = buffer[480] . . .  x[32] = buffer[511]; 
  
* x[31] = buffer[512] . . .  x[0] = buffer[543];       // 覆盖了一次 
  
* x[511] = buffer[544] . . .  x[480] = buffer[575]; // 覆盖了一次 
  
* 
  
* 这里是所无法理解的. 
  
 *  
  
 * **buffer     :   输入PCM信号 ---- 576个样本, 但是每执行一次该函数并不是对576个样本全加窗 
  
 * z[HAN_SIZE]      :   加窗后的信号 ----  
  
 * k                :   channels, 对左声道或右声道加窗. 
  
 */ 
  
void 
  
window_subband(short int **buffer, double z[HAN_SIZE], int k) 
  
{ 
  
    /* replace 32 oldest samples with 32 new samples*/ 
  
    int i; 
  
    for (i=0; i<32; i++) 
  
    { 
  
    /* 
  
     * #define SCALE 32768 -- 0x8000  将信号幅度值左移了8位. 
  
     *  
  
     * x[k][index] 在第2维是否永远只有前32个元素是有效的? 
  
     * 我觉得奇怪, 后面的(HAN_SIZE-32)个元素都是空的. 
  
     * 以上观点是错误的, 因为之前我没有注意到off[]的存在以及它的含义, 呵呵! 
  
     */ 
  
        x[ k ][ 31-i+off[k] ] = (double)( *( (*buffer)++ ) ) / SCALE; 
  
    } 
  
    /*  
  
     * shift samples into proper window positions 
  
     * 加窗 
  
     * HAN_SIZE == 512 
  
     *  
  
     * 这里的窗是怎么加的啊, 好奇怪啊 
  
     *  
  
     * 明白了, x信号采取右移的方式, 移32个单位 
  
     * 一共会加18次窗 
  
     *  
  
     *         | enwindow    0       1      2      3   ... 511 
  
     * -------|----------------------------------------------------------- 
  
     *       0 | x[][i]      0       1      2      3   ... 511 
  
     *       1 | x[][i]      480     481     482 483 ...479 
  
     *         |                     . 
  
     *         |                 . 
  
     *         |                     . 
  
     *       17| x[][i]     511       510     509     508 ... 0    
  
     *  
  
     */ 
  
    for (i=HAN_SIZE; i--; ) 
  
    { 
  
        z[ i ] = x[ k ][ (i+off[k]) & (HAN_SIZE-1) ] * enwindow[ i ]; 
  
    } 
  
     
  
    off[k] += 480;           /* offset is modulo (HAN_SIZE) */ 
  
    off[k] &= (HAN_SIZE - 1); 
  
} 
  
 
  
      
  
 
  
      
  
 
  
      
  
 
  
      
  
 
  
      
  
// ----------- 加窗以及滤波实现如下 ------------- 
   
/** 
  
 *  
  
 * Frame -- 台湾资料上把这个称为音框, 不知道我们怎么称 
  
 *  
  
 * 假令是2声道16位  
  
 * Q: CBR/VBR?  A: CBR 
  
 *  
  
 * 对于MP3来说, 一个音框有2组数据, 每组数据有1152笔取样. 那么每个声道就有(1152/2)=576笔样值. 
  
 *  
  
 *  
  
 */ 
  
#include <math.h> 
  
#include <stdlib.h> 
  
 
  
      
  
#include "types.h" 
  
#include "tables.h"  
  
#include "l3subband.h" 
  
 
  
      
  
// ------------------------------------------------------------------------------- 
  
static int      off[ 2 ] = {0, 0};          // ??? 
  
/* 
  
 * 存储PCM信号, x[0]是左声道, x[1]是右声道. 读一次WAVE流会读1152个样本点. 
  
 * 那么左右声道各512个采样点. 
  
 */ 
  
static double  x[ 2 ][ HAN_SIZE ]; 
  
/* 
  
 * 子带滤波器阵列的系数, 由32个子带滤波器组成一个滤波器阵列.  
  
 * filter[0][64]    第1个滤波器系数 
  
 * filter[1][64]    第2个滤波器系数 
  
 * filter[2][64]    第3个滤波器系数 
  
 *      .               . 
  
 *      .               . 
  
 *      .               . 
  
 * filter[31][64]   第31个滤波器系数 
  
 * 对于一个子带滤波器而言, 有64个系数, 为什么是64个系数呢?这里我觉得奇怪 
  
 * 按照规范是512个, 但是每隔64个就出现一次循环, 所以只需要64个就够了. 
  
 * 假设前64个称为A 
  
 * 则第2个64个为-A 
  
 * 则第3个64个为A 
  
 * 则第4个64个为-A 
  
 *  
  
 */ 
  
static double  filter[ SBLIMIT ][ 64 ];    // 32个频带, 每个频带64个系数 
  
 
  
      
  
// ----------------------------------------------------------------------------- 
  
/* create_ana_filter */ 
  
/** 
  
 * PURPOSE:  Calculates the analysis filter bank coefficients 
  
 * SEMANTICS: 
  
 * Calculates the analysis filterbank coefficients and rounds to the 
  
 * 9th decimal place accuracy of the filterbank tables in the ISO 
  
 * document.  The coefficients are stored in #filter# 
  
 *  
  
 * 计算 各个子带滤波器的系数.  
  
 */ 
  
void 
  
subband_initialise( void ) 
  
{ 
  
    int i, j, k; 
  
     
  
    /*  
  
     * 初始化输入信号. 两个声道x[0]--左声道, x[1]--右声道 
  
     * PCM信号会存储到x中 
  
     */ 
  
    for(i=2; i--; ) 
  
        for(j=HAN_SIZE; j--; ) 
  
            x[i][j] = 0; 
  
     
  
    /*  
  
     * 余弦调变, 用来载到低通滤波器上形成带通滤波器. PCM信号通过的就是 
  
     * 调变后的带通滤波器. 
  
     * 32个频带, 需要32个子频带滤波器组成一组滤波器 [0, 31] 
  
     */ 
  
    for (i=32; i--; ) 
  
    { 
  
        /*  
  
         * 按道理k是需要循环512次, 注意到余弦调变的式子 cos( PI*(2*i+1)*(k-16)/64 ) 
  
         * 是一个周期性表达式, 64为半周期, 因此只需要计算64次就够了. 
  
         */ 
  
         for (k=64; k--; )  
  
         { 
  
             /*  
  
              * filter[i][k]*(1E9) 是为了希望精确度更高.实际上我认为 filter[i][k]*(1E8) 
  
              * 也是可以的, 不知道这里乘上1E9是规范定的还是作者根据自己经验定的. 
  
              */ 
  
             if ( ( filter[i][k]=1e9*cos((double)((2*i+1)*(16-k)*PI64)) ) >= 0 ) 
  
             { 
  
                 modf( filter[i][k]+0.5, &filter[i][k] ); 
  
             } 
  
             else 
  
             { 
  
                 modf( filter[i][k]-0.5, &filter[i][k] ); 
  
             } 
  
             filter[i][k] *= 1e-9;   /* 还原(保证了高精确度) */ 
  
        } 
  
    } 
  
} 
  
 
  
      
  
long int logx = 0; 
  
/** 
  
 * PURPOSE:  Overlapping window on PCM samples 
  
 * SEMANTICS: 
  
 * 32 16-bit pcm samples are scaled to fractional 2's complement and 
  
 * concatenated to the end of the window buffer #x#. The updated window 
  
 * buffer #x# is then windowed by the analysis window #enwindow# to produce the 
  
 * windowed sample #z# 
  
 *  
  
 * 给PCM信号buffer加窗, 这个PCM信号具有576个样本. 输出加窗后的信号为z. 
  
 * 这个加窗函数会被循环调用18次, 因为每个样本空间具有576个样本值, 而加窗函数每次只处理32个样本,  
  
 * 所以必须连续处理18(576/32)次才能处理完成. 处理的顺序是   
  
 * 0-31, 511-480, 479-448, ... 63-32, 31-0, 511-480. 
  
 * 发现第512-575这64个样本被抛弃了, 对此我感到很疑惑, 不知道为什么会被抛弃掉. 
  
 *  
  
 *  
  
 * **buffer     :   输入PCM信号 ---- 576个样本, 但是每执行一次该函数并不是对576个样本全加窗 
  
 * z[HAN_SIZE]  :   加窗后的信号 ----  
  
 * k            :   channels, 对左声道或右声道加窗. 
  
 */ 
  
void 
  
window_subband(short int **buffer, double z[HAN_SIZE], int k) 
  
{ 
  
    /*  
  
     * replace 32 oldest samples with 32 new samples 
  
     */ 
  
    int i; 
  
    for (i=0; i<32; i++) 
  
    { 
  
    /* 
  
     * #define SCALE 32768 -- 0x8000  
  
     * 将信号幅度值左移了8位. 
  
     *  
  
     * x[k][index] 在第2维是否永远只有前32个元素是有效的? 
  
     * 我觉得奇怪, 后面的(HAN_SIZE-32)个元素都是空的. 
  
     * 以上观点是错误的, 因为之前我没有注意到off[]的存在以及它的含义, 呵呵! 
  
     *  
  
     * window_subband会被调用18次, 32*18=576. 576个样本要存储到512个里头,  
  
     * shine采取的策略是舍去后64个样本. 如果舍去前32后32是不是也会一样有效果呢? 
  
     * 也许应该试一试. 
  
     * 处理的顺序是   
  
     * 0-31, 511-480, 479-448, ... 63-32, 31-0, 511-480. 
  
     *  
  
     */ 
  
        x[ k ][ 31-i+off[k] ] = (double)( *( (*buffer)++ ) ) / SCALE; 
  
        /* 
  
         掐头去尾舍去前32后32, 而不是舍去后64. 
  
         */ 
  
    } 
  
    /*  
  
     * shift samples into proper window positions 
  
     * 加窗 
  
     * HAN_SIZE == 512 
  
     *  
  
     * 这里的窗是怎么加的啊, 好奇怪啊 
  
     *  
  
     * 明白了, x信号采取右移的方式, 移32个单位 
  
     * 一共会加18次窗 
  
     *  
  
     *         | enwindow    0       1      2      3   ... 511 
  
     * -------|----------------------------------------------------------- 
  
     *       0 | x[][i]      0       1      2      3   ... 511 
  
     *       1 | x[][i]      480     481     482 483 ... 479 
  
     *         |                     . 
  
     *         |                 . 
  
     *         |                     . 
  
     *       17| x[][i]       511     510     509     508 ... 0    
  
     *  
  
     */ 
  
    for (i=HAN_SIZE; i--; ) 
  
    { 
  
        z[ i ] = x[ k ][ (i+off[k]) & (HAN_SIZE-1) ] * enwindow[ i ]; 
  
    } 
  
     
  
    off[k] += 480;           /* offset is modulo (HAN_SIZE) */ 
  
    off[k] &= (HAN_SIZE - 1); 
  
} 
  
 
  
      
  
/** 
  
 * PURPOSE:  Calculates the analysis filter bank coefficients 
  
 * SEMANTICS: 
  
 *      The windowed samples #z# is filtered by the digital filter matrix #filter# 
  
 * to produce the subband samples #s#. This done by first selectively 
  
 * picking out values from the windowed samples, and then multiplying 
  
 * them by the filter matrix, producing 32 subband samples. 
  
 * 
  
 * 加窗信号#z#通过数字滤波器阵#filter#产生子带信号#s#. 
  
 * 首先会从窗信号中选择 
  
 * 然后用矩阵#filter#去乘它们, 产生32个子带信号. 
  
 *  
  
 * z[HAN_SIZE]  :   加窗后的PCM信号, 该信号具有 512 个采样点. 
  
 * s[SBLIMIT]   :   子带信号 
  
 *  
  
 */ 
  
void 
  
filter_subband(double z[HAN_SIZE], double s[SBLIMIT]) 
  
{ 
  
    doubley[ 64 ]; 
  
   int    i,  j; 
  
    
  
    /* 
  
     * z信号是经过加窗后的信号. 
  
     * y[63] = z[511]+z[447]+z[383]+z[319]+z[255]+z[191]+z[127]+z[63] 
  
     * y[62] = z[510]+z[446]+z[382]+z[318]+z[254]+z[190]+z[126]+z[62] 
  
     *       . 
  
     *       . 
  
     *       . 
  
     * y[0]  = z[448]+z[384]+z[320]+z[256]+z[192]+z[128]+z[64]+z[0] 
  
     *  
  
     * 这里是在提取公因式了,  
  
     * 因为这8个会乘共同的因式cos(......)==filter[i][j] 
  
     */ 
  
   for (i=64; i--; ) 
  
   {        
  
        for (j=8, y[i]=0; j--; ) 
  
        {  
  
            y[i] += z[ i + (j<<6) ]; 
  
        } 
  
    } 
  
   /* 
  
    * 通过余弦调变方式分为32路子带信号. 
  
    * SBLIMIT==32 
  
    */  
  
   for (i=SBLIMIT; i--; ) 
  
   { 
  
         for (j=64, s[i]= 0; j--; ) 
  
         { 
  
           s[i] += filter[i][j] * y[j]; 
  
        } 
  
    } 
  
    /* 
  
     伪代码为: 
  
     int i, j; 
  
     double s[ 32 ]; 
  
     for (i=SBLIMIT; i--; ) 
  
     { 
  
        for(j=512, s[i]=0.0; j--; ) 
  
        { 
  
            s[i] += filter[i][j] * z[j];     
  
        } 
  
     } 
  
     其中filter[i][j] 的数学表达式为 : 
  
     cos( (2*i+1)*(j-16)*PI64 ) 
  
     */ 
  
}

 苦无同道, 有达人赐教不甚感激!