📄 l3psy.c
字号:
/**********************************************************************
* ISO MPEG Audio Subgroup Software Simulation Group (1996)
* ISO 13818-3 MPEG-2 Audio Encoder - Lower Sampling Frequency Extension
*
* $Id: l3psy.c,v 1.2 1997/01/19 22:28:29 rowlands Exp $
*
* $Log: l3psy.c,v $
* Revision 1.2 1997/01/19 22:28:29 rowlands
* Layer 3 bug fixes from Seymour Shlien
*
* Revision 1.1 1996/02/14 04:04:23 rowlands
* Initial revision
*
* Received from Mike Coleman
**********************************************************************/
/**********************************************************************
* 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 "common.h"
#include "encoder.h"
#include "l3psy.h"
#include "l3side.h"
#include <assert.h>
#define maximum(x,y) ( (x>y) ? x : y )
#define minimum(x,y) ( (x<y) ? x : y )
void L3para_read( double sfreq, int numlines[CBANDS], int partition_l[HBLKSIZE],
double minval[CBANDS], double qthr_l[CBANDS], double norm_l[CBANDS],
double s3_l[CBANDS][CBANDS], int partition_s[HBLKSIZE_s], double qthr_s[CBANDS_s],
double norm_s[CBANDS_s], double SNR_s[CBANDS_s],
int cbw_l[SBMAX_l], int bu_l[SBMAX_l], int bo_l[SBMAX_l],
double w1_l[SBMAX_l], double w2_l[SBMAX_l],
int cbw_s[SBMAX_s], int bu_s[SBMAX_s], int bo_s[SBMAX_s],
double w1_s[SBMAX_s], double w2_s[SBMAX_s] );
void L3psycho_anal( short int *buffer, short int savebuf[1344], int chn, int lay, FLOAT snr32[32],
double sfreq, double ratio_d[21], double ratio_ds[12][3],
double *pe, gr_info *cod_info )
{
static double ratio[2][21];
static double ratio_s[2][12][3];
int blocktype;
unsigned int b, i, j, k;
double r_prime, phi_prime; /* not FLOAT */
FLOAT freq_mult, bval_lo, min_thres, sum_energy;
double tb, temp1,temp2,temp3;
/* nint(); Layer III */
double thr[CBANDS];
/* 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 FLOAT window_s[BLKSIZE_s] ;
static int new = 0, old = 1, oldest = 0;
static int init = 0, flush, sync_flush, syncsize, sfreq_idx;
static double cw[HBLKSIZE], eb[CBANDS];
static double ctb[CBANDS];
static double SNR_l[CBANDS], SNR_s[CBANDS_s];
static int init_L3;
static double minval[CBANDS],qthr_l[CBANDS],norm_l[CBANDS];
static double qthr_s[CBANDS_s],norm_s[CBANDS_s];
static double nb_1[2][CBANDS], nb_2[2][CBANDS];
static double s3_l[CBANDS][CBANDS]; /* s3_s[CBANDS_s][CBANDS_s]; */
/* Scale Factor Bands */
static int cbw_l[SBMAX_l],bu_l[SBMAX_l],bo_l[SBMAX_l] ;
static int cbw_s[SBMAX_s],bu_s[SBMAX_s],bo_s[SBMAX_s] ;
static double w1_l[SBMAX_l], w2_l[SBMAX_l];
static double w1_s[SBMAX_s], w2_s[SBMAX_s];
static double en[SBMAX_l], thm[SBMAX_l] ;
static int blocktype_old[2] ;
int sb,sblock;
static int partition_l[HBLKSIZE],partition_s[HBLKSIZE_s];
/* The following static variables are constants. */
static double nmt = 5.5;
static FLOAT crit_band[27] = {0, 100, 200, 300, 400, 510, 630, 770,
920, 1080, 1270,1480,1720,2000,2320, 2700,
3150, 3700, 4400,5300,6400,7700,9500,12000,
15500,25000,30000};
static FLOAT bmax[27] = {20.0, 20.0, 20.0, 20.0, 20.0, 17.0, 15.0,
10.0, 7.0, 4.4, 4.5, 4.5, 4.5, 4.5,
4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5,
4.5, 4.5, 4.5, 3.5, 3.5, 3.5};
/* The following pointer variables point to large areas of memory */
/* dynamically allocated by the mem_alloc() function. Dynamic memory */
/* allocation is used in order to avoid stack frame or data area */
/* overflow errors that otherwise would have occurred at compile time */
/* on the Macintosh computer. */
FLOAT *grouped_c, *grouped_e, *nb, *cb, *ecb, *bc;
FLOAT *wsamp_r, *wsamp_i, *phi, *energy;
static FLOAT energy_s[3][256];
static FLOAT phi_s[3][256] ; /* 256 samples not 129 */
FLOAT *c, *fthr;
F32 *snrtmp;
static int *numlines ;
static int *partition;
static FLOAT *cbval, *rnorm;
static FLOAT *window;
static FLOAT *absthr;
static double *tmn;
static FCB *s;
static FHBLK *lthr;
static F2HBLK *r, *phi_sav;
/* These dynamic memory allocations simulate "automatic" variables */
/* placed on the stack. For each mem_alloc() call here, there must be */
/* a corresponding mem_free() call at the end of this function. */
grouped_c = (FLOAT *) mem_alloc(sizeof(FCB), "grouped_c");
grouped_e = (FLOAT *) mem_alloc(sizeof(FCB), "grouped_e");
nb = (FLOAT *) mem_alloc(sizeof(FCB), "nb");
cb = (FLOAT *) mem_alloc(sizeof(FCB), "cb");
ecb = (FLOAT *) mem_alloc(sizeof(FCB), "ecb");
bc = (FLOAT *) mem_alloc(sizeof(FCB), "bc");
wsamp_r = (FLOAT *) mem_alloc(sizeof(FBLK), "wsamp_r");
wsamp_i = (FLOAT *) mem_alloc(sizeof(FBLK), "wsamp_i");
phi = (FLOAT *) mem_alloc(sizeof(FBLK), "phi");
energy = (FLOAT *) mem_alloc(sizeof(FBLK), "energy");
c = (FLOAT *) mem_alloc(sizeof(FHBLK), "c");
fthr = (FLOAT *) mem_alloc(sizeof(FHBLK), "fthr");
snrtmp = (F32 *) mem_alloc(sizeof(F2_32), "snrtmp");
assert( lay == 3 );
if(init==0){
/* These dynamic memory allocations simulate "static" variables placed */
/* in the data space. Each mem_alloc() call here occurs only once at */
/* initialization time. The mem_free() function must not be called. */
numlines = (int *) mem_alloc(sizeof(ICB), "numlines");
partition = (int *) mem_alloc(sizeof(IHBLK), "partition");
cbval = (FLOAT *) mem_alloc(sizeof(FCB), "cbval");
rnorm = (FLOAT *) mem_alloc(sizeof(FCB), "rnorm");
window = (FLOAT *) mem_alloc(sizeof(FBLK), "window");
absthr = (FLOAT *) mem_alloc(sizeof(FHBLK), "absthr");
tmn = (double *) mem_alloc(sizeof(DCB), "tmn");
s = (FCB *) mem_alloc(sizeof(FCBCB), "s");
lthr = (FHBLK *) mem_alloc(sizeof(F2HBLK), "lthr");
r = (F2HBLK *) mem_alloc(sizeof(F22HBLK), "r");
phi_sav = (F2HBLK *) mem_alloc(sizeof(F22HBLK), "phi_sav");
/*#if 0 */
i = sfreq + 0.5;
switch(i){
case 32000: sfreq_idx = 0; break;
case 44100: sfreq_idx = 1; break;
case 48000: sfreq_idx = 2; break;
default: printf("error, invalid sampling frequency: %d Hz\n",i);
exit(-1);
}
printf("absthr[][] sampling frequency index: %d\n",sfreq_idx);
read_absthr(absthr, sfreq_idx);
switch(lay){
case 1: sync_flush=576; flush=384; syncsize=1024; break;
case 2: sync_flush=480; flush=576; syncsize=1056; break;
case 3: sync_flush=768; flush=576; syncsize=1344; break;
default: printf("Bad lay value:(%d)",lay); exit(-1); break;
}
/* #endif */
/* calculate HANN window coefficients */
/* for(i=0;i<BLKSIZE;i++) window[i] =0.5*(1-cos(2.0*PI*i/(BLKSIZE-1.0)));*/
for(i=0;i<BLKSIZE;i++) window[i] =0.5*(1-cos(2.0*PI*(i-0.5)/BLKSIZE));
for(i=0;i<BLKSIZE_s;i++)window_s[i]=0.5*(1-cos(2.0*PI*(i-0.5)/BLKSIZE_s));
/* reset states used in unpredictability measure */
for(i=0;i<HBLKSIZE;i++){
r[0][0][i]=r[1][0][i]=r[0][1][i]=r[1][1][i]=0;
phi_sav[0][0][i]=phi_sav[1][0][i]=0;
phi_sav[0][1][i]=phi_sav[1][1][i]=0;
lthr[0][i] = 60802371420160.0;
lthr[1][i] = 60802371420160.0;
}
/*****************************************************************************
* Initialization: Compute the following constants for use later *
* partition[HBLKSIZE] = the partition number associated with each *
* frequency line *
* cbval[CBANDS] = the center (average) bark value of each *
* partition *
* numlines[CBANDS] = the number of frequency lines in each partition *
* tmn[CBANDS] = tone masking noise *
*****************************************************************************/
/* compute fft frequency multiplicand */
freq_mult = sfreq/BLKSIZE;
/* calculate fft frequency, then bval of each line (use fthr[] as tmp storage)*/
for(i=0;i<HBLKSIZE;i++){
temp1 = i*freq_mult;
j = 1;
while(temp1>crit_band[j])j++;
fthr[i]=j-1+(temp1-crit_band[j-1])/(crit_band[j]-crit_band[j-1]);
}
partition[0] = 0;
/* temp2 is the counter of the number of frequency lines in each partition */
temp2 = 1;
cbval[0]=fthr[0];
bval_lo=fthr[0];
for(i=1;i<HBLKSIZE;i++){
if((fthr[i]-bval_lo)>0.33){
partition[i]=partition[i-1]+1;
cbval[partition[i-1]] = cbval[partition[i-1]]/temp2;
cbval[partition[i]] = fthr[i];
bval_lo = fthr[i];
numlines[partition[i-1]] = temp2;
temp2 = 1;
}
else {
partition[i]=partition[i-1];
cbval[partition[i]] += fthr[i];
temp2++;
}
}
numlines[partition[i-1]] = temp2;
cbval[partition[i-1]] = cbval[partition[i-1]]/temp2;
/************************************************************************
* Now compute the spreading function, s[j][i], the value of the spread-*
* ing function, centered at band j, for band i, store for later use *
************************************************************************/
for(j=0;j<CBANDS;j++){
for(i=0;i<CBANDS;i++){
temp1 = (cbval[i] - cbval[j])*1.05;
if(temp1>=0.5 && temp1<=2.5){
temp2 = temp1 - 0.5;
temp2 = 8.0 * (temp2*temp2 - 2.0 * temp2);
}
else temp2 = 0;
temp1 += 0.474;
temp3 = 15.811389+7.5*temp1-17.5*sqrt((double) (1.0+temp1*temp1));
if(temp3 <= -100) s[i][j] = 0;
else {
temp3 = (temp2 + temp3)*LN_TO_LOG10;
s[i][j] = exp(temp3);
}
}
}
/* Calculate Tone Masking Noise values */
for(j=0;j<CBANDS;j++){
temp1 = 15.5 + cbval[j];
tmn[j] = (temp1>24.5) ? temp1 : 24.5;
/* Calculate normalization factors for the net spreading functions */
rnorm[j] = 0;
for(i=0;i<CBANDS;i++){
rnorm[j] += s[j][i];
}
}
init++;
}
/************************* End of Initialization *****************************/
switch(lay) {
case 1:
case 2:
for ( i=0; i<lay; i++)
{
/*****************************************************************************
* Net offset is 480 samples (1056-576) for layer 2; this is because one must*
* stagger input data by 256 samples to synchronize psychoacoustic model with*
* filter bank outputs, then stagger so that center of 1024 FFT window lines *
* up with center of 576 "new" audio samples. *
* *
* For layer 1, the input data still needs to be staggered by 256 samples, *
* then it must be staggered again so that the 384 "new" samples are centered*
* in the 1024 FFT window. The net offset is then 576 and you need 448 "new"*
* samples for each iteration to keep the 384 samples of interest centered *
*****************************************************************************/
for (j=0; j<syncsize; j++)
{
if (j < (sync_flush) )
savebuf[j] = savebuf[j+flush];
else
savebuf[j] = *buffer++;
/**window data with HANN window***********************************************/
if (j<BLKSIZE)
{
wsamp_r[j] = window[j]*((FLOAT) savebuf[j]);
wsamp_i[j] = 0;
}
}
/**Compute FFT****************************************************************/
fft(wsamp_r,wsamp_i,energy,phi,1024);
/*****************************************************************************
* calculate the unpredictability measure, given energy[f] and phi[f] *
*****************************************************************************/
for(j=0; j<HBLKSIZE; j++){
r_prime = 2.0 * r[chn][old][j] - r[chn][oldest][j];
phi_prime = 2.0 * phi_sav[chn][old][j] - phi_sav[chn][oldest][j];
r[chn][new][j] = sqrt((double) energy[j]);
phi_sav[chn][new][j] = phi[j];
temp1 = r[chn][new][j] * cos((double) phi[j])
- r_prime * cos(phi_prime);
temp2=r[chn][new][j] * sin((double) phi[j])
- r_prime * sin(phi_prime);
temp3=r[chn][new][j] + fabs(r_prime);
if(temp3 != 0)c[j]=sqrt(temp1*temp1+temp2*temp2)/temp3;
else c[j] = 0;
}
/*only update data "age" pointers after you are done with the second channel */
/*for layer 1 computations, for the layer 2 double computations, the pointers*/
/*are reset automatically on the second pass */
if(lay==2 || chn==1){
if(new==0){new = 1; oldest = 1;}
else {new = 0; oldest = 0;}
if(old==0)old = 1; else old = 0;
}
/*****************************************************************************
* Calculate the grouped, energy-weighted, unpredictability measure, *
* grouped_c[], and the grouped energy. grouped_e[] *
*****************************************************************************/
for(j=1;j<CBANDS;j++){
grouped_e[j] = 0;
grouped_c[j] = 0;
}
grouped_e[0] = energy[0];
grouped_c[0] = energy[0]*c[0];
for(j=1;j<HBLKSIZE;j++){
grouped_e[partition[j]] += energy[j];
grouped_c[partition[j]] += energy[j]*c[j];
}
/*****************************************************************************
* convolve the grouped energy-weighted unpredictability measure *
* and the grouped energy with the spreading function, s[j][k] *
*****************************************************************************/
for(j=0;j<CBANDS;j++){
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -