MPEG-1音频编码原理与实验

MPEG-1 Audio Layer-2简介

MPEG-1是为CD光盘介质定制的视频和音频压缩格式，输入为16位的PCM信号，采样率为32，44.1或 48kHz； 输出为32kbps到384kbps。其音频压缩部分分为三层：MPEG-1 Layer1，MPEG-Layer2以及MPEG-Layer3，高层兼容低层。

在数字电视领域中，通常使用第二层即MPEG-Layer2，采用共同频域和时域掩蔽效应的心理声学模型，并对高、中，低频段的比特分配进行限制，并对比特分配、比例因子，取样进行附加编码。

编码器流程图：

人类听觉系统的感知特性

MPEG-1音频编码采用心理声学模型，计算信号中不可感知的部分，并对这部分数据进行压缩处理或去除。

听觉阈值：

人耳对声音强度的感知不仅与声音的响度级有关，还与频率有关。
心理声学中的等响度曲线如图：

声压级越高，等响度曲线越平坦。人耳对 $3-4kHz$ 的声音感觉最灵敏。蓝色曲线之下为听阈，这部分的声音无法被人耳感知，在编码时可去除。

频域掩蔽：

一个较弱的信号会被另一个相近频率的较强信号 的声音所掩蔽 ，这种特性称为频域掩蔽，其随声压级变化曲线如图：

听力阈随频域掩蔽而改变，其下方区域为不能听见的信号，对这部分信号可以进行去除。对于掩蔽阈下的量化噪声也可不予处理。

临界频带：

人类听觉系统大致等效于一个在0Hz到20KHz频率范围内 由25个重叠的带通滤波器组成的滤波器组。

• 人耳不能区分同一频带内同时发生的不同声音；
• 人耳频带被称为临界频带（critical band）；
• 一个临界频带的带宽单位为1巴克（bark) ；
• 500Hz以下每个临界频带的带宽大约是100Hz，从500Hz起 ，临界频带带宽线性增加。

时域掩蔽：

在时间上相邻的声音之间也存在着掩蔽现象，被成为时域掩蔽。时域掩蔽又可分为：

• 强音掩蔽其后的弱音：超前掩蔽，掩蔽时间约为5-20ms
• 弱音被稍后的强音所掩蔽：滞后掩蔽，掩蔽时间约为50-200ms

根据时域掩蔽效应，可以制定编码策略：编码时将时间上相继的一些样值归并成块，并计算样值的比例因子，根据比例因子的大小分配不同的比特数。

编码原理：

1、设计思想：

（1）整个编码流程图分为两条线
第一条线是输入的PCM码流经过多相滤波器组 分成32个子带信号，经过块形成后，对每个子带数据进行线性量化，对部分量化级别采用颗粒优化以增大压缩比，最后装帧输出。

第二条线则是进行FFT变换，经过心理声学模型去除信号中被掩蔽的部分，提取出比例因子后进行动态比特分配：对人耳听觉敏感的低频成分分配较多比特，对高频噪声分配较少的比特，以增大压缩比。比例因子经过选择后，数据编码装帧输出。

（2）时-频分析的矛盾

在时-频分析中也存在“测不准原理”，即时域分辨力和频域分辨力不可兼得，增大频域分辨力的同时时域分辨力便相应减小，反之亦然。
根据这一矛盾，MPEG-1音频编码采用了一系列方法来尽可能地兼顾时域分辨力和频域分辨力：

• 通过子带分析滤波器组使信号具有高的时间分辨率，确保在短暂冲击信号情况下，编码的声音信号具有足够高的质量。
• 使信号通过FFT运算具有高的频率分辨率，因为掩蔽阈值是从功率谱密度推出来 的。

2、多相滤波器组

多相滤波器组(Polyphase Filter Bank)将PCM样本变换到32个子带的频域信号 。

多相滤波器组的缺点在于等带宽的滤波器组与人类听觉系统的临界频带不对应 ，在低频区域，单个子带会覆盖多个临界频带。在这种情况下，量化比特数不能兼每个临界频带 ：

由于一个单频正弦信号输入可能在两个子带中产生非零信号，这种现象被称为混叠现象。滤波后，相邻子带有频率混叠现象，一个子带中的信号可以影响相邻子带的输出。

3、心理声学模型：

MPEG-1中的心理声学模型有两种：

• 心理声学模型1：计算复杂度低，但对假设用户听不到的部分压缩太严重 。
• 心理声学模型2：提供了适合Layer III编码的更多特征 。

实际模型的复杂度取决于所需要的压缩比。

具体步骤如下：

1. 将样本变换到频域。32个子带信号进行FFT变换以补偿频率分辨力，采用Hann加权减少频域中的边界效应 。
2. 确定子带声压级别。
3. 考虑绝对安静时阈值。
4. 根据掩蔽能力的不同，将信号分为乐音（tune）和非乐音/噪声。
5. 音调和非音调掩蔽成分的消除 ：利用标准中给出的绝对阈值消除被掩蔽成分， 考虑在每个临界频带内，小于0.5Bark的距离中只保留最高功率的成分 。
6. 单个掩蔽阈值的计算：音调成分和非音调成分单个掩蔽阈值根据标 准中给出的算法求得。
7. 全局掩蔽阈值的计算：还要考虑别的临界频带的影响。一个掩蔽信号会对其它频带上的信号产生掩蔽效应。这种掩蔽效应称为掩蔽扩散
8. 每个子带的掩蔽阈值 ：选择出本子带中最小的阈值作为子带阈值
9. 计算每个子带信号掩蔽比(signal-to-mask ratio, SMR) $$SMR=信号能量/掩蔽阈值$$之后将SMR传递给编码单元。

4、码率分配、量化

MPEG-1 Layer1的码率分配如下：
在得到比例因子、比例因子选择信息、比例分配信息以及辅助数据比特数后，可以确定用于样值编码的可用比特数：

码率分配的算法如下：

1. 对每个子带计算噪声-掩蔽比NMR：$$NMR=SMR–SNR(dB)$$其中SNR为信噪比，SMR为信掩比。
2. 对最高NMR的子带分配比特，使获益最大的子带的量化级别增加一级 ’
3. 重新计算分配了更多比特子带的NMR。
4. 循环重复2、3步骤，直到没有比特分配。

循环分配算法的目的是使整帧和每个子带的总噪声—掩蔽比达到最小 。

接着进行量化：每个子带从相同的量化集合中选择 。每个子频带12个连续的样值除以比例因子进行归一化，得到的值用 $X$ 表示 ,进行量化计算：$$AX+B$$
其中 $A$ 和 $B$ 为量化系数 ，根据比特分配的信息得到量化级数，再根据量化级数查表得到 $A$ 和 $B$ 的值。量化表如下：

Layer2的码率分配与Layer1基本一致，其量化在Layer1的基础上进行了改进：

• 根据采样和码率量化，不同子带可以从不同的量化器集合中选择，由此以来某些（高频）子带的比特数可能为0 。
• 对量化级别在3、5、9级时，采用“颗粒” 优化

5、装帧

Layer2帧结构：

• 3 组/帧 x 12个样本/子带 x 32个子带/帧 = 1152个样本/帧

• 每个样本的overhead更少

缩放因子：每个子带的3个组尽可能共用缩放因子 ， Layer 2： 1个/（24/36）个样本

• 1/2/3个缩放因子和缩放因子选择信息(scale factor selection information, SCFSI) （每子带2比特）一起传送 .
• 若缩放因子和下一个只有很小的差别，就只传送大的一个
• 若三个缩放因子差别比较大，则全部传送
  

编码器分析与调试实验

1、程序基本框架

整个编码器程序的结构可以主要分为以下几个模块：

（1）滑动窗口滤波，将码流分成32个子带

{
      int gr, bl, ch;
      /* New polyphase filter
	 Combines windowing and filtering. Ricardo Feb'03 */
      for( gr = 0; gr < 3; gr++ )
	for ( bl = 0; bl < 12; bl++ )
	  for ( ch = 0; ch < nch; ch++ )
	    WindowFilterSubband( &buffer[ch][gr * 12 * 32 + 32 * bl], ch,
				 &(*sb_sample)[ch][gr][bl][0] );
    }

WindowFilterSubband（）为滑动窗口滤波的函数。

（2）将左右声道合并为立体声

 if (frame.actual_mode == MPG_MD_JOINT_STEREO) {
      /* this way we calculate more mono than we need */
      /* but it is cheap */
      combine_LR (*sb_sample, *j_sample, frame.sblimit);
      scale_factor_calc (j_sample, &j_scale, 1, frame.sblimit);
    }

（3）计算一帧的比例因子、计算比例因子选择信息

 scale_factor_calc (*sb_sample, scalar, nch, frame.sblimit);
 pick_scale (scalar, &frame, max_sc);

（4）选取心理声学模型，计算心理声学掩蔽阈

 if ((glopts.quickmode == TRUE) && (++psycount % glopts.quickcount != 0)) {
      /* We're using quick mode, so we're only calculating the model every
         'quickcount' frames. Otherwise, just copy the old ones across */
      for (ch = 0; ch < nch; ch++) {
	for (sb = 0; sb < SBLIMIT; sb++)
	  smr[ch][sb] = smrdef[ch][sb];
      }
    } else {
      /* calculate the psymodel */
      switch (model) {
      case -1:
	psycho_n1 (smr, nch);
	break;
      case 0:	/* Psy Model A */
	psycho_0 (smr, nch, scalar, (FLOAT) s_freq[header.version][header.sampling_frequency] * 1000);	
	break;
      case 1:
	psycho_1 (buffer, max_sc, smr, &frame);
	break;
      case 2:
	for (ch = 0; ch < nch; ch++) {
	  psycho_2 (&buffer[ch][0], &sam[ch][0], ch, &smr[ch][0], //snr32,
		     (FLOAT) s_freq[header.version][header.sampling_frequency] *
		     1000, &glopts);
	}
	break;
      case 3:
	/* Modified psy model 1 */
	psycho_3 (buffer, max_sc, smr, &frame, &glopts);
	break;
      case 4:
	/* Modified Psycho Model 2 */
	for (ch = 0; ch < nch; ch++) {
	  psycho_4 (&buffer[ch][0], &sam[ch][0], ch, &smr[ch][0], // snr32,
		     (FLOAT) s_freq[header.version][header.sampling_frequency] *
		     1000, &glopts);
	}
	break;	
      case 5:
	/* Model 5 comparse model 1 and 3 */
	psycho_1 (buffer, max_sc, smr, &frame);
	fprintf(stdout,"1 ");
	smr_dump(smr,nch);
	psycho_3 (buffer, max_sc, smr, &frame, &glopts);
	fprintf(stdout,"3 ");
	smr_dump(smr,nch);
	break;
      case 6:
	/* Model 6 compares model 2 and 4 */
	for (ch = 0; ch < nch; ch++) 
	  psycho_2 (&buffer[ch][0], &sam[ch][0], ch, &smr[ch][0], //snr32,
		    (FLOAT) s_freq[header.version][header.sampling_frequency] *
		    1000, &glopts);
	fprintf(stdout,"2 ");
	smr_dump(smr,nch);
	for (ch = 0; ch < nch; ch++) 
	  psycho_4 (&buffer[ch][0], &sam[ch][0], ch, &smr[ch][0], // snr32,
		     (FLOAT) s_freq[header.version][header.sampling_frequency] *
		     1000, &glopts);
	fprintf(stdout,"4 ");
	smr_dump(smr,nch);
	break;
      case 7:
	fprintf(stdout,"Frame: %i\n",frameNum);
	/* Dump the SMRs for all models */	
	psycho_1 (buffer, max_sc, smr, &frame);
	fprintf(stdout,"1");
	smr_dump(smr, nch);
	psycho_3 (buffer, max_sc, smr, &frame, &glopts);
	fprintf(stdout,"3");
	smr_dump(smr,nch);
	for (ch = 0; ch < nch; ch++) 
	  psycho_2 (&buffer[ch][0], &sam[ch][0], ch, &smr[ch][0], //snr32,
		    (FLOAT) s_freq[header.version][header.sampling_frequency] *
		    1000, &glopts);
	fprintf(stdout,"2");
	smr_dump(smr,nch);
	for (ch = 0; ch < nch; ch++) 
	  psycho_4 (&buffer[ch][0], &sam[ch][0], ch, &smr[ch][0], // snr32,
		     (FLOAT) s_freq[header.version][header.sampling_frequency] *
		     1000, &glopts);
	fprintf(stdout,"4");
	smr_dump(smr,nch);
	break;
      case 8:
	/* Compare 0 and 4 */	
	psycho_n1 (smr, nch);
	fprintf(stdout,"0");
	smr_dump(smr,nch);

	for (ch = 0; ch < nch; ch++) 
	  psycho_4 (&buffer[ch][0], &sam[ch][0], ch, &smr[ch][0], // snr32,
		     (FLOAT) s_freq[header.version][header.sampling_frequency] *
		     1000, &glopts);
	fprintf(stdout,"4");
	smr_dump(smr,nch);
	break;
      default:
	fprintf (stderr, "Invalid psy model specification: %i\n", model);
	exit (0);
      }

（5）根据模块4中得到的掩蔽值计算NMR，进行循环比特分配

transmission_pattern (scalar, scfsi, &frame);
main_bit_allocation (smr, scfsi, bit_alloc, &adb, &frame, &glopts);

（6）若开启了纠错模块，进行CRC检验

if (error_protection)
      CRC_calc (&frame, bit_alloc, scfsi, &crc);
    encode_info (&frame, &bs);
if (error_protection)
      encode_CRC (crc, &bs);   

（7）将经过比特分配信息、比例因子、比例因子选择信息打包至比特流

 transmission_pattern (scalar, scfsi, &frame);
 main_bit_allocation (smr, scfsi, bit_alloc, &adb, &frame, &glopts);
 、、、
 encode_bit_alloc (bit_alloc, &frame, &bs);
 encode_scale (bit_alloc, scfsi, scalar, &frame, &bs);

（8）对子带进行量化，装帧，送至比特流

subband_quantization (scalar, *sb_sample, j_scale, *j_sample, bit_alloc,
			  *subband, &frame);
sample_encoding (*subband, bit_alloc, &frame, &bs);
#endif


    /* If not all the bits were used, write out a stack of zeros */
    for (i = 0; i < adb; i++)
      put1bit (&bs, 0);
    if (header.dab_extension) {
      /* Reserve some bytes for X-PAD in DAB mode */
      putbits (&bs, 0, header.dab_length * 8);
      
      for (i = header.dab_extension - 1; i >= 0; i--) {
	CRC_calcDAB (&frame, bit_alloc, scfsi, scalar, &crc, i);
	/* this crc is for the previous frame in DAB mode  */
	if (bs.buf_byte_idx + lg_frame < bs.buf_size)
	  bs.buf[bs.buf_byte_idx + lg_frame] = crc;
	/* reserved 2 bytes for F-PAD in DAB mode  */
	putbits (&bs, crc, 8);
      }
      putbits (&bs, 0, 16);
    }

    frameBits = sstell (&bs) - sentBits;

    if (frameBits % 8) {	/* a program failure */
      fprintf (stderr, "Sent %ld bits = %ld slots plus %ld\n", frameBits,
	       frameBits / 8, frameBits % 8);
      fprintf (stderr, "If you are reading this, the program is broken\n");
      fprintf (stderr, "email [mfc at NOTplanckenerg.com] without the NOT\n");
      fprintf (stderr, "with the command line arguments and other info\n");
      exit (0);
    }

    sentBits += frameBits;
  }

  close_bit_stream_w (&bs);

2、输出音频的采样率和目标码率

在void parse_args（）函数中可以看到关于输入参数的代码：

while (++i < argc && err == 0) {
    char c, *token, *arg, *nextArg;
    int argUsed;

    token = argv[i];
    if (*token++ == '-') {
      if (i + 1 < argc)
	nextArg = argv[i + 1];
      else
	nextArg = "";
      argUsed = 0;
      if (!*token) {
	/* The user wants to use stdin and/or stdout. */
	if (inPath[0] == '\0')
	  strncpy (inPath, argv[i], MAX_NAME_SIZE);
	else if (outPath[0] == '\0')
	  strncpy (outPath, argv[i], MAX_NAME_SIZE);
      }

可知需要输入的命令参数为：输入文件名 输出文件名
运行代码，得到输出音频的采样率为16kHz，目标码率为96kbs

3、输出某个数据帧的比例因子、比特分配结果、该帧所分配的比特数

由编码原理和代码分析得出

• 每一帧所分配的比特数存放在变量adb中
• 子带数SBLIMIT设定为32，声道数为nch，每个声道都有SBLIMIT个子带；
• 比例因子存放在数组scalar中，每个子带有3个比例因子；
• 比特分配结果存放在bit_alloc数组中

在main函数的开头和末尾设定好文件流打开和关闭后，在帧循环中添加代码：

 while (get_audio (musicin, buffer, num_samples, nch, &header) > 0) {
    if (glopts.verbosity > 1)
、、、、
adb = available_bits (&header, &glopts);

	if (frameNum == 1)
	{
		fprintf(info,"当前帧:Frame[%d]\n",frameNum);
		fprintf(info,"分配的比特数：%d",adb);
	}
、、、、
transmission_pattern (scalar, scfsi, &frame);
main_bit_allocation (smr, scfsi, bit_alloc, &adb, &frame, &glopts);

	if (frameNum == 1)
	{
		fprintf(info, "\n比例因子:\n");
		for (int k = 0; k < nch; k++)
		{
			fprintf(info, "声道%d：\n", k + 1);
			for (int i = 0; i < SBLIMIT; i++)
			{
				fprintf(info, "子带[%d]：\t", i);
				for (int j = 0; j < 3; j++)
				{
					fprintf(info, "%d\t", scalar[k][j][i]);
				}
				fprintf(info, "\n");
			}


		}
		fprintf(info, "\n比特分配结果\n");
		for (int k = 0; k < nch; k++)
		{
			fprintf(info, "声道%d：\n", k + 1);
			for (int i = 0; i < SBLIMIT; i++)
			{
				fprintf(info, "子带[%d]：\t", i);
				fprintf(info, "%d\n", bit_alloc[k][i]);
			}
		}

	}

输出结果：