📄 psymodel.c
字号:
/**********************************************************************
* date programmers comment *
* 2/25/91 Davis Pan start of version 1.0 records *
* 5/10/91 W. Joseph Carter Ported to Macintosh and Unix. *
* 7/10/91 Earle Jennings Ported to MsDos. *
* replace of floats with FLOAT *
* 2/11/92 W. Joseph Carter Fixed mem_alloc() arg for "absthr". *
* 3/16/92 Masahiro Iwadare Modification for Layer III *
* 17/4/93 Masahiro Iwadare Updated for IS Modification *
**********************************************************************/
#include "util.h"
#include "encoder.h"
#include "psymodel.h"
#include "l3side.h"
#include <assert.h>
#ifdef HAVEGTK
#include "gtkanal.h"
#endif
#include "tables.h"
#include "fft.h"
#ifdef M_LN10
#define LN_TO_LOG10 (M_LN10/10)
#else
#define LN_TO_LOG10 0.2302585093
#endif
void L3para_read( FLOAT8 sfreq, int numlines[CBANDS],int numlines_s[CBANDS], int partition_l[HBLKSIZE],
FLOAT8 minval[CBANDS], FLOAT8 qthr_l[CBANDS],
FLOAT8 s3_l[CBANDS + 1][CBANDS + 1],
FLOAT8 s3_s[CBANDS + 1][CBANDS + 1],
FLOAT8 qthr_s[CBANDS],
FLOAT8 SNR_s[CBANDS],
int bu_l[SBPSY_l], int bo_l[SBPSY_l],
FLOAT8 w1_l[SBPSY_l], FLOAT8 w2_l[SBPSY_l],
int bu_s[SBPSY_s], int bo_s[SBPSY_s],
FLOAT8 w1_s[SBPSY_s], FLOAT8 w2_s[SBPSY_s] );
void L3psycho_anal( lame_global_flags *gfp,
short int *buffer[2],int gr_out ,
FLOAT8 *ms_ratio,
FLOAT8 *ms_ratio_next,
FLOAT8 *ms_ener_ratio,
III_psy_ratio masking_ratio[2][2],
III_psy_ratio masking_MS_ratio[2][2],
FLOAT8 percep_entropy[2],FLOAT8 percep_MS_entropy[2],
int blocktype_d[2])
{
/* to get a good cache performance, one has to think about
* the sequence, in which the variables are used
*/
/* The static variables "r", "phi_sav", "new", "old" and "oldest" have */
/* to be remembered for the unpredictability measure. For "r" and */
/* "phi_sav", the first index from the left is the channel select and */
/* the second index is the "age" of the data. */
static FLOAT8 minval[CBANDS],qthr_l[CBANDS];
static FLOAT8 qthr_s[CBANDS];
static FLOAT8 nb_1[4][CBANDS], nb_2[4][CBANDS];
static FLOAT8 s3_s[CBANDS + 1][CBANDS + 1];
static FLOAT8 s3_l[CBANDS + 1][CBANDS + 1];
static III_psy_xmin thm[4];
static III_psy_xmin en[4];
/* unpredictability calculation
*/
static int cw_upper_index;
static int cw_lower_index;
static FLOAT ax_sav[4][2][HBLKSIZE];
static FLOAT bx_sav[4][2][HBLKSIZE];
static FLOAT rx_sav[4][2][HBLKSIZE];
static FLOAT cw[HBLKSIZE];
/* fft and energy calculation
*/
FLOAT (*wsamp_l)[BLKSIZE];
FLOAT (*wsamp_s)[3][BLKSIZE_s];
FLOAT tot_ener[4];
static FLOAT wsamp_L[2][BLKSIZE];
static FLOAT energy[HBLKSIZE];
static FLOAT wsamp_S[2][3][BLKSIZE_s];
static FLOAT energy_s[3][HBLKSIZE_s];
/* convolution
*/
static FLOAT8 eb[CBANDS];
static FLOAT8 cb[CBANDS];
static FLOAT8 thr[CBANDS];
/* Scale Factor Bands
*/
static FLOAT8 w1_l[SBPSY_l], w2_l[SBPSY_l];
static FLOAT8 w1_s[SBPSY_s], w2_s[SBPSY_s];
static FLOAT8 mld_l[SBPSY_l],mld_s[SBPSY_s];
static int bu_l[SBPSY_l],bo_l[SBPSY_l] ;
static int bu_s[SBPSY_s],bo_s[SBPSY_s] ;
static int npart_l,npart_s;
static int npart_l_orig,npart_s_orig;
static int s3ind[CBANDS][2];
static int s3ind_s[CBANDS][2];
static int numlines_s[CBANDS] ;
static int numlines_l[CBANDS];
static int partition_l[HBLKSIZE];
/* frame analyzer
*/
#ifdef HAVEGTK
static FLOAT energy_save[4][HBLKSIZE];
static FLOAT8 pe_save[4];
static FLOAT8 ers_save[4];
#endif
/* ratios
*/
static FLOAT8 pe[4]={0,0,0,0};
static FLOAT8 ms_ratio_s_old=0,ms_ratio_l_old=0;
static FLOAT8 ms_ener_ratio_old=.25;
FLOAT8 ms_ratio_l=0,ms_ratio_s=0;
/* block type
*/
static int blocktype_old[2];
int blocktype[2],uselongblock[2];
/* usual variables like loop indices, etc..
*/
int numchn, chn;
int b, i, j, k;
int sb,sblock;
FLOAT cwlimit;
/* initialization of static variables
*/
if((gfp->frameNum==0) && (gr_out==0)){
FLOAT8 SNR_s[CBANDS];
blocktype_old[0]=STOP_TYPE;
blocktype_old[1]=STOP_TYPE;
i = gfp->out_samplerate;
switch(i){
case 32000: break;
case 44100: break;
case 48000: break;
case 16000: break;
case 22050: break;
case 24000: break;
default: fprintf(stderr,"error, invalid sampling frequency: %d Hz\n",i);
exit(-1);
}
/* reset states used in unpredictability measure */
memset (rx_sav,0, sizeof(rx_sav));
memset (ax_sav,0, sizeof(ax_sav));
memset (bx_sav,0, sizeof(bx_sav));
memset (en,0, sizeof(en));
memset (thm,0, sizeof(thm));
/* gfp->cwlimit = sfreq*j/1024.0; */
cw_lower_index=6;
if (gfp->cwlimit>0)
cwlimit=gfp->cwlimit;
else
cwlimit=8.8717;
cw_upper_index = cwlimit*1000.0*1024.0/((FLOAT8) gfp->out_samplerate);
cw_upper_index=Min(HBLKSIZE-4,cw_upper_index); /* j+3 < HBLKSIZE-1 */
cw_upper_index=Max(6,cw_upper_index);
for ( j = 0; j < HBLKSIZE; j++ )
cw[j] = 0.4;
/* setup stereo demasking thresholds */
/* formula reverse enginerred from plot in paper */
for ( sb = 0; sb < SBPSY_s; sb++ ) {
FLOAT8 mld = 1.25*(1-cos(PI*sb/SBPSY_s))-2.5;
mld_s[sb] = pow(10.0,mld);
}
for ( sb = 0; sb < SBPSY_l; sb++ ) {
FLOAT8 mld = 1.25*(1-cos(PI*sb/SBPSY_l))-2.5;
mld_l[sb] = pow(10.0,mld);
}
for (i=0;i<HBLKSIZE;i++) partition_l[i]=-1;
L3para_read( (FLOAT8) gfp->out_samplerate,numlines_l,numlines_s,partition_l,minval,qthr_l,s3_l,s3_s,
qthr_s,SNR_s,
bu_l,bo_l,w1_l,w2_l, bu_s,bo_s,w1_s,w2_s );
/* npart_l_orig = number of partition bands before convolution */
/* npart_l = number of partition bands after convolution */
npart_l_orig=0; npart_s_orig=0;
for (i=0;i<HBLKSIZE;i++)
if (partition_l[i]>npart_l_orig) npart_l_orig=partition_l[i];
npart_l_orig++;
for (i=0;numlines_s[i]>=0;i++)
;
npart_s_orig = i;
npart_l=bo_l[SBPSY_l-1]+1;
npart_s=bo_s[SBPSY_s-1]+1;
/* MPEG2 tables are screwed up
* the mapping from paritition bands to scalefactor bands will use
* more paritition bands than we have.
* So we will not compute these fictitious partition bands by reducing
* npart_l below. */
if (npart_l > npart_l_orig) {
npart_l=npart_l_orig;
bo_l[SBPSY_l-1]=npart_l-1;
w2_l[SBPSY_l-1]=1.0;
}
if (npart_s > npart_s_orig) {
npart_s=npart_s_orig;
bo_s[SBPSY_s-1]=npart_s-1;
w2_s[SBPSY_s-1]=1.0;
}
for (i=0; i<npart_l; i++) {
for (j = 0; j < npart_l_orig; j++) {
if (s3_l[i][j] != 0.0)
break;
}
s3ind[i][0] = j;
for (j = npart_l_orig - 1; j > 0; j--) {
if (s3_l[i][j] != 0.0)
break;
}
s3ind[i][1] = j;
}
for (i=0; i<npart_s; i++) {
for (j = 0; j < npart_s_orig; j++) {
if (s3_s[i][j] != 0.0)
break;
}
s3ind_s[i][0] = j;
for (j = npart_s_orig - 1; j > 0; j--) {
if (s3_s[i][j] != 0.0)
break;
}
s3ind_s[i][1] = j;
}
/*
#include "debugscalefac.c"
*/
#define AACS3
#define NEWS3XX
#ifdef AACS3
/* AAC values, results in more masking over MP3 values */
# define TMN 18
# define NMT 6
#else
/* MP3 values */
# define TMN 29
# define NMT 6
#endif
#define rpelev 2
#define rpelev2 16
/* compute norm_l, norm_s instead of relying on table data */
for ( b = 0;b < npart_l; b++ ) {
FLOAT8 norm=0;
for ( k = s3ind[b][0]; k <= s3ind[b][1]; k++ ) {
norm += s3_l[b][k];
}
for ( k = s3ind[b][0]; k <= s3ind[b][1]; k++ ) {
s3_l[b][k] *= exp(-LN_TO_LOG10 * NMT) / norm;
}
/*printf("%i norm=%f norm_l=%f \n",b,1/norm,norm_l[b]);*/
}
/* MPEG1 SNR_s data is given in db, convert to energy */
if (gfp->version == 1) {
for ( b = 0;b < npart_s; b++ ) {
SNR_s[b]=exp( (FLOAT8) SNR_s[b] * LN_TO_LOG10 );
}
}
for ( b = 0;b < npart_s; b++ ) {
FLOAT8 norm=0;
for ( k = s3ind_s[b][0]; k <= s3ind_s[b][1]; k++ ) {
norm += s3_s[b][k];
}
for ( k = s3ind_s[b][0]; k <= s3ind_s[b][1]; k++ ) {
s3_s[b][k] *= SNR_s[b] / norm;
}
/*printf("%i norm=%f norm_s=%f \n",b,1/norm,norm_l[b]);*/
}
init_fft();
}
/************************* End of Initialization *****************************/
numchn = gfp->stereo;
/* chn=2 and 3 = Mid and Side channels */
if (gfp->mode == MPG_MD_JOINT_STEREO) numchn=4;
for (chn=0; chn<numchn; chn++) {
wsamp_s = wsamp_S+(chn & 1);
wsamp_l = wsamp_L+(chn & 1);
if (chn<2) {
/**********************************************************************
* compute FFTs
**********************************************************************/
fft_long ( *wsamp_l, chn, buffer);
fft_short( *wsamp_s, chn, buffer);
/* LR maskings */
percep_entropy[chn] = pe[chn];
masking_ratio[gr_out][chn].thm = thm[chn];
masking_ratio[gr_out][chn].en = en[chn];
}else{
/* MS maskings */
percep_MS_entropy[chn-2] = pe[chn];
masking_MS_ratio[gr_out][chn-2].en = en[chn];
masking_MS_ratio[gr_out][chn-2].thm = thm[chn];
if (chn == 2)
{
for (j = BLKSIZE-1; j >=0 ; --j)
{
FLOAT l = wsamp_L[0][j];
FLOAT r = wsamp_L[1][j];
wsamp_L[0][j] = (l+r)*(FLOAT)(SQRT2*0.5);
wsamp_L[1][j] = (l-r)*(FLOAT)(SQRT2*0.5);
}
for (b = 2; b >= 0; --b)
{
for (j = BLKSIZE_s-1; j >= 0 ; --j)
{
FLOAT l = wsamp_S[0][b][j];
FLOAT r = wsamp_S[1][b][j];
wsamp_S[0][b][j] = (l+r)*(FLOAT)(SQRT2*0.5);
wsamp_S[1][b][j] = (l-r)*(FLOAT)(SQRT2*0.5);
}
}
}
}
/**********************************************************************
* compute energies
**********************************************************************/
energy[0] = (*wsamp_l)[0];
energy[0] *= energy[0];
tot_ener[chn] = energy[0]; /* sum total energy at nearly no extra cost */
for (j=BLKSIZE/2-1; j >= 0; --j)
{
FLOAT re = (*wsamp_l)[BLKSIZE/2-j];
FLOAT im = (*wsamp_l)[BLKSIZE/2+j];
energy[BLKSIZE/2-j] = (re * re + im * im) * (FLOAT)0.5;
tot_ener[chn] += energy[BLKSIZE/2-j];
}
for (b = 2; b >= 0; --b)
{
energy_s[b][0] = (*wsamp_s)[b][0];
energy_s[b][0] *= energy_s [b][0];
for (j=BLKSIZE_s/2-1; j >= 0; --j)
{
FLOAT re = (*wsamp_s)[b][BLKSIZE_s/2-j];
FLOAT im = (*wsamp_s)[b][BLKSIZE_s/2+j];
energy_s[b][BLKSIZE_s/2-j] = (re * re + im * im) * (FLOAT)0.5;
}
}
#ifdef HAVEGTK
if(gfp->gtkflag) {
for (j=0; j<HBLKSIZE ; j++) {
pinfo->energy[gr_out][chn][j]=energy_save[chn][j];
energy_save[chn][j]=energy[j];
}
}
#endif
/**********************************************************************
* compute unpredicatability of first six spectral lines *
**********************************************************************/
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -