diff --git a/libcelt/Makefile.am b/libcelt/Makefile.am
index 5c0780c3e8d6e16656cb798b84baa0b924714596..63616e949104d17af72cb2fb867fac7e763e1e5a 100644
--- a/libcelt/Makefile.am
+++ b/libcelt/Makefile.am
@@ -10,14 +10,14 @@ lib_LTLIBRARIES = libcelt.la
# Sources for compilation in the library
libcelt_la_SOURCES = bands.c celt.c cwrs.c fftwrap.c mdct.c modes.c pitch.c \
- quant_bands.c quant_pitch.c smallft.c vq.c
+ psy.c quant_bands.c quant_pitch.c smallft.c vq.c
#noinst_HEADERS =
libcelt_la_LDFLAGS = -version-info @CELT_LT_CURRENT@:@CELT_LT_REVISION@:@CELT_LT_AGE@
noinst_HEADERS = arch.h bands.h celt.h cwrs.h fftwrap.h mdct.h modes.h \
- os_support.h pitch.h quant_bands.h quant_pitch.h smallft.h vq.h
+ os_support.h pitch.h psy.h quant_bands.h quant_pitch.h smallft.h vq.h
noinst_PROGRAMS = testcelt
testcelt_SOURCES = testcelt.c
diff --git a/libcelt/pitch.c b/libcelt/pitch.c
index 44ff08ad11e6938f8cb6d337adb061ce6d1af669..decbac9b23ba35c05b8e844d36c3afc5e11bca4b 100644
--- a/libcelt/pitch.c
+++ b/libcelt/pitch.c
@@ -25,6 +25,7 @@
#include <math.h>
#include "fftwrap.h"
#include "pitch.h"
+#include "psy.h"
void find_spectral_pitch(void *fft, float *x, float *y, int lag, int len, int *pitch)
{
@@ -42,14 +43,8 @@ void find_spectral_pitch(void *fft, float *x, float *y, int lag, int len, int *p
spx_fft(fft, xx, X);
spx_fft(fft, y, Y);
- curve[0] = 1;
- for (i=1;i<n2;i++)
- {
- curve[i] = sqrt((X[2*i-1]*X[2*i-1] + X[2*i ]*X[2*i ])*(Y[2*i-1]*Y[2*i-1] + Y[2*i ]*Y[2*i ]));
- curve[i] = curve[i]+.7*curve[i];
- }
- for (i=n2-2;i>=0;i--)
- curve[i] = curve[i] + .7*curve[i+1];
+
+ compute_masking(X, curve, lag, 44100);
X[0] = 0;
for (i=1;i<lag/2;i++)
@@ -57,9 +52,8 @@ void find_spectral_pitch(void *fft, float *x, float *y, int lag, int len, int *p
float n;
//n = 1.f/(1e1+sqrt(sqrt((X[2*i-1]*X[2*i-1] + X[2*i ]*X[2*i ])*(Y[2*i-1]*Y[2*i-1] + Y[2*i ]*Y[2*i ]))));
//n = 1;
- n = 1.f/pow(1+curve[i],.8)/(i+60);
- //if (i>lag/6)
- // n *= .5;
+ n = 1.f/pow(1+curve[i],.5)/(i+60);
+ //n = 1.f/(1+curve[i]);
float tmp = X[2*i-1];
X[2*i-1] = (X[2*i-1]*Y[2*i-1] + X[2*i ]*Y[2*i ])*n;
X[2*i ] = (- X[2*i ]*Y[2*i-1] + tmp*Y[2*i ])*n;
diff --git a/libcelt/psy.c b/libcelt/psy.c
new file mode 100644
index 0000000000000000000000000000000000000000..51da3665ffa1449c87f66571db1c5fadeeacaa0d
--- /dev/null
+++ b/libcelt/psy.c
@@ -0,0 +1,87 @@
+/* (C) 2007 Jean-Marc Valin, CSIRO
+*/
+/*
+ Redistribution and use in source and binary forms, with or without
+ modification, are permitted provided that the following conditions
+ are met:
+
+ - Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+ - Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer in the
+ documentation and/or other materials provided with the distribution.
+
+ - Neither the name of the Xiph.org Foundation nor the names of its
+ contributors may be used to endorse or promote products derived from
+ this software without specific prior written permission.
+
+ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR
+ CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+ LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+#include "psy.h"
+#include <math.h>
+
+/* Psychoacoustic spreading function. The idea here is compute a first order
+ recursive smoothing. 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. */
+static void spreading_func(float *psd, float *mask, int len, int Fs)
+{
+ int i;
+ float decayL[len], decayR[len];
+ float mem;
+ //for (i=0;i<len;i++) printf ("%f ", psd[i]);
+ /* This can easily be tabulated, which makes the function very fast. */
+ for (i=0;i<len;i++)
+ {
+ float f;
+ float deriv;
+ f = Fs*i*(1/(2.f*len));
+ deriv = (8.288e-8 * f)/(3.4225e-16 *f*f*f*f + 1) + .009694/(5.476e-7 *f*f + 1) + 1e-4;
+ deriv *= Fs*(1/(2.f*len));
+ decayR[i] = pow(.1f, deriv);
+ decayL[i] = pow(0.0031623f, deriv);
+ }
+ /* Compute right slope (-10 dB/Bark) */
+ mem=psd[0];
+ for (i=0;i<len;i++)
+ {
+ mask[i] = (1-decayR[i])*psd[i] + decayR[i]*mem;
+ mem = mask[i];
+ }
+ /* Compute left slope (-25 dB/Bark) */
+ mem=mask[len-1];
+ for (i=len-1;i>=0;i--)
+ {
+ mask[i] = (1-decayR[i])*mask[i] + decayL[i]*mem;
+ mem = mask[i];
+ }
+ //for (i=0;i<len;i++) printf ("%f ", mask[i]); printf ("\n");
+}
+
+/* Compute a marking threshold from the spectrum X. */
+void compute_masking(float *X, float *mask, int len, int Fs)
+{
+ int i;
+ int N=len/2;
+ float psd[N];
+ psd[0] = X[0]*X[0];
+ for (i=1;i<N;i++)
+ psd[i] = X[i*2-1]*X[i*2-1] + X[i*2]*X[i*2];
+ /* TODO: Do tone masking */
+ /* Noise masking */
+ spreading_func(psd, mask, N, Fs);
+
+}
+
diff --git a/libcelt/psy.h b/libcelt/psy.h
new file mode 100644
index 0000000000000000000000000000000000000000..bf2c9ffb37b19ee032dae630fefe8f4db26612b3
--- /dev/null
+++ b/libcelt/psy.h
@@ -0,0 +1,36 @@
+/* (C) 2007 Jean-Marc Valin, CSIRO
+*/
+/*
+ Redistribution and use in source and binary forms, with or without
+ modification, are permitted provided that the following conditions
+ are met:
+
+ - Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+ - Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer in the
+ documentation and/or other materials provided with the distribution.
+
+ - Neither the name of the Xiph.org Foundation nor the names of its
+ contributors may be used to endorse or promote products derived from
+ this software without specific prior written permission.
+
+ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR
+ CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+ LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+#ifndef PSY_H
+#define PSY_H
+
+void compute_masking(float *X, float *mask, int len, int Fs);
+
+#endif /* PSY_H */