aacquant.c

《Visual C++小波变换技术与工程实践》作者：靳济芳。书上的代码。第3章：语音的去噪处理

/*
 * FAAC - Freeware Advanced Audio Coder
 * Copyright (C) 2001 Menno Bakker
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.

 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 * $Id: aacquant.c,v 1.3 2001/06/04 23:02:24 wmay Exp $
 */

#include <math.h>
#include <stdlib.h>

#include "aacquant.h"
#include "coder.h"
#include "huffman.h"
#include "psych.h"
#include "util.h"


#define XRPOW_FTOI(src,dest) ((dest) = (int)(src))
#define QUANTFAC(rx)  adj43[rx]
#define ROUNDFAC 0.4054


double *pow43;
double *adj43;
double *adj43asm;


void AACQuantizeInit(CoderInfo *coderInfo, unsigned int numChannels)
{
	unsigned int channel, i;

	pow43 = (double*)AllocMemory(PRECALC_SIZE*sizeof(double));
	adj43 = (double*)AllocMemory(PRECALC_SIZE*sizeof(double));
	adj43asm = (double*)AllocMemory(PRECALC_SIZE*sizeof(double));

	pow43[0] = 0.0;
	for(i=1;i<PRECALC_SIZE;i++)
		pow43[i] = pow((double)i, 4.0/3.0);

	adj43asm[0] = 0.0;
	for (i = 1; i < PRECALC_SIZE; i++)
		adj43asm[i] = i - 0.5 - pow(0.5 * (pow43[i - 1] + pow43[i]),0.75);
	for (i = 0; i < PRECALC_SIZE-1; i++)
		adj43[i] = (i + 1) - pow(0.5 * (pow43[i] + pow43[i + 1]), 0.75);
	adj43[i] = 0.5;

	for (channel = 0; channel < numChannels; channel++) {
		coderInfo[channel].old_value = 0;
		coderInfo[channel].CurrentStep = 4;

		coderInfo[channel].requantFreq = (double*)AllocMemory(BLOCK_LEN_LONG*sizeof(double));
	}
}

void AACQuantizeEnd(CoderInfo *coderInfo, unsigned int numChannels)
{
	unsigned int channel;

	if (pow43) FreeMemory(pow43);
	if (adj43) FreeMemory(adj43);
	if (adj43asm) FreeMemory(adj43asm);

	for (channel = 0; channel < numChannels; channel++) {
		if (coderInfo[channel].requantFreq) FreeMemory(coderInfo[channel].requantFreq);
	}
}

int AACQuantize(CoderInfo *coderInfo,
				PsyInfo *psyInfo,
				ChannelInfo *channelInfo,
				int *cb_width,
				int num_cb,
				double *xr,
				int desired_rate)
{
	int sb, i, do_q = 0;
	int bits, sign;
	double *xr_pow, *xmin;
	int *xi;

	/* Use local copy's */
	int *scale_factor = coderInfo->scale_factor;


	xr_pow = (double*)AllocMemory(FRAME_LEN*sizeof(double));
	xmin = (double*)AllocMemory(MAX_SCFAC_BANDS*sizeof(double));
	xi = (int*)AllocMemory(FRAME_LEN*sizeof(int));


	if (coderInfo->block_type == ONLY_SHORT_WINDOW) {
		SortForGrouping(coderInfo, psyInfo, channelInfo, cb_width, xr);
	} else {
		for (sb = 0; sb < coderInfo->nr_of_sfb; sb++) {
			if (channelInfo->msInfo.is_present && channelInfo->msInfo.ms_used[sb]) {
				psyInfo->maskThr[sb] = psyInfo->maskThrMS[sb];
				psyInfo->maskEn[sb] = psyInfo->maskEnMS[sb];
			}
		}
	}


	/* Set all scalefactors to 0 */
	coderInfo->global_gain = 0;
	for (sb = 0; sb < coderInfo->nr_of_sfb; sb++)
		scale_factor[sb] = 0;

	/* Compute xr_pow */
	for (i = 0; i < FRAME_LEN; i++) {
        double temp = fabs(xr[i]);
        xr_pow[i] = sqrt(temp * sqrt(temp));
        do_q += (temp > 1E-20);
    }

	if (do_q) {
		bits = SearchStepSize(coderInfo, desired_rate, xr_pow, xi);
		coderInfo->old_value = coderInfo->global_gain;

		CalcAllowedDist(psyInfo, cb_width, num_cb, xr, xmin);
		OuterLoop(coderInfo, xr, xr_pow, xi, xmin, desired_rate);

		for ( i = 0; i < FRAME_LEN; i++ )  {
			sign = (xr[i] < 0) ? -1 : 1;
			xi[i] *= sign;
		}
	} else {
		coderInfo->global_gain = 0;
		SetMemory(xi, 0, FRAME_LEN*sizeof(int));
	}

	CountBitsLong(coderInfo, xi);

	/* offset the difference of common_scalefac and scalefactors by SF_OFFSET  */
	for (i = 0; i < coderInfo->nr_of_sfb; i++) {
		if ((coderInfo->book_vector[i]!=INTENSITY_HCB)&&(coderInfo->book_vector[i]!=INTENSITY_HCB2)) {
			scale_factor[i] = coderInfo->global_gain - scale_factor[i] + SF_OFFSET;
		}
	}
	coderInfo->global_gain = scale_factor[0];

	/* place the codewords and their respective lengths in arrays data[] and len[] respectively */
	/* there are 'counter' elements in each array, and these are variable length arrays depending on the input */
	coderInfo->spectral_count = 0;
	for(i = 0; i < coderInfo->nr_of_sfb; i++) {
		OutputBits(
			coderInfo,
			coderInfo->book_vector[i],
			xi,
			coderInfo->sfb_offset[i],
			coderInfo->sfb_offset[i+1]-coderInfo->sfb_offset[i]);
	}


	if (xmin) FreeMemory(xmin);
	if (xr_pow) FreeMemory(xr_pow);
	if (xi) FreeMemory(xi);

	return bits;
}

static int SearchStepSize(CoderInfo *coderInfo,
						  const int desired_rate,
						  const double *xr,
						  int *xi)
{
	int flag_GoneOver = 0;
	int CurrentStep = coderInfo->CurrentStep;
	int nBits;
	int StepSize = coderInfo->old_value;
	int Direction = 0;

	do
	{
		coderInfo->global_gain = StepSize;
		nBits = CountBits(coderInfo, xi, xr);

		if (CurrentStep == 1 ) {
			break; /* nothing to adjust anymore */
		}
		if (flag_GoneOver) {
			CurrentStep /= 2;
		}
		if (nBits > desired_rate) { /* increase Quantize_StepSize */
			if (Direction == -1 && !flag_GoneOver) {
				flag_GoneOver = 1;
				CurrentStep /= 2; /* late adjust */
			}
			Direction = 1;
			StepSize += CurrentStep;
		} else if (nBits < desired_rate) {
			if (Direction == 1 && !flag_GoneOver) {
				flag_GoneOver = 1;
				CurrentStep /= 2; /* late adjust */
			}
			Direction = -1;
			StepSize -= CurrentStep;
		} else break;
	} while (1);

	CurrentStep = coderInfo->old_value - StepSize;

	coderInfo->CurrentStep = CurrentStep/4 != 0 ? 4 : 2;
	coderInfo->old_value = coderInfo->global_gain;

	return nBits;
}

#if 1 /* TAKEHIRO_IEEE754_HACK */

#pragma warning( disable : 4244 4307 )

typedef union {
    float f;
    int i;
} fi_union;

#define MAGIC_FLOAT (65536*(128))
#define MAGIC_INT 0x4b000000

static void Quantize(const double *xp, int *pi, double istep)
{
	int j;
	fi_union *fi;

	fi = (fi_union *)pi;
	for (j = FRAME_LEN/4 - 1; j >= 0; --j) {
		double x0 = istep * xp[0];
		double x1 = istep * xp[1];
		double x2 = istep * xp[2];
		double x3 = istep * xp[3];

		x0 += MAGIC_FLOAT; fi[0].f = x0;
		x1 += MAGIC_FLOAT; fi[1].f = x1;
		x2 += MAGIC_FLOAT; fi[2].f = x2;
		x3 += MAGIC_FLOAT; fi[3].f = x3;

		fi[0].f = x0 + (adj43asm - MAGIC_INT)[fi[0].i];
		fi[1].f = x1 + (adj43asm - MAGIC_INT)[fi[1].i];
		fi[2].f = x2 + (adj43asm - MAGIC_INT)[fi[2].i];
		fi[3].f = x3 + (adj43asm - MAGIC_INT)[fi[3].i];

		fi[0].i -= MAGIC_INT;
		fi[1].i -= MAGIC_INT;
		fi[2].i -= MAGIC_INT;
		fi[3].i -= MAGIC_INT;
		fi += 4;
		xp += 4;
	}
}
#else
static void Quantize(const double *xr, int *ix, double istep)
{
	int j;

    for (j = FRAME_LEN/8; j > 0; --j) {
		double x1, x2, x3, x4, x5, x6, x7, x8;
		int rx1, rx2, rx3, rx4, rx5, rx6, rx7, rx8;

		x1 = *xr++ * istep;
		x2 = *xr++ * istep;
		XRPOW_FTOI(x1, rx1);
		x3 = *xr++ * istep;
		XRPOW_FTOI(x2, rx2);
		x4 = *xr++ * istep;
		XRPOW_FTOI(x3, rx3);
		x5 = *xr++ * istep;
		XRPOW_FTOI(x4, rx4);
		x6 = *xr++ * istep;
		XRPOW_FTOI(x5, rx5);
		x7 = *xr++ * istep;
		XRPOW_FTOI(x6, rx6);
		x8 = *xr++ * istep;
		XRPOW_FTOI(x7, rx7);
		x1 += QUANTFAC(rx1);
		XRPOW_FTOI(x8, rx8);
		x2 += QUANTFAC(rx2);
		XRPOW_FTOI(x1,*ix++);
		x3 += QUANTFAC(rx3);
		XRPOW_FTOI(x2,*ix++);
		x4 += QUANTFAC(rx4);
		XRPOW_FTOI(x3,*ix++);
		x5 += QUANTFAC(rx5);
		XRPOW_FTOI(x4,*ix++);
		x6 += QUANTFAC(rx6);
		XRPOW_FTOI(x5,*ix++);
		x7 += QUANTFAC(rx7);
		XRPOW_FTOI(x6,*ix++);
		x8 += QUANTFAC(rx8);
		XRPOW_FTOI(x7,*ix++);
		XRPOW_FTOI(x8,*ix++);
    }
}
#endif

static int CountBitsLong(CoderInfo *coderInfo, int *xi)
{
	int i, bits = 0;

	/* find a good method to section the scalefactor bands into huffman codebook sections */
	BitSearch(coderInfo, xi);

	/* calculate the amount of bits needed for encoding the huffman codebook numbers */
	bits += SortBookNumbers(coderInfo, NULL, 0);

	/* calculate the amount of bits needed for the spectral values */
	coderInfo->spectral_count = 0;
	for(i = 0; i < coderInfo->nr_of_sfb; i++) {  
		bits += CalcBits(coderInfo,
			coderInfo->book_vector[i],
			xi,
			coderInfo->sfb_offset[i], 
			coderInfo->sfb_offset[i+1] - coderInfo->sfb_offset[i]);
	}

	/* the number of bits for the scalefactors */
	bits += WriteScalefactors(coderInfo, NULL, 0);

	/* the total amount of bits required */
	return bits;
}

static int CountBits(CoderInfo *coderInfo, int *ix, const double *xr)
{
	int bits = 0, i;

	/* since quantize uses table lookup, we need to check this first: */
	double w = (IXMAX_VAL) / IPOW20(coderInfo->global_gain);
	for ( i = 0; i < FRAME_LEN; i++ )  {
		if (xr[i] > w)
			return LARGE_BITS;
	}

	Quantize(xr, ix, IPOW20(coderInfo->global_gain));

	bits = CountBitsLong(coderInfo, ix);

	return bits;
}