A bunch of pointers marked as "restrict" to ease the job of the compiler
[opus.git] / libcelt / psy.c
index 0f7ffaf..7c7164b 100644 (file)
    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
 
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
 #include "psy.h"
 #include <math.h>
 #include "os_support.h"
+#include "arch.h"
 
 /* The Vorbis freq<->Bark mapping */
 #define toBARK(n)   (13.1f*atan(.00074f*(n))+2.24f*atan((n)*(n)*1.85e-8f)+1e-4f*(n))
 #define fromBARK(z) (102.f*(z)-2.f*pow(z,2.f)+.4f*pow(z,3.f)+pow(1.46f,z)-1.f)
 
-
+#ifndef STATIC_MODES
 /* Psychoacoustic spreading function. The idea here is compute a first order 
    recursive filter. The filter coefficient is frequency dependent and 
    chosen such that we have a -10dB/Bark slope on the right side and a -25dB/Bark
    slope on the left side. */
-void psydecay_init(struct PsyDecay *decay, int len, int Fs)
+void psydecay_init(struct PsyDecay *decay, int len, celt_int32_t Fs)
 {
    int i;
-   decay->decayR = celt_alloc(sizeof(float)*len);
-   decay->decayL = celt_alloc(sizeof(float)*len);
+   celt_word16_t *decayR = (celt_word16_t*)celt_alloc(sizeof(celt_word16_t)*len);
+   /*decay->decayL = celt_alloc(sizeof(celt_word16_t)*len);*/
    for (i=0;i<len;i++)
    {
       float f;
@@ -58,45 +63,52 @@ void psydecay_init(struct PsyDecay *decay, int len, int Fs)
       /* Back to FFT bin units */
       deriv *= Fs*(1/(2.f*len));
       /* decay corresponding to -10dB/Bark */
-      decay->decayR[i] = pow(.1f, deriv);
+      decayR[i] = Q15ONE*pow(.1f, deriv);
       /* decay corresponding to -25dB/Bark */
-      decay->decayL[i] = pow(0.0031623f, deriv);
-      //printf ("%f %f\n", decayL[i], decayR[i]);
+      /*decay->decayL[i] = Q15ONE*pow(0.0031623f, deriv);*/
+      /*printf ("%f %f\n", decayL[i], decayR[i]);*/
    }
+   decay->decayR = decayR;
 }
 
 void psydecay_clear(struct PsyDecay *decay)
 {
-   celt_free(decay->decayR);
-   celt_free(decay->decayL);
+   celt_free((celt_word16_t *)decay->decayR);
+   /*celt_free(decay->decayL);*/
 }
+#endif
 
-static void spreading_func(struct PsyDecay *d, float *psd, float *mask, int len, int Fs)
+static void spreading_func(const struct PsyDecay *d, celt_word32_t * restrict psd, int len)
 {
    int i;
-   float mem;
-   //for (i=0;i<len;i++) printf ("%f ", psd[i]);
+   celt_word32_t mem;
+   /*for (i=0;i<len;i++) printf ("%f ", psd[i]);*/
    /* Compute right slope (-10 dB/Bark) */
    mem=psd[0];
    for (i=0;i<len;i++)
    {
-      mask[i] = (1-d->decayR[i])*psd[i] + d->decayR[i]*mem;
-      mem = mask[i];
+      /* psd = (1-decay)*psd + decay*mem */
+      psd[i] = EPSILON + psd[i] + MULT16_32_Q15(d->decayR[i],mem-psd[i]);
+      mem = psd[i];
    }
    /* Compute left slope (-25 dB/Bark) */
-   mem=mask[len-1];
+   mem=psd[len-1];
    for (i=len-1;i>=0;i--)
    {
-      mask[i] = (1-d->decayR[i])*mask[i] + d->decayL[i]*mem;
-      mem = mask[i];
+      /* Left side has around twice the slope as the right side, so we just
+         square the coef instead of storing two sets of decay coefs */
+      celt_word16_t decayL = MULT16_16_Q15(d->decayR[i], d->decayR[i]);
+      /* psd = (1-decay)*psd + decay*mem */
+      psd[i] = EPSILON + psd[i] + MULT16_32_Q15(decayL,mem-psd[i]);
+      mem = psd[i];
    }
-   //for (i=0;i<len;i++) printf ("%f ", mask[i]); printf ("\n");
+   /*for (i=0;i<len;i++) printf ("%f ", mask[i]); printf ("\n");*/
 #if 0 /* Prints signal and mask energy per critical band */
    for (i=0;i<25;i++)
    {
       int start,end;
       int j;
-      float Esig=0, Emask=0;
+      celt_word32_t Esig=0, Emask=0;
       start = (int)floor(fromBARK((float)i)*(2*len)/Fs);
       if (start<0)
          start = 0;
@@ -117,34 +129,36 @@ static void spreading_func(struct PsyDecay *d, float *psd, float *mask, int len,
 }
 
 /* Compute a marking threshold from the spectrum X. */
-void compute_masking(struct PsyDecay *decay, float *X, float *mask, int len, int Fs)
+void compute_masking(const struct PsyDecay *decay, celt_word16_t *X, celt_mask_t * restrict mask, int len)
 {
    int i;
-   int N=len/2;
-   float psd[N];
-   psd[0] = X[0]*X[0];
+   int N;
+   N=len>>1;
+   mask[0] = MULT16_16(X[0], X[0]);
    for (i=1;i<N;i++)
-      psd[i] = X[i*2]*X[i*2] + X[i*2+1]*X[i*2+1];
+      mask[i] = ADD32(MULT16_16(X[i*2], X[i*2]), MULT16_16(X[i*2+1], X[i*2+1]));
    /* TODO: Do tone masking */
    /* Noise masking */
-   spreading_func(decay, psd, mask, N, Fs);
-   
+   spreading_func(decay, mask, N);
 }
 
-void compute_mdct_masking(struct PsyDecay *decay, float *X, float *mask, int len, int Fs)
+#if 0 /* Not needed for now, but will be useful in the future */
+void compute_mdct_masking(const struct PsyDecay *decay, celt_word32_t *X, celt_mask_t *mask, int len)
 {
    int i;
-   float psd[len];
-   float mem;
+   VARDECL(float *psd);
+   SAVE_STACK;
+   ALLOC(psd, len, float);
    for (i=0;i<len;i++)
       mask[i] = X[i]*X[i];
    for (i=1;i<len-1;i++)
       psd[i] = .5*mask[i] + .25*(mask[i-1]+mask[i+1]);
-   //psd[0] = .5*mask[0]+.25*(mask[1]+mask[2]);
+   /*psd[0] = .5*mask[0]+.25*(mask[1]+mask[2]);*/
    psd[0] = .5*mask[0]+.5*mask[1];
    psd[len-1] = .5*(mask[len-1]+mask[len-2]);
    /* TODO: Do tone masking */
    /* Noise masking */
-   spreading_func(decay, psd, mask, len, Fs);
-   
+   spreading_func(decay, psd, mask, len);
+   RESTORE_STACK;  
 }
+#endif