576样本的PCM信号加窗之疑惑
从WAV文件的语音数据中读出1152*2个语音数据, 对应到一个音框(frame). 这个1152*2个语音数据在一个frame中是有如图的分析
分组, 分声道处理的伪代码如下:
/* 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 )
*/
}
苦无同道, 有达人赐教不甚感激!