1 /********************************************************************
3 * THIS FILE IS PART OF THE Ogg Vorbis SOFTWARE CODEC SOURCE CODE. *
4 * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
5 * THE GNU PUBLIC LICENSE 2, WHICH IS INCLUDED WITH THIS SOURCE. *
6 * PLEASE READ THESE TERMS DISTRIBUTING. *
8 * THE OggSQUISH SOURCE CODE IS (C) COPYRIGHT 1994-1999 *
9 * by 1999 Monty <monty@xiph.org> and The XIPHOPHORUS Company *
10 * http://www.xiph.org/ *
12 ********************************************************************
14 function: random psychoacoustics (not including preecho)
15 author: Monty <xiphmont@mit.edu>
16 modifications by: Monty
17 last modification date: Oct 15 1999
19 ********************************************************************/
31 /* Set up decibel threshhold 'curves'. Actually, just set a level at
32 log frequency intervals, interpolate, and call it a curve. */
34 static void set_curve(double *points,
35 double *ref,int rn,double rrate,
36 double *c,int n, double crate){
40 int endpos=points[i+1]*n*rrate/crate;
42 double delta=(ref[i+1]-base)/(endpos-j);
43 for(;j<endpos && j<n;j++){
50 void _vp_psy_init(psy_lookup *p,vorbis_info *vi,int n){
51 memset(p,0,sizeof(psy_lookup));
52 p->noisethresh=malloc(n*sizeof(double));
53 p->maskthresh=malloc(n*sizeof(double));
57 /* set up the curves for a given blocksize and sample rate */
59 set_curve(vi->threshhold_points,
60 vi->noisethresh,THRESH_POINTS,48000,p->noisethresh,n,vi->rate);
61 set_curve(vi->threshhold_points,
62 vi->maskthresh, THRESH_POINTS,48000,p->maskthresh, n,vi->rate);
70 sprintf(buffer,"noise_threshhold_%d.m",n);
71 out=fopen(buffer,"w+");
73 fprintf(out,"%g\n",p->noisethresh[j]);
75 sprintf(buffer,"mask_threshhold_%d.m",n);
76 out=fopen(buffer,"w+");
78 fprintf(out,"%g\n",p->maskthresh[j]);
85 void _vp_psy_clear(psy_lookup *p){
87 if(p->noisethresh)free(p->noisethresh);
88 if(p->maskthresh)free(p->maskthresh);
89 memset(p,0,sizeof(psy_lookup));
93 /* Find the mean log energy of a given 'band'; used to evaluate tones
94 against background noise */
96 /* This is faster than a real convolution, gives us roughly the log f
97 scale we seek, and gives OK results. So, that means it's a good
100 void _vp_noise_floor(psy_lookup *p, double *f, double *m){
102 vorbis_info *vi=p->vi;
109 double bias=n*vi->noisebias;
112 long newlo=i*vi->lnoise-bias;
113 long newhi=i*vi->hnoise+bias;
119 for(j=hi;j<newhi;j++){
123 for(j=lo;j<newlo;j++){
131 /* attenuate by the noise threshhold curve */
132 temp=fromdB(acc/div+p->noisethresh[i]);
133 if(m[i]<temp)m[i]=temp;
137 /* Masking curve: linear rolloff on a dB scale, adjusted by octave,
138 attenuated by maskthresh */
140 void _vp_mask_floor(psy_lookup *p,double *f, double *m){
142 double hroll=p->vi->hroll;
143 double lroll=p->vi->lroll;
144 double ocSCALE=1./log(2);
145 double curmask=-9.e40;
146 double maskbias=n*p->vi->maskbias;
147 double curoc=log(maskbias)*ocSCALE;
150 /* run mask forward then backward */
152 double newmask=todB(f[i])+p->maskthresh[i];
153 double newoc=log(i+maskbias)*ocSCALE;
154 double roll=curmask-(newoc-curoc)*hroll;
157 roll=curmask=newmask;
161 if(m[i]<troll)m[i]=troll;
165 curoc=log(n+maskbias)*ocSCALE;
167 double newmask=todB(f[i])+p->maskthresh[i];
168 double newoc=log(i+maskbias)*ocSCALE;
169 double roll=curmask-(curoc-newoc)*lroll;
172 roll=curmask=newmask;
176 if(m[i]<troll)m[i]=troll;
180 /* s must be padded at the end with m-1 zeroes */
181 static void time_convolve(double *s,double *r,int n,int m){
189 acc+=s[i+j]*r[m-j-1];
195 /*************************************************************************/
196 /* Continuous balance analysis/synthesis *********************************/
199 /* Compute the average continual spectral balance of the given vectors
200 (in radians); encode into LPC coefficients */
202 double _vp_balance_compute(double *A, double *B, double *lpc,lpc_lookup *vb){
203 /* correlate in time (the response function is small). Log
204 frequency scale, small mapping */
209 /* 256/15 are arbitrary but unimportant to decoding */
213 /* This is encode side. Don't think too hard about it */
215 double workA[mapped+resp];
216 double workB[mapped+resp];
220 memset(workA,0,sizeof(workA));
221 memset(workB,0,sizeof(workB));
222 memset(workC,0,sizeof(workC));
225 int j=vb->bscale[i]+(resp>>1);
226 double mag_sq=A[i]*A[i]+B[i]*B[i];
233 if((A[i]<0 && B[i]>0)||
235 /* rotate II and IV into the first quadrant. III is already there */
240 workA[j]+=mag_sq*sin(phi);
241 workB[j]+=mag_sq*cos(phi);
244 /* prepare convolution vector. Play with a few different shapes */
247 workC[i]=sin(M_PI*(i+1)/(float)(resp+1));
251 time_convolve(workA,workC,mapped,resp);
252 time_convolve(workB,workC,mapped,resp);
257 for(i=0;i<mapped;i++){
258 p[i]=atan2(workA[i],workB[i]);
261 amp1=sqrt(vorbis_gen_lpc(p,lpc,vb));
268 /*void _vp_balance_apply(double *A, double *B, double *lpc, double amp,
269 lpc_lookup *vb,int divp){
271 for(i=0;i<vb->n;i++){
272 double mag=sqrt(A[i]*A[i]+B[i]*B[i]);
273 double del=vorbis_lpc_magnitude(vb->dscale[i],lpc,vb->m)*amp;
274 double phi=atan2(A[i],B[i]);