From 46014ca49a91e345cd4d3a9d70a67a1cbe73c2ac Mon Sep 17 00:00:00 2001
From: Jean-Marc Valin <Jean-Marc.Valin@csiro.au>
Date: Fri, 14 Dec 2007 13:47:04 +1100
Subject: [PATCH] Propagating perceptual weighting around (not used yet).

---
 libcelt/bands.c    |  6 +++---
 libcelt/bands.h    |  2 +-
 libcelt/celt.c     | 11 ++++++++++-
 libcelt/testcelt.c |  4 ++--
 libcelt/vq.c       |  6 +++---
 libcelt/vq.h       |  4 ++--
 6 files changed, 21 insertions(+), 12 deletions(-)

diff --git a/libcelt/bands.c b/libcelt/bands.c
index b1a75803..f6b460e2 100644
--- a/libcelt/bands.c
+++ b/libcelt/bands.c
@@ -214,7 +214,7 @@ void pitch_quant_bands(const CELTMode *m, float *X, float *P, float *gains)
       P[i] = 0;
 }
 
-void quant_bands(const CELTMode *m, float *X, float *P, ec_enc *enc)
+void quant_bands(const CELTMode *m, float *X, float *P, float *W, ec_enc *enc)
 {
    int i, j, B;
    const int *eBands = m->eBands;
@@ -227,14 +227,14 @@ void quant_bands(const CELTMode *m, float *X, float *P, ec_enc *enc)
       q = m->nbPulses[i];
       if (q>0) {
          float n = sqrt(B*(eBands[i+1]-eBands[i]));
-         alg_quant(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, P+B*eBands[i], 0.7, enc);
+         alg_quant(X+B*eBands[i], W+B*eBands[i], B*(eBands[i+1]-eBands[i]), q, P+B*eBands[i], 0.7, enc);
          for (j=B*eBands[i];j<B*eBands[i+1];j++)
             norm[j] = X[j] * n;
          //printf ("%f ", log2(ncwrs64(B*(eBands[i+1]-eBands[i]), q))/(B*(eBands[i+1]-eBands[i])));
          //printf ("%f ", log2(ncwrs64(B*(eBands[i+1]-eBands[i]), q)));
       } else {
          float n = sqrt(B*(eBands[i+1]-eBands[i]));
-         copy_quant(X+B*eBands[i], B*(eBands[i+1]-eBands[i]), -q, norm, B, eBands[i], enc);
+         copy_quant(X+B*eBands[i], W+B*eBands[i], B*(eBands[i+1]-eBands[i]), -q, norm, B, eBands[i], enc);
          for (j=B*eBands[i];j<B*eBands[i+1];j++)
             norm[j] = X[j] * n;
          //printf ("%f ", (1+log2(eBands[i]-(eBands[i+1]-eBands[i]))+log2(ncwrs64(B*(eBands[i+1]-eBands[i]), -q)))/(B*(eBands[i+1]-eBands[i])));
diff --git a/libcelt/bands.h b/libcelt/bands.h
index 01d69fcd..5a195764 100644
--- a/libcelt/bands.h
+++ b/libcelt/bands.h
@@ -47,7 +47,7 @@ void compute_pitch_gain(const CELTMode *m, float *X, float *P, float *gains, flo
 
 void pitch_quant_bands(const CELTMode *m, float *X, float *P, float *gains);
 
-void quant_bands(const CELTMode *m, float *X, float *P, ec_enc *enc);
+void quant_bands(const CELTMode *m, float *X, float *P, float *W, ec_enc *enc);
 
 void unquant_bands(const CELTMode *m, float *X, float *P, ec_dec *dec);
 
diff --git a/libcelt/celt.c b/libcelt/celt.c
index 0bcf53be..ae7f306f 100644
--- a/libcelt/celt.c
+++ b/libcelt/celt.c
@@ -40,6 +40,7 @@
 #include "probenc.h"
 #include "quant_pitch.h"
 #include "quant_bands.h"
+#include "psy.h"
 
 #define MAX_PERIOD 1024
 
@@ -50,6 +51,7 @@ struct CELTEncoder {
    int block_size;
    int nb_blocks;
    int channels;
+   int Fs;
    
    ec_byte_buffer buf;
    ec_enc         enc;
@@ -83,6 +85,7 @@ CELTEncoder *celt_encoder_new(const CELTMode *mode)
    st->frame_size = B*N;
    st->block_size = N;
    st->nb_blocks  = B;
+   st->Fs = 44100;
    
    ec_byte_writeinit(&st->buf);
    ec_enc_init(&st->enc,&st->buf);
@@ -207,6 +210,7 @@ int celt_encode(CELTEncoder *st, short *pcm)
    
    float X[B*C*N];         /**< Interleaved signal MDCTs */
    float P[B*C*N];         /**< Interleaved pitch MDCTs*/
+   float mask[B*C*N];      /**< Masking curve */
    float bandE[st->mode->nbEBands];
    float gains[st->mode->nbPBands];
    int pitch_index;
@@ -228,6 +232,11 @@ int celt_encode(CELTEncoder *st, short *pcm)
    /* Compute MDCTs */
    compute_mdcts(&st->mdct_lookup, st->window, in, X, N, B, C);
    
+   compute_masking(X, mask, B*C*N, st->Fs);
+   /* Invert and stretch the mask to length of X */
+   for (i=B*C*N-1;i>=0;i--)
+      mask[i] = 1/(1+mask[i>>1]);
+   
    /* Pitch analysis */
    for (c=0;c<C;c++)
    {
@@ -286,7 +295,7 @@ int celt_encode(CELTEncoder *st, short *pcm)
       sum += X[i]*X[i];
    printf ("%f\n", sum);*/
    /* Residual quantisation */
-   quant_bands(st->mode, X, P, &st->enc);
+   quant_bands(st->mode, X, P, mask, &st->enc);
    
    if (0) {//This is just for debugging
       ec_enc_done(&st->enc);
diff --git a/libcelt/testcelt.c b/libcelt/testcelt.c
index a4f50ec8..943c7e5b 100644
--- a/libcelt/testcelt.c
+++ b/libcelt/testcelt.c
@@ -52,8 +52,8 @@ int main(int argc, char *argv[])
    outFile = argv[2];
    fout = fopen(outFile, "wb+");
    
-   enc = celt_encoder_new(celt_mode2);
-   dec = celt_decoder_new(celt_mode2);
+   enc = celt_encoder_new(celt_mode1);
+   dec = celt_decoder_new(celt_mode1);
    
    while (!feof(fin))
    {
diff --git a/libcelt/vq.c b/libcelt/vq.c
index 40735875..46899c02 100644
--- a/libcelt/vq.c
+++ b/libcelt/vq.c
@@ -38,7 +38,7 @@
 /* Improved algebraic pulse-base quantiser. The signal x is replaced by the sum of the pitch 
    a combination of pulses such that its norm is still equal to 1. The only difference with 
    the quantiser above is that the search is more complete. */
-void alg_quant(float *x, int N, int K, float *p, float alpha, ec_enc *enc)
+void alg_quant(float *x, float *W, int N, int K, float *p, float alpha, ec_enc *enc)
 {
    int L = 3;
    //float tata[200];
@@ -223,7 +223,7 @@ void alg_quant(float *x, int N, int K, float *p, float alpha, ec_enc *enc)
 static const float pg[5] = {1.f, .6f, .45f, 0.35f, 0.25f};
 
 /* Finds the right offset into Y and copy it */
-void copy_quant(float *x, int N, int K, float *Y, int B, int N0, ec_enc *enc)
+void copy_quant(float *x, float *W, int N, int K, float *Y, int B, int N0, ec_enc *enc)
 {
    int i,j;
    int best=0;
@@ -287,7 +287,7 @@ void copy_quant(float *x, int N, int K, float *Y, int B, int N0, ec_enc *enc)
       E = .8/sqrt(E);
       for (j=0;j<N;j++)
          P[j] *= E;
-      alg_quant(x, N, K, P, 0, enc);
+      alg_quant(x, W, N, K, P, 0, enc);
    }
 }
 
diff --git a/libcelt/vq.h b/libcelt/vq.h
index ae75f7e9..2e4f5e01 100644
--- a/libcelt/vq.h
+++ b/libcelt/vq.h
@@ -39,12 +39,12 @@
 /* Algebraic pulse-base quantiser. The signal x is replaced by the sum of the pitch 
    a combination of pulses such that its norm is still equal to 1. The only difference with 
    the quantiser above is that the search is more complete. */
-void alg_quant(float *x, int N, int K, float *p, float alpha, ec_enc *enc);
+void alg_quant(float *x, float *W, int N, int K, float *p, float alpha, ec_enc *enc);
 
 void alg_unquant(float *x, int N, int K, float *p, float alpha, ec_dec *dec);
 
 /* Finds the right offset into Y and copy it */
-void copy_quant(float *x, int N, int K, float *Y, int B, int N0, ec_enc *enc);
+void copy_quant(float *x, float *W, int N, int K, float *Y, int B, int N0, ec_enc *enc);
 
 void copy_unquant(float *x, int N, int K, float *Y, int B, int N0, ec_dec *dec);
 #endif /* VQ_H */