psych.c

AAC音频解码算法程序

/* CREATED BY :  Bernhard Grill -- August-96  */
#include <stdlib.h>
#include <math.h>
#include <memory.h>
#include "tf_main.h"
#include "psych.h"
#include "transfo.h"


SR_INFO sr_info_aac[MAX_SAMPLING_RATES+1] =
{
	{ 8000, 40, 15,
		{
			12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 16, 
			16, 16, 16, 16, 16, 16, 20, 20, 20, 20, 24, 24, 24, 28, 
			28, 32, 32, 36, 40, 44, 48, 52, 56, 60, 64, 80
		}, {
			4, 4, 4, 4, 4, 4, 4, 8, 8, 8, 8, 12, 16, 20, 20
		}
	}, { 11025, 43, 15,
		{
			8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 12, 12, 12, 
			12, 12, 12, 12, 12, 12, 16, 16, 16, 16, 20, 20, 20, 24, 
			24, 28, 28, 32, 36, 40, 40, 44, 48, 52, 56, 60, 64, 64, 64
		}, {
			4, 4, 4, 4, 4, 4, 4, 4, 8, 8, 12, 12, 16, 20, 20
		}
	}, { 12000, 43, 15,
		{
			8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 12, 12, 12, 
			12, 12, 12, 12, 12, 12, 16, 16, 16, 16, 20, 20, 20, 24, 
			24, 28, 28, 32, 36, 40, 40, 44, 48, 52, 56, 60, 64, 64, 64
		}, {
			4, 4, 4, 4, 4, 4, 4, 4, 8, 8, 12, 12, 16, 20, 20
		}
	}, { 16000, 43, 15,
		{
			8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 12, 12, 12, 
			12, 12, 12, 12, 12, 12, 16, 16, 16, 16, 20, 20, 20, 24, 
			24, 28, 28, 32, 36, 40, 40, 44, 48, 52, 56, 60, 64, 64, 64
		}, {
			4, 4, 4, 4, 4, 4, 4, 4, 8, 8, 12, 12, 16, 20, 20
		}
	}, { 22050, 47, 15,
		{
			4,  4,  4,  4,  4,  4,  4,  4,  4,  4,  4,  8,  8,  8,  8,  8,  8,  8,
			8,  8,  8,  12, 12, 12, 12, 16, 16, 16, 20, 20, 24, 24, 28, 28, 32,
			36, 36, 40, 44, 48, 52, 52, 64, 64, 64, 64, 64
		}, {
			4,  4,  4,  4,  4,  4,  4,  8,  8,  8,  8, 12, 16, 16, 20
		}
	},{ 24000, 47, 15,
		{
			4,  4,  4,  4,  4,  4,  4,  4,  4,  4,  4,  8,  8,  8,  8,  8,  8,  8,
			8,  8,  8,  12, 12, 12, 12, 16, 16, 16, 20, 20, 24, 24, 28, 28, 32,
			36, 36, 40, 44, 48, 52, 52, 64, 64, 64, 64, 64
		}, {
			4,  4,  4,  4,  4,  4,  4,  8,  8,  8,  8, 12, 16, 16, 20
		}
	}, { 32000, 51, 14,
		{
			4,	4,	4,	4,	4,	4,	4,	4,	4,	4,	8,	8,	8,	8,	
			8,	8,	8,	12,	12,	12,	12,	16,	16,	20,	20,	24,	24,	28,	
			28,	32,	32,	32,	32,	32,	32,	32,	32,	32,	32,	32,	32,	32,
			32,	32,	32,	32,	32,	32,	32,	32,	32
		},{
			4,	4,	4,	4,	4,	8,	8,	8,	12,	12,	12,	16,	16,	16
		}
	}, { 44100, 49, 14,
		{
			4,  4,  4,  4,  4,  4,  4,  4,  4,  4,  8,  8,  8,  8,  8,  8,  8, 
			12, 12, 12, 12, 16, 16, 20, 20, 24, 24, 28, 28, 32, 32, 32, 32, 32, 32,
			32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 96 
		}, {
			4,  4,  4,  4,  4,  8,  8,  8, 12, 12, 12, 16, 16, 16 
		}
	}, { 48000, 49, 14,
		{
			4,  4,  4,  4,  4,  4,  4,  4,  4,  4,  8,  8,  8,  8,  8,  8,  8, 
			12, 12, 12, 12, 16, 16, 20, 20, 24, 24, 28, 28, 32, 32, 32, 32, 32, 32,
			32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 96
		}, {
			4,  4,  4,  4,  4,  8,  8,  8, 12, 12, 12, 16, 16, 16 
		}
	}, {64000, 47, 12,
		{
			4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
			8, 8, 8, 8, 12, 12, 12, 16, 16, 16, 20, 24, 24, 28,
			36, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40,
			40, 40, 40, 40, 40
		},{
			4, 4, 4, 4, 4, 4, 8, 8, 8, 16, 28, 32
		}
	}, { 88200, 41, 12,
		{
			4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 
			8, 8, 8, 8, 8, 12, 12, 12, 12, 12, 16, 16, 24, 28, 
			36, 44, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64
		},{
			4, 4, 4, 4, 4, 4, 8, 8, 8, 16, 28, 36
		}
	}, { 96000, 41, 12,
		{
			4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 
			8, 8, 8, 8, 8, 12, 12, 12, 12, 12, 16, 16, 24, 28, 
			36, 44, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64, 64
		},{
			4, 4, 4, 4, 4, 4, 8, 8, 8, 16, 28, 36
		}
	},
	{ -1 }
};

double          sample[MAX_TIME_CHANNELS+2][BLOCK_LEN_LONG*2];
                               /* sample value */

FFT_TABLE_LONG    fft_tbl_long;  /* table for long fft */
FFT_TABLE_SHORT    fft_tbl_short;  /* table for short fft */
PARTITION_TABLE_LONG    part_tbl_long;  
PARTITION_TABLE_SHORT    part_tbl_short;
DYN_PART_TABLE_LONG     dyn_long;  
DYN_PART_TABLE_SHORT    dyn_short;
PSY_STATVARIABLE_LONG    psy_stvar_long[MAX_TIME_CHANNELS+2];
                               /* variables for long block */
PSY_STATVARIABLE_SHORT    psy_stvar_short[MAX_TIME_CHANNELS+2];
                               /* variables for short block */
/* added by T. Araki (1997.10.16) end */

void EncTf_psycho_acoustic_init( void )
{
	int chanNum;

	/* added by T. Araki (1997.10.16) */
	psy_fft_table_init(&fft_tbl_long, &fft_tbl_short);
	/* initializing fft table */
	for (chanNum=0;chanNum<MAX_TIME_CHANNELS+2;chanNum++) {
		psy_calc_init(&sample[chanNum], &psy_stvar_long[chanNum], &psy_stvar_short[chanNum]);
		/* initializing static variables */
	}
	/* added by T. Araki (1997.10.16) end */
}

/* added by T. Okada (1997.07.10) */
void psy_fft_table_init(FFT_TABLE_LONG *fft_tbl_long,
			FFT_TABLE_SHORT *fft_tbl_short
			)
{
	int i;

	/* generating Hann window */
	for(i = 0; i < BLOCK_LEN_LONG*2; i++)
		fft_tbl_long->hw[i] = 0.5 * (1-cos(2.0*M_PI*(i+0.5)/(BLOCK_LEN_LONG*2)));
	for(i = 0; i < BLOCK_LEN_SHORT*2; i++)
		fft_tbl_short->hw[i] = 0.5 * (1-cos(2.0*M_PI*(i+0.5)/(BLOCK_LEN_SHORT*2)));

	MakeFFTOrder();
}
/* added by T. Okada (1997.07.10) end */

/*
 * This Function calculates the Frequency in Hertz given a 
 * Bark-value. It uses the Traunmueller-formula for bark>2
 * and a linear inerpolation below. 
 * KAF
 */
double bark2hz (double bark)
{
    double hz;

    if(bark>2.0)
	hz = 1960 * (bark + 0.53) / (26.28 - bark);
    else
	hz = bark * 102.9;

    return (hz);
}

/*
 * This Function calculates the Frequency in Bark given a 
 * Frequency-value in Hertz. It uses the Traunmueller-formula 
 * for frequency>200Hz and a linear inerpolation below. 
 * KAF
 */
double hz2bark (double hz)
{
    double bark;

    if(hz>200.0)
	bark = 26.81 * hz / (1960 + hz) - 0.53; 
    else
	bark = hz / 102.9;

    return (bark);
}

/* added by T. Araki (1997.07.10) */
void psy_part_table_init(double sampling_rate,
						 PARTITION_TABLE_LONG *part_tbl_long,
						 PARTITION_TABLE_SHORT *part_tbl_short
						 )
{
    int b,bb; /* Jul 10 */
    double tmp;
	int partition[1024], j, w;
	int cbands, prev_cbound, crit_bands, cbound;

	cbands = (int)hz2bark(sampling_rate/2.0) + 1;
	cbands *= 3;
	part_tbl_long->sampling_rate = (int)sampling_rate;
	part_tbl_long->w_low[0] = 0;
	part_tbl_long->w_high[0] = 0;
	part_tbl_long->width[0] = 1;
    prev_cbound = 0;
	crit_bands = 0;
	for(j = 1; j <= cbands; j++)
	{
		cbound = (int)(bark2hz((double)j/3) * (double)BLOCK_LEN_LONG * 2.0 / sampling_rate + 0.5);
		if(cbound > prev_cbound) {
			crit_bands++;
			part_tbl_long->w_low[crit_bands] = min(prev_cbound,BLOCK_LEN_LONG-1);
			part_tbl_long->w_high[crit_bands] = min(cbound-1,BLOCK_LEN_LONG-1);
			part_tbl_long->width[crit_bands] = 
				part_tbl_long->w_high[crit_bands] - part_tbl_long->w_low[crit_bands] + 1;
			prev_cbound = cbound;
			if (part_tbl_long->w_high[crit_bands] == (BLOCK_LEN_LONG-1))
				break;
		}
	}
	part_tbl_long->len = crit_bands+1;
//	printf("%d %d\t",part_tbl_long->len, part_tbl_long->w_high[crit_bands]);

	part_tbl_short->sampling_rate = (int)sampling_rate;
	part_tbl_short->w_low[0] = 0;
	part_tbl_short->w_high[0] = 0;
	part_tbl_short->width[0] = 1;
    prev_cbound = 0;
	crit_bands = 0;
	for(j = 1; j <= cbands; j++)
	{
		cbound = (int)(bark2hz((double)j/3) * (double)BLOCK_LEN_SHORT * 2.0 / sampling_rate +0.5);
		if(cbound > prev_cbound) {
			crit_bands++;
			part_tbl_short->w_low[crit_bands] = min(prev_cbound,BLOCK_LEN_SHORT-1);
			part_tbl_short->w_high[crit_bands] = min(cbound-1,BLOCK_LEN_SHORT-1);
			part_tbl_short->width[crit_bands] = 
				part_tbl_short->w_high[crit_bands] - part_tbl_short->w_low[crit_bands] + 1;
			prev_cbound = cbound;
			if (part_tbl_short->w_high[crit_bands] == (BLOCK_LEN_SHORT-1))
				break;
		}
	}
	part_tbl_short->len = crit_bands+1;
//	printf("%d %d\n",part_tbl_short->len, part_tbl_short->w_high[crit_bands]);

	for (b = 0; b < part_tbl_long->len; b++) {
		for(w = part_tbl_long->w_low[b]; w <= part_tbl_long->w_high[b]; ++w){
			partition[w] = b;
		}
	}

	for(b = 0; b < part_tbl_long->len ; b++) {
		for (j=0;(b != partition[j]);j++);
		{
			double ji = j + (part_tbl_long->width[b]-1)/2.0;
			double freq = part_tbl_long->sampling_rate*ji/2048;
			double bark = 13*atan(0.00076*freq)+3.5*atan((freq/7500)*(freq/7500));
			dyn_long.bval[b] = bark;
		}
	}

	for (b = 0; b < part_tbl_short->len; b++) {
		for(w = part_tbl_short->w_low[b]; w <= part_tbl_short->w_high[b]; ++w){
			partition[w] = b;
		}
	}

	for(b = 0; b < part_tbl_short->len ; b++) {
		for (j=0;(b != partition[j]);j++);
		{
			double ji = j + (part_tbl_short->width[b]-1)/2.0;
			double freq = part_tbl_short->sampling_rate*ji/256;
			double bark = 13*atan(0.00076*freq) + 3.5*atan((freq/7500)*(freq/7500));
			dyn_short.bval[b]=bark;
		}
	}

	// Calculate the spreading function
	{
		double tmpx,tmpy,tmpz,b1,b2;
		int b, bb;

		for( b = 0; b < part_tbl_long->len; b++) {
			b2 = dyn_long.bval[b];
			for( bb = 0; bb < part_tbl_long->len; bb++) {
				b1 = dyn_long.bval[bb];

				//tmpx = (b2 >= b1) ? 3.0*(b2-b1) : 1.5*(b2-b1);
				tmpx = (b >= bb) ? 3.0*(b2-b1) : 1.5*(b2-b1);

				tmpz = 8.0 * psy_min( (tmpx-0.5)*(tmpx-0.5) - 2.0*(tmpx-0.5),0.0 );

				tmpy = 15.811389 + 7.5*(tmpx + 0.474)-17.5 *sqrt(1.0 + (tmpx+0.474)*(tmpx+0.474));

				dyn_long.spreading[b][bb] = ( tmpy < -100.0 ? 0.0 : pow(10.0, (tmpz + tmpy)/10.0) );
			}
		}

		for( b = 0; b < part_tbl_short->len; b++) {
			b2 = dyn_short.bval[b];
			for( bb = 0; bb < part_tbl_short->len; bb++) {
				b1 = dyn_short.bval[bb];

				//tmpx = (b2 >= b1) ? 3.0*(b2-b1) : 1.5*(b2-b1);
				tmpx = (b >= bb) ? 3.0*(b2-b1) : 1.5*(b2-b1);
				
				tmpz = 8.0 * psy_min( (tmpx-0.5)*(tmpx-0.5) - 2.0*(tmpx-0.5),0.0 );

				tmpy = 15.811389 + 7.5*(tmpx + 0.474)-17.5 *sqrt(1.0 + (tmpx+0.474)*(tmpx+0.474));

				dyn_short.spreading[b][bb] = ( tmpy < -100.0 ? 0.0 : pow(10.0, (tmpz + tmpy)/10.0) );
			}
		}
	}

    /* added by T. Okada (1997.07.10) */
    for( b = 0; b < part_tbl_long->len; b++){
		tmp = 0.0;
		for( bb = 0; bb < part_tbl_long->len; bb++)
			//tmp += sprdngf( (*part_tbl_long),(*part_tbl_short), bb, b, 0);
			tmp += dyn_long.spreading[bb][b];
		dyn_long.rnorm[b] = 1.0/tmp;
    }
    /* added by T. Okada (1997.07.10) end */

    /* added by T. Araki (1997.10.16) */
    for( b = 0; b < part_tbl_short->len; b++){
		tmp = 0.0;
		for( bb = 0; bb < part_tbl_short->len; bb++)
			//tmp += sprdngf( (*part_tbl_long), (*part_tbl_short), bb, b, 1);
			tmp += dyn_short.spreading[bb][b];
		dyn_short.rnorm[b] = 1.0/tmp;
    }
    /* added by T. Araki (1997.10.16) end */

	for(b = 0; b < part_tbl_long->len; b++) {
		dyn_long.bmax[b] = pow(10, -3*(0.5+0.5*(M_PI*(min(dyn_long.bval[b], 15.5)/15.5))));
	}
	for(b = 0; b < part_tbl_short->len; b++) {
		dyn_short.bmax[b] = pow(10, -3*(0.5+0.5*(M_PI*(min(dyn_short.bval[b], 15.5)/15.5))));
	}

	part_tbl_long->dyn = &dyn_long;
	part_tbl_short->dyn = &dyn_short;
}


void psy_calc_init( 
		   double sample[][BLOCK_LEN_LONG*2],
		   PSY_STATVARIABLE_LONG *psy_stvar_long,
		   PSY_STATVARIABLE_SHORT *psy_stvar_short
		   )
{
	int ch = 0;
	int i;

	for(i = 0; i < BLOCK_LEN_LONG*2; i++){
		sample[ch][i] = 0.0;
	}

    for(i = 0; i < BLOCK_LEN_LONG*3; i++){
		psy_stvar_long->fft_r[i] = 0.0;
		psy_stvar_long->fft_f[i] = 0.0;
    }

	psy_stvar_long->p_fft = 0;

    for(i = 0; i < NPART_LONG*2; i++){
		psy_stvar_long->nb[i] = 90.0;
    }

	psy_stvar_long->p_nb = NPART_LONG;
	/* added by T. Araki (1997.07.10) end */

	/* added by T. Araki (1997.10.16) */
	for(i = 0; i < BLOCK_LEN_SHORT; i++) {
		psy_stvar_short->last6_fft_r[i] = 0.0;
		psy_stvar_short->last6_fft_f[i] = 0.0;
		psy_stvar_short->last7_fft_r[i] = 0.0;
		psy_stvar_short->last7_fft_f[i] = 0.0;
    }

    for(i = 0; i < NPART_SHORT; i++){
		psy_stvar_short->last7_nb[i] = 90.0;
    }
	/* added by T. Araki (1997.10.16) end */
}

void psy_fill_lookahead(double *p_time_signal[], int no_of_chan)
{
	int i, ch;

	for (ch = 0; ch < no_of_chan; ch++) {
		for(i = 0; i < BLOCK_LEN_LONG; i++){
			sample[ch][i+BLOCK_LEN_LONG] = p_time_signal[ch][i];
		}
	}
}

/* main */
void EncTf_psycho_acoustic(