📄 tonal.c
字号:
/**********************************************************************
* ISO MPEG Audio Subgroup Software Simulation Group (1996)
* ISO 13818-3 MPEG-2 Audio Multichannel Encoder
*
* $Id: tonal.c 1.8 1996/02/12 07:13:35 rowlands Exp $
*
* $Log: tonal.c $
* Revision 1.8 1996/02/12 07:13:35 rowlands
* Release following Munich meeting
*
* Revision 1.5.2.1 1995/11/06 04:19:12 rowlands
* Received from Uwe Felderhoff (IRT)
*
* Revision 1.7 1995/08/14 08:06:37 tenkate
* ML-LSF added Warner ten Kate 7/8/95 (Philips)
* II_psycho_one() split into II_psycho_one() for audio data and
* II_psycho_one_ml() for MultiLingual data.
* Variables crit_band, cbound and sub_size has been made local to
* these functions.
* Tables "cb" and "th" are copied from LSF-directory.
* ltg has found its ml counterpart.
*
* Revision 1.6 1995/07/31 07:47:50 tenkate
* bugs correction (if center==1 in Psycho_One), 25/07/95 WtK
*
* Revision 1.4.2.1 1995/06/16 03:46:42 rowlands
* Input from Susanne Ritscher (IRT)
*
**********************************************************************/
/**********************************************************************
* date programmers comment *
* 2/25/91 Douglas Wong start of version 1.1 records *
* 3/06/91 Douglas Wong rename: setup.h to endef.h *
* updated I_psycho_one and II_psycho_one*
* 3/11/91 W. J. Carter Added Douglas Wong's updates dated *
* 3/9/91 for I_Psycho_One() and for *
* II_Psycho_One(). *
* 5/10/91 W. Joseph Carter Ported to Macintosh and Unix. *
* Located and fixed numerous software *
* bugs and table data errors. *
* 6/11/91 Davis Pan corrected several bugs *
* based on comments from H. Fuchs *
* 01jul91 dpwe (Aware Inc.) Made pow() args float *
* Removed logical bug in I_tonal_label: *
* Sometimes *tone returned == STOP *
* 7/10/91 Earle Jennings no change necessary in port to MsDos *
* 11sep91 dpwe@aware.com Subtracted 90.3dB from II_f_f_t peaks *
* 10/1/91 Peter W. Farrett Updated II_Psycho_One(),I_Psycho_One()*
* to include comments. *
*11/29/91 Masahiro Iwadare Bug fix regarding POWERNORM *
* fixed several other miscellaneous bugs*
* 2/11/92 W. Joseph Carter Ported new code to Macintosh. Most *
* important fixes involved changing *
* 16-bit ints to long or unsigned in *
* bit alloc routines for quant of 65535 *
* and passing proper function args. *
* Removed "Other Joint Stereo" option *
* and made bitrate be total channel *
* bitrate, irrespective of the mode. *
* Fixed many small bugs & reorganized. *
* 2/12/92 Masahiro Iwadare Fixed some potential bugs in *
* Davis Pan subsampling() *
* 2/25/92 Masahiro Iwadare Fixed some more potential bugs *
* 92-11-06 Soren H. Nielsen Corrected power calculation in I_ and *
* II_f_f_t. *
**********************************************************************
* *
* *
* MPEG/audio Phase 2 coding/decoding multichannel *
* *
* 7/27/93 Susanne Ritscher, IRT Munich *
* 8/27/93 Susanne Ritscher, IRT Munich *
* Channel-Switching is working *
* 9/1/93 Susanne Ritscher, IRT Munich *
* all channels normalized *
* *
* 9/20/93 channel-switching is only performed at a *
* certain limit of TC_ALLOC dB, which is included *
* in encoder.h *
* *
* 10/18/93 seperated smr and ltmin *
* *
* Version 1.0 *
* *
* 07/12/94 Susanne Ritscher, IRT Munich *
* Tel: +49 89 32399 458 *
* Fax: +49 89 32399 415 *
* *
* Version 1.1 *
* *
* 02/23/95 Susanne Ritscher, IRT Munich *
* corrected some bugs *
* extension bitstream is working *
* *
**********************************************************************/
#define TONAL_WIDTH 0.1 /* When more than one tonal component is within
this width in Bark, the weaker one(s) are
eliminated */
#include "common.h"
#include "encoder.h"
/**********************************************************************/
/*
/* This module implements the psychoacoustic model I for the
/* MPEG encoder layer II. It uses simplified tonal and noise masking
/* threshold analysis to generate SMR for the encoder bit allocation
/* routine.
/*
/**********************************************************************/
int read_crit_band(int lay, int freq)
{
int crit_band;
FILE *fp;
char r[16], t[80];
strcpy(r, "2cb1");
r[0] = (char) lay + '0';
r[3] = (char) freq + '0';
if( !(fp = OpenTableFile(r)) ){ /* check boundary values */
printf("Please check %s boundary table\n",r);
exit(0);
}
fgets(t,80,fp);
sscanf(t,"%d\n",&crit_band);
fclose(fp);
return(crit_band);
}
void read_cbound(int lay, int freq, int crit_band, int *cbound)
/* this function reads in critical band boundaries */
{
int i,j,k;
FILE *fp;
char r[16], t[80];
strcpy(r, "2cb1");
r[0] = (char) lay + '0';
r[3] = (char) freq + '0';
if( !(fp = OpenTableFile(r)) ){ /* check boundary values */
printf("Please check %s boundary table\n",r);
exit(0);
}
fgets(t,80,fp); /* skip input for critical bands */
sscanf(t,"%d\n",&k);
for(i=0;i<crit_band;i++){
fgets(t,80,fp);
sscanf(t,"%d %d\n",&j, &k);
if(i==j) cbound[j] = k;
else { /* error */
printf("Please check index %d in cbound table %s\n",i,r);
exit(0);
}
}
fclose(fp);
}
void read_freq_band(int *sub_size, g_ptr *ltg, int lay, int freq)
/* this function reads in frequency bands and bark values */
{
int i,j, k;
double a,b,c;
FILE *fp;
char r[16], t[80];
strcpy(r, "2th1");
r[0] = (char) lay + '0';
r[3] = (char) freq + '0';
if( !(fp = OpenTableFile(r)) ){ /* check freq. values */
printf("Please check frequency and cband table %s\n",r);
exit(0);
}
fgets(t,80,fp); /* read input for freq. subbands */
sscanf(t,"%d\n",sub_size);
*ltg = (g_ptr) mem_alloc (sizeof(g_thres) * (*sub_size), "ltg");
(*ltg)[0].line = 0; /* initialize global masking threshold */
(*ltg)[0].bark = 0;
(*ltg)[0].hear = 0;
for(i=1;i<(*sub_size);i++){ /* continue to read freq. subband */
fgets(t,80,fp); /* and assign */
sscanf(t,"%d %d %lf %lf\n",&j, &k, &b, &c);
if(i == j){
(*ltg)[j].line = k;
(*ltg)[j].bark = b;
(*ltg)[j].hear = c;
}
else { /* error */
printf("Please check index %d in freq-cb table %s\n",i,r);
exit(0);
}
}
fclose(fp);
}
void make_map (int sub_size, mask *power, g_thres *ltg)
/* this function calculates the global masking threshold */
{
int i, j;
for (i = 1; i < sub_size; i++)
for (j = ltg[i-1].line; j <= ltg[i].line; j++)
power[j].map = i;
}
double add_db(double a, double b)
{
a = pow(10.0,a/10.0);
b = pow(10.0,b/10.0);
return 10 * log10(a+b);
}
/****************************************************************/
/*
/* Fast Fourier transform of the input samples.
/*
/****************************************************************/
void II_f_f_t (double *sample, mask *power) /* this function calculates an */
/* FFT analysis for the freq. */
/* domain */
{
int i,j,k,ll,l=0;
int ip, le, le1;
double t_r, t_i, u_r, u_i;
static int M, MM1, init = 0, N, NV2, NM1;
double *x_r, *x_i, *energy;
static int *rev;
static double *w_r, *w_i;
x_r = (double *) mem_alloc (sizeof (DFFT), "x_r");
x_i = (double *) mem_alloc (sizeof (DFFT), "x_i");
energy = (double *) mem_alloc (sizeof (DFFT), "energy");
if (!init)
{
rev = (int *) mem_alloc (sizeof (IFFT), "rev");
w_r = (double *) mem_alloc (sizeof (D10), "w_r");
w_i = (double *) mem_alloc (sizeof (D10), "w_i");
M = 10;
MM1 = 9;
N = FFT_SIZE;
NV2 = FFT_SIZE >> 1;
NM1 = FFT_SIZE - 1;
for (ll = 0; ll < M; ll++)
{
le = 1 << (M-ll);
le1 = le >> 1;
w_r[ll] = cos (PI / le1);
w_i[ll] = -sin (PI / le1);
}
for (i = 0; i < FFT_SIZE; rev[i] = l, i++)
for (j = 0, l = 0; j < 10; j++)
{
k = (i >> j) & 1;
l |= (k << 9 - j);
}
init = 1;
}
for (i = 0; i < FFT_SIZE; i++)
{
x_r[i] = sample[i];
x_i[i] = energy[i] = 0;
}
/*
memcpy ((char *) x_r, (char *) sample, sizeof (double) * FFT_SIZE);
*/
for (ll = 0; ll < MM1; ll++)
{
le = 1 << (M - ll);
le1 = le >> 1;
u_r = 1;
u_i = 0;
for (j = 0; j < le1; j++)
{
for (i = j; i < N; i += le)
{
ip = i + le1;
t_r = x_r[i] + x_r[ip];
t_i = x_i[i] + x_i[ip];
x_r[ip] = x_r[i] - x_r[ip];
x_i[ip] = x_i[i] - x_i[ip];
x_r[i] = t_r;
x_i[i] = t_i;
t_r = x_r[ip];
x_r[ip] = x_r[ip] * u_r - x_i[ip] * u_i;
x_i[ip] = x_i[ip] * u_r + t_r * u_i;
}
t_r = u_r;
u_r = u_r * w_r[ll] - u_i * w_i[ll];
u_i = u_i * w_r[ll] + t_r * w_i[ll];
}
}
for (i = 0; i < N; i += 2)
{
ip = i + 1;
t_r = x_r[i] + x_r[ip];
t_i = x_i[i] + x_i[ip];
x_r[ip] = x_r[i] - x_r[ip];
x_i[ip] = x_i[i] - x_i[ip];
x_r[i] = t_r;
x_i[i] = t_i;
energy[i] = x_r[i] * x_r[i] + x_i[i] * x_i[i];
}
for (i = 0; i < FFT_SIZE; i++)
if (i < rev[i])
{
t_r = energy[i];
energy[i] = energy[rev[i]];
energy[rev[i]] = t_r;
}
for (i = 0; i < HAN_SIZE; i++)
{
/* calculate power density spectrum */
if (energy[i] < 1E-20)
energy[i] = 1E-20;
/* power calculation corrected with a factor 4, both positive
and negative frequencies exist, 1992-11-06 shn */
power[i].x = 10 * log10 (energy[i] * 4.0) + POWERNORM;
power[i].next = STOP;
power[i].type = FALSE;
}
mem_free ((void **) &x_r);
mem_free ((void **) &x_i);
mem_free ((void **) &energy);
}
/****************************************************************/
/*
/* Window the incoming audio signal.
/*
/****************************************************************/
void II_hann_win (double *sample) /* this function calculates a */
/* Hann window for PCM (input) */
{ /* samples for a 1024-pt. FFT */
register int i;
register double sqrt_8_over_3;
static int init = 0;
static double *window;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -