多频带滤波器移植

 

发了狠了， 买过一个扳子。

多频带滤波器移植

 

// ------------------------------ main.c ------------------------------

#include <stdio.h>

#include <stdlib.h>

#include <reg51.h>

#include "l3subband.h"

 

#define uchar unsigned char

#define uint  unsigned int

#define ulong unsigned long

 

//L3_window_subband(&buffer_window[channel], &win_que[channel][0], channel);

//L3_filter_subband(&win_que[channel][0], &l3_sb_sample[channel][gr+1][i][0]);

 

static short    buffer[2][1152]; // 存储PCM数据, 该数据直接来自WAV文件.[2][1152];

short           *buffer_window[2];// buffer是未加窗信号. 在经过心理模型逻辑的时候会对PCM数据加Hann Window.

 

double          win_que[2][HAN_SIZE];//

double          l3_sb_sample[2][3][18][SBLIMIT];// 经过多相滤波后的子带信号

 

//char code dx516[3] _at_ 0x003b;//这是为了仿真设置的

void main( void )

{  

    uchar gr, channel, i;

    unsigned int j;

    /*

    SCON = 0x52;

    TMOD = 0x20;

    TH1 = 0xf3;

    TR1 = 1;

    */

 

    /* 心里模型初始化. */

    //L3_psycho_initialise( );

    /* 初始化32个子频带滤波器. PCM信号穿过32个频带, 子滤波器的冲激响应 */

    L3_subband_initialise( );

    /* 初始化MDCT */

    //L3_mdct_initialise( );

    /* 为外循环和内循环初始化数据, loop是压缩过程 */

    //L3_loop_initialise( );

 

        // 造假数据

       for(j=0; j<1152; j++)

       {

           buffer[0][j] = j;

           buffer[1][j] = j;

       }

       buffer_window[0] = buffer[0];

       buffer_window[1] = buffer[1];

 

       // 一个Frame有2个颗粒, 排列方式为 gr0, gr1

       for(gr=0; gr<2; gr++)

        //for(gr=0; gr<config.mpeg.mode_gr; gr++)

        {

        // 一个颗粒有两个channel, 排列方式为channel0, channel1

            for(channel=0; channel<2; channel++)

           //for(channel=0; channel<config.wave.channels; channel++)

            {

              /*

               * 一个 channel 具有 576个采样点.

               * 根据规范576个采样点会穿过多相滤波器生成32路子频带信号.

               *

               * 每次处理576个样本 18 == ( 576 / 32 )

               * 18次循环针对一个channel, 循环完毕则处理完毕一个颗粒

               */

                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[2][HAN_SIZE] = win_que[2][512]

                   * 得到加窗后的信号为win_que[channel][0]

                   */

                  L3_window_subband(&buffer_window[channel], &win_que[channel][0], channel);

                   

                    /*

                     * 对加窗信号win_que进行滤波, 得到子带信号l3_sb_sample

                     * l3_sb_sample[2][2][18][SBLIMIT] = l3_sb_sample[2][2][18][32]

                     * l3_sb_sample[2][2][18][32]

                     */

                    L3_filter_subband(&win_que[channel][0], &l3_sb_sample[channel][gr+1][i][0]);

              }

           }

       }/* 多相滤波器结束, 得到了32路子带信号, 时域相卷, 频域相乘 */

 

    while( 1 );

}

 

 

 

// ------------------------------ l3subband.h ------------------------------

#ifndef L3SUBBAND_H

#define L3SUBBAND_H

 

//#include "types.h"

 

// 定义PI, PI4=PI/4, PI64=PI/64

#define PI          3.14159265358979

//#define PI4           PI/4

#define PI64        PI/64

 

#define SBLIMIT 32    // 32路子带信号.

#define HAN_SIZE   512      // Hanning Window

 

 

//#define SCALE_BLOCK    12    //

//#define SCALE_RANGE    64    //

#define SCALE       32768 //

 

void L3_subband_initialise();

void L3_window_subband(short **buffer, double z[HAN_SIZE], int k);

void L3_filter_subband(double z[HAN_SIZE], double s[SBLIMIT]);

 

#endif

 

// ------------------------------ l3subband.c ------------------------------

/*

分析子带滤波器

 

多重相位滤波器

H_i[n] = h[n] * cos( PI64*(2*i+1)*(n-16) )  where i=0~31,n=0~511

其中h[n]是一个低通滤波器,

if(the integer part of n/64 is odd)

   h[n] = -C[n];

else

   h[n] = C[n];

C[n]是MP3规范所定义的一组窗框系数.

 

 

将输入的PCM信号通过32个滤波器(H_i[n], i=0~31)后, 得到P_i[n]

 

P_i[n] = 0;

for(int m=0; m<512; m++)

   P_i[n] += x[n-m] * H_i[m];  where i=0~31

  

最后, 将P_i[n] down-sampling by 32则可得到32个子频带信号.

S_i[n] = P_i[32*n];

for(int m=0; m<512; m++)

   S_i[n] = x[32*n-m]*H_i[m];  where i=0~31

S_0[n]~S_31[n]是滤波器组分析后的结果

 

;;

 

1. 取得512个PCM数据, 不妨令其为pcm[512]

*/

#include <math.h>

#include<string.h>

//#include "types.h"

//#include "tables.h"

#include "l3subband.h"

 

double enwindow[] = {

   0.000000, -0.000000, -0.000000, -0.000000, -0.000000, -0.000000, -0.000000, -0.000001, -0.000001, -0.000001,

   -0.000001, -0.000001, -0.000001, -0.000002, -0.000002, -0.000002, -0.000002, -0.000003, -0.000003, -0.000003,

   -0.000004, -0.000004, -0.000005, -0.000005, -0.000006, -0.000007, -0.000008, -0.000008, -0.000009, -0.000010,

   -0.000011, -0.000012, -0.000014, -0.000015, -0.000017, -0.000018, -0.000020, -0.000021, -0.000023, -0.000025,

   -0.000028, -0.000030, -0.000032, -0.000035, -0.000038, -0.000041, -0.000043, -0.000046, -0.000050, -0.000053,

   -0.000056, -0.000060, -0.000063, -0.000066, -0.000070, -0.000073, -0.000077, -0.000081, -0.000084, -0.000087,

   -0.000091, -0.000093, -0.000096, -0.000099, 0.000102, 0.000104, 0.000106, 0.000107, 0.000108, 0.000109,

   0.000109, 0.000108, 0.000107, 0.000105, 0.000103, 0.000099, 0.000095, 0.000090, 0.000084, 0.000078,

   0.000070, 0.000061, 0.000051, 0.000040, 0.000027, 0.000014, -0.000001, -0.000017, -0.000034, -0.000053,

   -0.000073, -0.000094, -0.000116, -0.000140, -0.000165, -0.000191, -0.000219, -0.000247, -0.000277, -0.000308,

   -0.000339, -0.000371, -0.000404, -0.000438, -0.000473, -0.000507, -0.000542, -0.000577, -0.000612, -0.000647,

   -0.000681, -0.000714, -0.000747, -0.000779, -0.000810, -0.000839, -0.000866, -0.000892, -0.000915, -0.000936,

   -0.000954, -0.000969, -0.000981, -0.000989, -0.000994, -0.000995, -0.000992, -0.000984, 0.000971, 0.000954,

   0.000931, 0.000903, 0.000869, 0.000829, 0.000784, 0.000732, 0.000674, 0.000610, 0.000539, 0.000463,

   0.000379, 0.000288, 0.000192, 0.000088, -0.000021, -0.000137, -0.000260, -0.000388, -0.000522, -0.000662,

   -0.000807, -0.000957, -0.001111, -0.001270, -0.001432, -0.001598, -0.001767, -0.001937, -0.002110, -0.002283,

   -0.002457, -0.002631, -0.002803, -0.002974, -0.003142, -0.003307, -0.003467, -0.003623, -0.003772, -0.003914,

   -0.004049, -0.004175, -0.004291, -0.004396, -0.004490, -0.004570, -0.004638, -0.004691, -0.004728, -0.004749,

   -0.004752, -0.004737, -0.004703, -0.004649, -0.004574, -0.004477, -0.004358, -0.004215, -0.004049, -0.003859,

   -0.003643, -0.003402, 0.003135, 0.002841, 0.002522, 0.002175, 0.001801, 0.001400, 0.000971, 0.000516,

   0.000033, -0.000476, -0.001012, -0.001574, -0.002162, -0.002774, -0.003411, -0.004072, -0.004756, -0.005462,

   -0.006189, -0.006937, -0.007703, -0.008487, -0.009288, -0.010104, -0.010933, -0.011775, -0.012628, -0.013489,

   -0.014359, -0.015234, -0.016113, -0.016994, -0.017876, -0.018757, -0.019634, -0.020507, -0.021372, -0.022229,

   -0.023074, -0.023907, -0.024725, -0.025527, -0.026311, -0.027074, -0.027815, -0.028533, -0.029225, -0.029890,

   -0.030527, -0.031133, -0.031707, -0.032248, -0.032755, -0.033226, -0.033660, -0.034056, -0.034413, -0.034730,

  -0.035007, -0.035242, -0.035435, -0.035586, -0.035694, -0.035759, 0.035781, 0.035759, 0.035694, 0.035586,

   0.035435, 0.035242, 0.035007, 0.034730, 0.034413, 0.034056, 0.033660, 0.033226, 0.032755, 0.032248,

   0.031707, 0.031133, 0.030527, 0.029890, 0.029225, 0.028533, 0.027815, 0.027074, 0.026311, 0.025527,

   0.024725, 0.023907, 0.023074, 0.022229, 0.021372, 0.020507, 0.019634, 0.018757, 0.017876, 0.016994,

   0.016113, 0.015234, 0.014359, 0.013489, 0.012628, 0.011775, 0.010933, 0.010104, 0.009288, 0.008487,

   0.007703, 0.006937, 0.006189, 0.005462, 0.004756, 0.004072, 0.003411, 0.002774, 0.002162, 0.001574,

  0.001012, 0.000476, -0.000033, -0.000516, -0.000971, -0.001400, -0.001801, -0.002175, -0.002522, -0.002841,

   0.003135, 0.003402, 0.003643, 0.003859, 0.004049, 0.004215, 0.004358, 0.004477, 0.004574, 0.004649,

   0.004703, 0.004737, 0.004752, 0.004749, 0.004728, 0.004691, 0.004638, 0.004570, 0.004490, 0.004396,

   0.004291, 0.004175, 0.004049, 0.003914, 0.003772, 0.003623, 0.003467, 0.003307, 0.003142, 0.002974,

   0.002803, 0.002631, 0.002457, 0.002283, 0.002110, 0.001937, 0.001767, 0.001598, 0.001432, 0.001270,

   0.001111, 0.000957, 0.000807, 0.000662, 0.000522, 0.000388, 0.000260, 0.000137, 0.000021, -0.000088,

   -0.000192, -0.000288, -0.000379, -0.000463, -0.000539, -0.000610, -0.000674, -0.000732, -0.000784, -0.000829,

   -0.000869, -0.000903, -0.000931, -0.000954, 0.000971, 0.000984, 0.000992, 0.000995, 0.000994, 0.000989,

   0.000981, 0.000969, 0.000954, 0.000936, 0.000915, 0.000892, 0.000866, 0.000839, 0.000810, 0.000779,

   0.000747, 0.000714, 0.000681, 0.000647, 0.000612, 0.000577, 0.000542, 0.000507, 0.000473, 0.000438,

   0.000404, 0.000371, 0.000339, 0.000308, 0.000277, 0.000247, 0.000219, 0.000191, 0.000165, 0.000140,

   0.000116, 0.000094, 0.000073, 0.000053, 0.000034, 0.000017, 0.000001, -0.000014, -0.000027, -0.000040,

   -0.000051, -0.000061, -0.000070, -0.000078, -0.000084, -0.000090, -0.000095, -0.000099, -0.000103, -0.000105,

   -0.000107, -0.000108, -0.000109, -0.000109, -0.000108, -0.000107, -0.000106, -0.000104, 0.000102, 0.000099,

   0.000096, 0.000093, 0.000091, 0.000087, 0.000084, 0.000081, 0.000077, 0.000073, 0.000070, 0.000066,

   0.000063, 0.000060, 0.000056, 0.000053, 0.000050, 0.000046, 0.000043, 0.000041, 0.000038, 0.000035,

   0.000032, 0.000030, 0.000028, 0.000025, 0.000023, 0.000021, 0.000020, 0.000018, 0.000017, 0.000015,

   0.000014, 0.000012, 0.000011, 0.000010, 0.000009, 0.000008, 0.000008, 0.000007, 0.000006, 0.000005,

   0.000005, 0.000004, 0.000004, 0.000003, 0.000003, 0.000003, 0.000002, 0.000002, 0.000002, 0.000002,

   0.000001, 0.000001, 0.000001, 0.000001, 0.000001, 0.000001, 0.000000, 0.000000, 0.000000, 0.000000,

   0.000000, 0.000000 };

 

/* 这是个非常好的亮点, off[]的出现避免了搬移数据. 达到循环加窗的目的. */

static int off[2] = {0, 0};

/*

 * 存储PCM信号, x[0]是左声道, x[1]是右声道. 读一次WAVE流会读1152=2*576个样本点.

 * HAN_SIZE = Hann Window Size = 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

 *

 * 呵呵, 注意查看<<MP3编解码系统的研究与开发.kdh>>

 * filter[i][k] 实际就是 M[i][k], 也就是H[n]的cos部分,

 * 注意公式的推导哦：

 * M[i][k] = cos( PI*(2*i+1)*(k-16) / 64 );

 * 令k = k + 64*j, 则

 * M[i][k+64*j] = cos( PI*(2*i+1)*((k+64*j)-16) / 64 );

 * =cos( PI*(2*i+1)*(k-16)/64 + PI*(2*i+1)*j );

 * =cos(PI*(2*i+1)*(k-16)/64)*cos(PI*(2*i+1)*j) - sin(PI*(2*i+1)*(k-16)/64)*sin(PI*(2*i+1)*j)

 * =pow((-1), j) * cos(PI*(2*i+1)*(k-16)/64)

 */

static double filter[SBLIMIT][64];// 32个频带, 每个频带64个系数

/**

 * 初始化滤波器组中的一些系数, 主要是cos部分的初始化.

 */

void L3_subband_initialise( void )

{

    int i, j, k;

   

    /* 初始化输入信号. 两个声道x[0]--左声道, x[1]--右声道

    * PCM信号会存储到x中 */

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

        for(j=0; j<HAN_SIZE; j++)

           x[i][j] = 0.0;    

   //直接set为0就好了，不用循环. 优化.

   //memset (x, 0, sizeof(double)*2*HAN_SIZE);

  

   /* 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#

    *

    * 余弦调变, 用来载到低通滤波器上形成带通滤波器. PCM信号通过的就是

    * 调变后的带通滤波器.

    * 32个频带, 需要32个子频带滤波器组成一组滤波器 [0, 31]

    */

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

   {

      /* 按道理k是需要循环512次, 注意到余弦调变的式子 cos( PI*(2*i+1)*(k-16)/64 )

       * 是一个周期性表达式, 64为半周期, 因此只需要计算64次就够了.

       *

       * H_i[n] = h[n] * cos( PI64*(2*i+1)*(n-16) )  where i=0~31,n=0~511

       */

      for (k=0; k<64; k++)

        {

           /*

            * filter[i][k]*(1E9) 是为了希望精确度更高.实际上我认为 filter[i][k]*(1E8)

            * 也是可以的, 不知道这里乘上1E9是规范定的还是作者根据自己经验定的.

            *

            * 我靠, 这里不是未了高精度哦, 是为了跟cos比较.

            */

            if ( ( filter[i][k] = 1e9*cos((double)((2*i+1)*(16-k)*PI64)) ) >= 0 )

            {

               modf( filter[i][k]+0.5, &filter[i][k] );

           }

           else

            {  /* extern float modf(float num, float *i);

                * 功能：将浮点数num分解成整数部分和小数部分

                * 说明：返回小数部分，将整数部分存入*i所指内存中 */

               modf( filter[i][k]-0.5, &filter[i][k] );

            }

         filter[i][k] *= 1e-9;/* 这里是失去了精度, 而不是提高精度. */

        }

   }

}

 

/**

 * 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加Hann Window, 这个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 L3_window_subband(short **buffer, double z[HAN_SIZE], int k)

{

    int i;

 

    /* replace 32 oldest samples with 32 new samples */

    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个里头?

        * 处理的顺序是 :

        *    x ---> 31-0,   511-480, 479-448, ... 63-32,   31-0,    511-480.

        *               |     |        |

        * buffer---> 0-31, 32-63,   64-95,   ... 480-511, 512-543, 544-575

        */

       x[k][31-i+off[k]] = (double)*(*buffer)++ / SCALE;

    }

   

    /*

     * shift samples into proper window positions

     *

     * 这里有很大的优化空间 .....

     *

     * 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=0; 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 L3_filter_subband(double z[HAN_SIZE], double s[SBLIMIT])

{

   double          y[64];

   register 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=0; i<64; i++)

   {

       for (j=0, y[i] = 0; j<8; j++)

       {

           y[i] += z[i+(j<<6)];

       }

   }

  

    /* 通过余弦调变方式分为32路子带信号. SBLIMIT==32 */

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

   {

       for (j=0, s[i]= 0; j<64; 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 )

   */

}