Directories
-----------
+script - shell scripts to simply common operations
+speex - LSP quantisation code borrowed from Speex for testing
src - C source code
octave - Matlab/Octave scripts
pitch - pitch estimator output files
--- /dev/null
+/* Copyright (C) 2003 Jean-Marc Valin */
+/**
+ @file arch.h
+ @brief Various architecture definitions Speex
+*/
+/*
+ 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 ARCH_H
+#define ARCH_H
+
+#ifndef SPEEX_VERSION
+#define SPEEX_MAJOR_VERSION 1 /**< Major Speex version. */
+#define SPEEX_MINOR_VERSION 1 /**< Minor Speex version. */
+#define SPEEX_MICRO_VERSION 15 /**< Micro Speex version. */
+#define SPEEX_EXTRA_VERSION "" /**< Extra Speex version. */
+#define SPEEX_VERSION "speex-1.2beta3" /**< Speex version string. */
+#endif
+
+/* A couple test to catch stupid option combinations */
+#ifdef FIXED_POINT
+
+#ifdef FLOATING_POINT
+#error You cannot compile as floating point and fixed point at the same time
+#endif
+#ifdef _USE_SSE
+#error SSE is only for floating-point
+#endif
+#if ((defined (ARM4_ASM)||defined (ARM4_ASM)) && defined(BFIN_ASM)) || (defined (ARM4_ASM)&&defined(ARM5E_ASM))
+#error Make up your mind. What CPU do you have?
+#endif
+#ifdef VORBIS_PSYCHO
+#error Vorbis-psy model currently not implemented in fixed-point
+#endif
+
+#else
+
+#ifndef FLOATING_POINT
+#error You now need to define either FIXED_POINT or FLOATING_POINT
+#endif
+#if defined (ARM4_ASM) || defined(ARM5E_ASM) || defined(BFIN_ASM)
+#error I suppose you can have a [ARM4/ARM5E/Blackfin] that has float instructions?
+#endif
+#ifdef FIXED_POINT_DEBUG
+#error "Don't you think enabling fixed-point is a good thing to do if you want to debug that?"
+#endif
+
+
+#endif
+
+#ifndef OUTSIDE_SPEEX
+#include "speex/speex_types.h"
+#endif
+
+#define ABS(x) ((x) < 0 ? (-(x)) : (x)) /**< Absolute integer value. */
+#define ABS16(x) ((x) < 0 ? (-(x)) : (x)) /**< Absolute 16-bit value. */
+#define MIN16(a,b) ((a) < (b) ? (a) : (b)) /**< Maximum 16-bit value. */
+#define MAX16(a,b) ((a) > (b) ? (a) : (b)) /**< Maximum 16-bit value. */
+#define ABS32(x) ((x) < 0 ? (-(x)) : (x)) /**< Absolute 32-bit value. */
+#define MIN32(a,b) ((a) < (b) ? (a) : (b)) /**< Maximum 32-bit value. */
+#define MAX32(a,b) ((a) > (b) ? (a) : (b)) /**< Maximum 32-bit value. */
+
+#ifdef FIXED_POINT
+
+typedef spx_int16_t spx_word16_t;
+typedef spx_int32_t spx_word32_t;
+typedef spx_word32_t spx_mem_t;
+typedef spx_word16_t spx_coef_t;
+typedef spx_word16_t spx_lsp_t;
+typedef spx_word32_t spx_sig_t;
+
+#define Q15ONE 32767
+
+#define LPC_SCALING 8192
+#define SIG_SCALING 16384
+#define LSP_SCALING 8192.
+#define GAMMA_SCALING 32768.
+#define GAIN_SCALING 64
+#define GAIN_SCALING_1 0.015625
+
+#define LPC_SHIFT 13
+#define LSP_SHIFT 13
+#define SIG_SHIFT 14
+#define GAIN_SHIFT 6
+
+#define VERY_SMALL 0
+#define VERY_LARGE32 ((spx_word32_t)2147483647)
+#define VERY_LARGE16 ((spx_word16_t)32767)
+#define Q15_ONE ((spx_word16_t)32767)
+
+
+#ifdef FIXED_DEBUG
+#include "fixed_debug.h"
+#else
+
+#include "fixed_generic.h"
+
+#ifdef ARM5E_ASM
+#include "fixed_arm5e.h"
+#elif defined (ARM4_ASM)
+#include "fixed_arm4.h"
+#elif defined (BFIN_ASM)
+#include "fixed_bfin.h"
+#endif
+
+#endif
+
+
+#else
+
+typedef float spx_mem_t;
+typedef float spx_coef_t;
+typedef float spx_lsp_t;
+typedef float spx_sig_t;
+typedef float spx_word16_t;
+typedef float spx_word32_t;
+
+#define Q15ONE 1.0f
+#define LPC_SCALING 1.f
+#define SIG_SCALING 1.f
+#define LSP_SCALING 1.f
+#define GAMMA_SCALING 1.f
+#define GAIN_SCALING 1.f
+#define GAIN_SCALING_1 1.f
+
+
+#define VERY_SMALL 1e-15f
+#define VERY_LARGE32 1e15f
+#define VERY_LARGE16 1e15f
+#define Q15_ONE ((spx_word16_t)1.f)
+
+#define QCONST16(x,bits) (x)
+#define QCONST32(x,bits) (x)
+
+#define NEG16(x) (-(x))
+#define NEG32(x) (-(x))
+#define EXTRACT16(x) (x)
+#define EXTEND32(x) (x)
+#define SHR16(a,shift) (a)
+#define SHL16(a,shift) (a)
+#define SHR32(a,shift) (a)
+#define SHL32(a,shift) (a)
+#define PSHR16(a,shift) (a)
+#define PSHR32(a,shift) (a)
+#define VSHR32(a,shift) (a)
+#define SATURATE16(x,a) (x)
+#define SATURATE32(x,a) (x)
+
+#define PSHR(a,shift) (a)
+#define SHR(a,shift) (a)
+#define SHL(a,shift) (a)
+#define SATURATE(x,a) (x)
+
+#define ADD16(a,b) ((a)+(b))
+#define SUB16(a,b) ((a)-(b))
+#define ADD32(a,b) ((a)+(b))
+#define SUB32(a,b) ((a)-(b))
+#define MULT16_16_16(a,b) ((a)*(b))
+#define MULT16_16(a,b) ((spx_word32_t)(a)*(spx_word32_t)(b))
+#define MAC16_16(c,a,b) ((c)+(spx_word32_t)(a)*(spx_word32_t)(b))
+
+#define MULT16_32_Q11(a,b) ((a)*(b))
+#define MULT16_32_Q13(a,b) ((a)*(b))
+#define MULT16_32_Q14(a,b) ((a)*(b))
+#define MULT16_32_Q15(a,b) ((a)*(b))
+#define MULT16_32_P15(a,b) ((a)*(b))
+
+#define MAC16_32_Q11(c,a,b) ((c)+(a)*(b))
+#define MAC16_32_Q15(c,a,b) ((c)+(a)*(b))
+
+#define MAC16_16_Q11(c,a,b) ((c)+(a)*(b))
+#define MAC16_16_Q13(c,a,b) ((c)+(a)*(b))
+#define MAC16_16_P13(c,a,b) ((c)+(a)*(b))
+#define MULT16_16_Q11_32(a,b) ((a)*(b))
+#define MULT16_16_Q13(a,b) ((a)*(b))
+#define MULT16_16_Q14(a,b) ((a)*(b))
+#define MULT16_16_Q15(a,b) ((a)*(b))
+#define MULT16_16_P15(a,b) ((a)*(b))
+#define MULT16_16_P13(a,b) ((a)*(b))
+#define MULT16_16_P14(a,b) ((a)*(b))
+
+#define DIV32_16(a,b) (((spx_word32_t)(a))/(spx_word16_t)(b))
+#define PDIV32_16(a,b) (((spx_word32_t)(a))/(spx_word16_t)(b))
+#define DIV32(a,b) (((spx_word32_t)(a))/(spx_word32_t)(b))
+#define PDIV32(a,b) (((spx_word32_t)(a))/(spx_word32_t)(b))
+
+
+#endif
+
+
+#if defined (CONFIG_TI_C54X) || defined (CONFIG_TI_C55X)
+
+/* 2 on TI C5x DSP */
+#define BYTES_PER_CHAR 2
+#define BITS_PER_CHAR 16
+#define LOG2_BITS_PER_CHAR 4
+
+#else
+
+#define BYTES_PER_CHAR 1
+#define BITS_PER_CHAR 8
+#define LOG2_BITS_PER_CHAR 3
+
+#endif
+
+
+
+#ifdef FIXED_DEBUG
+extern long long spx_mips;
+#endif
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2002 Jean-Marc Valin
+ File: speex_bits.c
+
+ Handles bit packing/unpacking
+
+ 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.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <speex/speex_bits.h>
+#include "arch.h"
+#include "os_support.h"
+
+/* Maximum size of the bit-stream (for fixed-size allocation) */
+#ifndef MAX_CHARS_PER_FRAME
+#define MAX_CHARS_PER_FRAME (2000/BYTES_PER_CHAR)
+#endif
+
+void speex_bits_init(SpeexBits *bits)
+{
+ bits->chars = (char*)speex_alloc(MAX_CHARS_PER_FRAME);
+ if (!bits->chars)
+ return;
+
+ bits->buf_size = MAX_CHARS_PER_FRAME;
+
+ bits->owner=1;
+
+ speex_bits_reset(bits);
+}
+
+void speex_bits_init_buffer(SpeexBits *bits, void *buff, int buf_size)
+{
+ bits->chars = (char*)buff;
+ bits->buf_size = buf_size;
+
+ bits->owner=0;
+
+ speex_bits_reset(bits);
+}
+
+void speex_bits_set_bit_buffer(SpeexBits *bits, void *buff, int buf_size)
+{
+ bits->chars = (char*)buff;
+ bits->buf_size = buf_size;
+
+ bits->owner=0;
+
+ bits->nbBits=buf_size<<LOG2_BITS_PER_CHAR;
+ bits->charPtr=0;
+ bits->bitPtr=0;
+ bits->overflow=0;
+
+}
+
+void speex_bits_destroy(SpeexBits *bits)
+{
+ if (bits->owner)
+ speex_free(bits->chars);
+ /* Will do something once the allocation is dynamic */
+}
+
+void speex_bits_reset(SpeexBits *bits)
+{
+ /* We only need to clear the first byte now */
+ bits->chars[0]=0;
+ bits->nbBits=0;
+ bits->charPtr=0;
+ bits->bitPtr=0;
+ bits->overflow=0;
+}
+
+void speex_bits_rewind(SpeexBits *bits)
+{
+ bits->charPtr=0;
+ bits->bitPtr=0;
+ bits->overflow=0;
+}
+
+void speex_bits_read_from(SpeexBits *bits, char *chars, int len)
+{
+ int i;
+ int nchars = len / BYTES_PER_CHAR;
+ if (nchars > bits->buf_size)
+ {
+ speex_notify("Packet is larger than allocated buffer");
+ if (bits->owner)
+ {
+ char *tmp = (char*)speex_realloc(bits->chars, nchars);
+ if (tmp)
+ {
+ bits->buf_size=nchars;
+ bits->chars=tmp;
+ } else {
+ nchars=bits->buf_size;
+ speex_warning("Could not resize input buffer: truncating input");
+ }
+ } else {
+ speex_warning("Do not own input buffer: truncating oversize input");
+ nchars=bits->buf_size;
+ }
+ }
+#if (BYTES_PER_CHAR==2)
+/* Swap bytes to proper endian order (could be done externally) */
+#define HTOLS(A) ((((A) >> 8)&0xff)|(((A) & 0xff)<<8))
+#else
+#define HTOLS(A) (A)
+#endif
+ for (i=0;i<nchars;i++)
+ bits->chars[i]=HTOLS(chars[i]);
+
+ bits->nbBits=nchars<<LOG2_BITS_PER_CHAR;
+ bits->charPtr=0;
+ bits->bitPtr=0;
+ bits->overflow=0;
+}
+
+static void speex_bits_flush(SpeexBits *bits)
+{
+ int nchars = ((bits->nbBits+BITS_PER_CHAR-1)>>LOG2_BITS_PER_CHAR);
+ if (bits->charPtr>0)
+ SPEEX_MOVE(bits->chars, &bits->chars[bits->charPtr], nchars-bits->charPtr);
+ bits->nbBits -= bits->charPtr<<LOG2_BITS_PER_CHAR;
+ bits->charPtr=0;
+}
+
+void speex_bits_read_whole_bytes(SpeexBits *bits, char *chars, int nbytes)
+{
+ int i,pos;
+ int nchars = nbytes/BYTES_PER_CHAR;
+
+ if (((bits->nbBits+BITS_PER_CHAR-1)>>LOG2_BITS_PER_CHAR)+nchars > bits->buf_size)
+ {
+ /* Packet is larger than allocated buffer */
+ if (bits->owner)
+ {
+ char *tmp = (char*)speex_realloc(bits->chars, (bits->nbBits>>LOG2_BITS_PER_CHAR)+nchars+1);
+ if (tmp)
+ {
+ bits->buf_size=(bits->nbBits>>LOG2_BITS_PER_CHAR)+nchars+1;
+ bits->chars=tmp;
+ } else {
+ nchars=bits->buf_size-(bits->nbBits>>LOG2_BITS_PER_CHAR)-1;
+ speex_warning("Could not resize input buffer: truncating oversize input");
+ }
+ } else {
+ speex_warning("Do not own input buffer: truncating oversize input");
+ nchars=bits->buf_size;
+ }
+ }
+
+ speex_bits_flush(bits);
+ pos=bits->nbBits>>LOG2_BITS_PER_CHAR;
+ for (i=0;i<nchars;i++)
+ bits->chars[pos+i]=HTOLS(chars[i]);
+ bits->nbBits+=nchars<<LOG2_BITS_PER_CHAR;
+}
+
+int speex_bits_write(SpeexBits *bits, char *chars, int max_nbytes)
+{
+ int i;
+ int max_nchars = max_nbytes/BYTES_PER_CHAR;
+ int charPtr, bitPtr, nbBits;
+
+ /* Insert terminator, but save the data so we can put it back after */
+ bitPtr=bits->bitPtr;
+ charPtr=bits->charPtr;
+ nbBits=bits->nbBits;
+ speex_bits_insert_terminator(bits);
+ bits->bitPtr=bitPtr;
+ bits->charPtr=charPtr;
+ bits->nbBits=nbBits;
+
+ if (max_nchars > ((bits->nbBits+BITS_PER_CHAR-1)>>LOG2_BITS_PER_CHAR))
+ max_nchars = ((bits->nbBits+BITS_PER_CHAR-1)>>LOG2_BITS_PER_CHAR);
+
+ for (i=0;i<max_nchars;i++)
+ chars[i]=HTOLS(bits->chars[i]);
+ return max_nchars*BYTES_PER_CHAR;
+}
+
+int speex_bits_write_whole_bytes(SpeexBits *bits, char *chars, int max_nbytes)
+{
+ int max_nchars = max_nbytes/BYTES_PER_CHAR;
+ int i;
+ if (max_nchars > ((bits->nbBits)>>LOG2_BITS_PER_CHAR))
+ max_nchars = ((bits->nbBits)>>LOG2_BITS_PER_CHAR);
+ for (i=0;i<max_nchars;i++)
+ chars[i]=HTOLS(bits->chars[i]);
+
+ if (bits->bitPtr>0)
+ bits->chars[0]=bits->chars[max_nchars];
+ else
+ bits->chars[0]=0;
+ bits->charPtr=0;
+ bits->nbBits &= (BITS_PER_CHAR-1);
+ return max_nchars*BYTES_PER_CHAR;
+}
+
+void speex_bits_pack(SpeexBits *bits, int data, int nbBits)
+{
+ unsigned int d=data;
+
+ if (bits->charPtr+((nbBits+bits->bitPtr)>>LOG2_BITS_PER_CHAR) >= bits->buf_size)
+ {
+ speex_notify("Buffer too small to pack bits");
+ if (bits->owner)
+ {
+ int new_nchars = ((bits->buf_size+5)*3)>>1;
+ char *tmp = (char*)speex_realloc(bits->chars, new_nchars);
+ if (tmp)
+ {
+ bits->buf_size=new_nchars;
+ bits->chars=tmp;
+ } else {
+ speex_warning("Could not resize input buffer: not packing");
+ return;
+ }
+ } else {
+ speex_warning("Do not own input buffer: not packing");
+ return;
+ }
+ }
+
+ while(nbBits)
+ {
+ int bit;
+ bit = (d>>(nbBits-1))&1;
+ bits->chars[bits->charPtr] |= bit<<(BITS_PER_CHAR-1-bits->bitPtr);
+ bits->bitPtr++;
+
+ if (bits->bitPtr==BITS_PER_CHAR)
+ {
+ bits->bitPtr=0;
+ bits->charPtr++;
+ bits->chars[bits->charPtr] = 0;
+ }
+ bits->nbBits++;
+ nbBits--;
+ }
+}
+
+int speex_bits_unpack_signed(SpeexBits *bits, int nbBits)
+{
+ unsigned int d=speex_bits_unpack_unsigned(bits,nbBits);
+ /* If number is negative */
+ if (d>>(nbBits-1))
+ {
+ d |= (-1)<<nbBits;
+ }
+ return d;
+}
+
+unsigned int speex_bits_unpack_unsigned(SpeexBits *bits, int nbBits)
+{
+ unsigned int d=0;
+ if ((bits->charPtr<<LOG2_BITS_PER_CHAR)+bits->bitPtr+nbBits>bits->nbBits)
+ bits->overflow=1;
+ if (bits->overflow)
+ return 0;
+ while(nbBits)
+ {
+ d<<=1;
+ d |= (bits->chars[bits->charPtr]>>(BITS_PER_CHAR-1 - bits->bitPtr))&1;
+ bits->bitPtr++;
+ if (bits->bitPtr==BITS_PER_CHAR)
+ {
+ bits->bitPtr=0;
+ bits->charPtr++;
+ }
+ nbBits--;
+ }
+ return d;
+}
+
+unsigned int speex_bits_peek_unsigned(SpeexBits *bits, int nbBits)
+{
+ unsigned int d=0;
+ int bitPtr, charPtr;
+ char *chars;
+
+ if ((bits->charPtr<<LOG2_BITS_PER_CHAR)+bits->bitPtr+nbBits>bits->nbBits)
+ bits->overflow=1;
+ if (bits->overflow)
+ return 0;
+
+ bitPtr=bits->bitPtr;
+ charPtr=bits->charPtr;
+ chars = bits->chars;
+ while(nbBits)
+ {
+ d<<=1;
+ d |= (chars[charPtr]>>(BITS_PER_CHAR-1 - bitPtr))&1;
+ bitPtr++;
+ if (bitPtr==BITS_PER_CHAR)
+ {
+ bitPtr=0;
+ charPtr++;
+ }
+ nbBits--;
+ }
+ return d;
+}
+
+int speex_bits_peek(SpeexBits *bits)
+{
+ if ((bits->charPtr<<LOG2_BITS_PER_CHAR)+bits->bitPtr+1>bits->nbBits)
+ bits->overflow=1;
+ if (bits->overflow)
+ return 0;
+ return (bits->chars[bits->charPtr]>>(BITS_PER_CHAR-1 - bits->bitPtr))&1;
+}
+
+void speex_bits_advance(SpeexBits *bits, int n)
+{
+ if (((bits->charPtr<<LOG2_BITS_PER_CHAR)+bits->bitPtr+n>bits->nbBits) || bits->overflow){
+ bits->overflow=1;
+ return;
+ }
+ bits->charPtr += (bits->bitPtr+n) >> LOG2_BITS_PER_CHAR; /* divide by BITS_PER_CHAR */
+ bits->bitPtr = (bits->bitPtr+n) & (BITS_PER_CHAR-1); /* modulo by BITS_PER_CHAR */
+}
+
+int speex_bits_remaining(SpeexBits *bits)
+{
+ if (bits->overflow)
+ return -1;
+ else
+ return bits->nbBits-((bits->charPtr<<LOG2_BITS_PER_CHAR)+bits->bitPtr);
+}
+
+int speex_bits_nbytes(SpeexBits *bits)
+{
+ return ((bits->nbBits+BITS_PER_CHAR-1)>>LOG2_BITS_PER_CHAR);
+}
+
+void speex_bits_insert_terminator(SpeexBits *bits)
+{
+ if (bits->bitPtr)
+ speex_bits_pack(bits, 0, 1);
+ while (bits->bitPtr)
+ speex_bits_pack(bits, 1, 1);
+}
--- /dev/null
+/* Copyright (C) 2003 Jean-Marc Valin */
+/**
+ @file fixed_generic.h
+ @brief Generic fixed-point operations
+*/
+/*
+ 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 FIXED_GENERIC_H
+#define FIXED_GENERIC_H
+
+#define QCONST16(x,bits) ((spx_word16_t)(.5+(x)*(((spx_word32_t)1)<<(bits))))
+#define QCONST32(x,bits) ((spx_word32_t)(.5+(x)*(((spx_word32_t)1)<<(bits))))
+
+#define NEG16(x) (-(x))
+#define NEG32(x) (-(x))
+#define EXTRACT16(x) ((spx_word16_t)(x))
+#define EXTEND32(x) ((spx_word32_t)(x))
+#define SHR16(a,shift) ((a) >> (shift))
+#define SHL16(a,shift) ((a) << (shift))
+#define SHR32(a,shift) ((a) >> (shift))
+#define SHL32(a,shift) ((a) << (shift))
+#define PSHR16(a,shift) (SHR16((a)+((1<<((shift))>>1)),shift))
+#define PSHR32(a,shift) (SHR32((a)+((EXTEND32(1)<<((shift))>>1)),shift))
+#define VSHR32(a, shift) (((shift)>0) ? SHR32(a, shift) : SHL32(a, -(shift)))
+#define SATURATE16(x,a) (((x)>(a) ? (a) : (x)<-(a) ? -(a) : (x)))
+#define SATURATE32(x,a) (((x)>(a) ? (a) : (x)<-(a) ? -(a) : (x)))
+
+#define SHR(a,shift) ((a) >> (shift))
+#define SHL(a,shift) ((spx_word32_t)(a) << (shift))
+#define PSHR(a,shift) (SHR((a)+((EXTEND32(1)<<((shift))>>1)),shift))
+#define SATURATE(x,a) (((x)>(a) ? (a) : (x)<-(a) ? -(a) : (x)))
+
+
+#define ADD16(a,b) ((spx_word16_t)((spx_word16_t)(a)+(spx_word16_t)(b)))
+#define SUB16(a,b) ((spx_word16_t)(a)-(spx_word16_t)(b))
+#define ADD32(a,b) ((spx_word32_t)(a)+(spx_word32_t)(b))
+#define SUB32(a,b) ((spx_word32_t)(a)-(spx_word32_t)(b))
+
+
+/* result fits in 16 bits */
+#define MULT16_16_16(a,b) ((((spx_word16_t)(a))*((spx_word16_t)(b))))
+
+/* (spx_word32_t)(spx_word16_t) gives TI compiler a hint that it's 16x16->32 multiply */
+#define MULT16_16(a,b) (((spx_word32_t)(spx_word16_t)(a))*((spx_word32_t)(spx_word16_t)(b)))
+
+#define MAC16_16(c,a,b) (ADD32((c),MULT16_16((a),(b))))
+#define MULT16_32_Q12(a,b) ADD32(MULT16_16((a),SHR((b),12)), SHR(MULT16_16((a),((b)&0x00000fff)),12))
+#define MULT16_32_Q13(a,b) ADD32(MULT16_16((a),SHR((b),13)), SHR(MULT16_16((a),((b)&0x00001fff)),13))
+#define MULT16_32_Q14(a,b) ADD32(MULT16_16((a),SHR((b),14)), SHR(MULT16_16((a),((b)&0x00003fff)),14))
+
+#define MULT16_32_Q11(a,b) ADD32(MULT16_16((a),SHR((b),11)), SHR(MULT16_16((a),((b)&0x000007ff)),11))
+#define MAC16_32_Q11(c,a,b) ADD32(c,ADD32(MULT16_16((a),SHR((b),11)), SHR(MULT16_16((a),((b)&0x000007ff)),11)))
+
+#define MULT16_32_P15(a,b) ADD32(MULT16_16((a),SHR((b),15)), PSHR(MULT16_16((a),((b)&0x00007fff)),15))
+#define MULT16_32_Q15(a,b) ADD32(MULT16_16((a),SHR((b),15)), SHR(MULT16_16((a),((b)&0x00007fff)),15))
+#define MAC16_32_Q15(c,a,b) ADD32(c,ADD32(MULT16_16((a),SHR((b),15)), SHR(MULT16_16((a),((b)&0x00007fff)),15)))
+
+
+#define MAC16_16_Q11(c,a,b) (ADD32((c),SHR(MULT16_16((a),(b)),11)))
+#define MAC16_16_Q13(c,a,b) (ADD32((c),SHR(MULT16_16((a),(b)),13)))
+#define MAC16_16_P13(c,a,b) (ADD32((c),SHR(ADD32(4096,MULT16_16((a),(b))),13)))
+
+#define MULT16_16_Q11_32(a,b) (SHR(MULT16_16((a),(b)),11))
+#define MULT16_16_Q13(a,b) (SHR(MULT16_16((a),(b)),13))
+#define MULT16_16_Q14(a,b) (SHR(MULT16_16((a),(b)),14))
+#define MULT16_16_Q15(a,b) (SHR(MULT16_16((a),(b)),15))
+
+#define MULT16_16_P13(a,b) (SHR(ADD32(4096,MULT16_16((a),(b))),13))
+#define MULT16_16_P14(a,b) (SHR(ADD32(8192,MULT16_16((a),(b))),14))
+#define MULT16_16_P15(a,b) (SHR(ADD32(16384,MULT16_16((a),(b))),15))
+
+#define MUL_16_32_R15(a,bh,bl) ADD32(MULT16_16((a),(bh)), SHR(MULT16_16((a),(bl)),15))
+
+#define DIV32_16(a,b) ((spx_word16_t)(((spx_word32_t)(a))/((spx_word16_t)(b))))
+#define PDIV32_16(a,b) ((spx_word16_t)(((spx_word32_t)(a)+((spx_word16_t)(b)>>1))/((spx_word16_t)(b))))
+#define DIV32(a,b) (((spx_word32_t)(a))/((spx_word32_t)(b)))
+#define PDIV32(a,b) (((spx_word32_t)(a)+((spx_word16_t)(b)>>1))/((spx_word32_t)(b)))
+
+#endif
--- /dev/null
+/* Copyright (C) 2002 Jean-Marc Valin
+ File: high_lsp_tables.c
+ Codebooks for high-band LSPs in SB-CELP mode
+
+ Redistribution and use in source and binary forms, with or without
+ modification, are permitted provided that the following conditions are
+ met:
+
+ 1. Redistributions of source code must retain the above copyright notice,
+ this list of conditions and the following disclaimer.
+
+ 2. 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.
+
+ 3. The name of the author may not be used to endorse or promote products
+ derived from this software without specific prior written permission.
+
+ THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``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 AUTHOR 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.
+*/
+
+const signed char high_lsp_cdbk[512]={
+39,12,-14,-20,-29,-61,-67,-76,
+-32,-71,-67,68,77,46,34,5,
+-13,-48,-46,-72,-81,-84,-60,-58,
+-40,-28,82,93,68,45,29,3,
+-19,-47,-28,-43,-35,-30,-8,-13,
+-39,-91,-91,-123,-96,10,10,-6,
+-18,-55,-60,-91,-56,-36,-27,-16,
+-48,-75,40,28,-10,-28,35,9,
+37,19,1,-20,-31,-41,-18,-25,
+-35,-68,-80,45,27,-1,47,13,
+0,-29,-35,-57,-50,-79,-73,-38,
+-19,5,35,14,-10,-23,16,-8,
+5,-24,-40,-62,-23,-27,-22,-16,
+-18,-46,-72,-77,43,21,33,1,
+-80,-70,-70,-64,-56,-52,-39,-33,
+-31,-38,-19,-19,-15,32,33,-2,
+7,-15,-15,-24,-23,-33,-41,-56,
+-24,-57,5,89,64,41,27,5,
+-9,-47,-60,-97,-97,-124,-20,-9,
+-44,-73,31,29,-4,64,48,7,
+-35,-57,0,-3,-26,-47,-3,-6,
+-40,-76,-79,-48,12,81,55,10,
+9,-24,-43,-73,-57,-69,16,5,
+-28,-53,18,29,20,0,-4,-11,
+6,-13,23,7,-17,-35,-37,-37,
+-30,-68,-63,6,24,-9,-14,3,
+21,-13,-27,-57,-49,-80,-24,-41,
+-5,-16,-5,1,45,25,12,-7,
+3,-15,-6,-16,-15,-8,6,-13,
+-42,-81,-80,-87,14,1,-10,-3,
+-43,-69,-46,-24,-28,-29,36,6,
+-43,-56,-12,12,54,79,43,9,
+54,22,2,8,-12,-43,-46,-52,
+-38,-69,-89,-5,75,38,33,5,
+-13,-53,-62,-87,-89,-113,-99,-55,
+-34,-37,62,55,33,16,21,-2,
+-17,-46,-29,-38,-38,-48,-39,-42,
+-36,-75,-72,-88,-48,-30,21,2,
+-15,-57,-64,-98,-84,-76,25,1,
+-46,-80,-12,18,-7,3,34,6,
+38,31,23,4,-1,20,14,-15,
+-43,-78,-91,-24,14,-3,54,16,
+0,-27,-28,-44,-56,-83,-92,-89,
+-3,34,56,41,36,22,20,-8,
+-7,-35,-42,-62,-49,3,12,-10,
+-50,-87,-96,-66,92,70,38,9,
+-70,-71,-62,-42,-39,-43,-11,-7,
+-50,-79,-58,-50,-31,32,31,-6,
+-4,-25,7,-17,-38,-70,-58,-27,
+-43,-83,-28,59,36,20,31,2,
+-27,-71,-80,-109,-98,-75,-33,-32,
+-31,-2,33,15,-6,43,33,-5,
+0,-22,-10,-27,-34,-49,-11,-20,
+-41,-91,-100,-121,-39,57,41,10,
+-19,-50,-38,-59,-60,-70,-18,-20,
+-8,-31,-8,-15,1,-14,-26,-25,
+33,21,32,17,1,-19,-19,-26,
+-58,-81,-35,-22,45,30,11,-11,
+3,-26,-48,-87,-67,-83,-58,3,
+-1,-26,-20,44,10,25,39,5,
+-9,-35,-27,-38,7,10,4,-9,
+-42,-85,-102,-127,52,44,28,10,
+-47,-61,-40,-39,-17,-1,-10,-33,
+-42,-74,-48,21,-4,70,52,10};
+
+
+const signed char high_lsp_cdbk2[512]={
+-36,-62,6,-9,-10,-14,-56,23,
+1,-26,23,-48,-17,12,8,-7,
+23,29,-36,-28,-6,-29,-17,-5,
+40,23,10,10,-46,-13,36,6,
+4,-30,-29,62,32,-32,-1,22,
+-14,1,-4,-22,-45,2,54,4,
+-30,-57,-59,-12,27,-3,-31,8,
+-9,5,10,-14,32,66,19,9,
+2,-25,-37,23,-15,18,-38,-31,
+5,-9,-21,15,0,22,62,30,
+15,-12,-14,-46,77,21,33,3,
+34,29,-19,50,2,11,9,-38,
+-12,-37,62,1,-15,54,32,6,
+2,-24,20,35,-21,2,19,24,
+-13,55,4,9,39,-19,30,-1,
+-21,73,54,33,8,18,3,15,
+6,-19,-47,6,-3,-48,-50,1,
+26,20,8,-23,-50,65,-14,-55,
+-17,-31,-37,-28,53,-1,-17,-53,
+1,57,11,-8,-25,-30,-37,64,
+5,-52,-45,15,23,31,15,14,
+-25,24,33,-2,-44,-56,-18,6,
+-21,-43,4,-12,17,-37,20,-10,
+34,15,2,15,55,21,-11,-31,
+-6,46,25,16,-9,-25,-8,-62,
+28,17,20,-32,-29,26,30,25,
+-19,2,-16,-17,26,-51,2,50,
+42,19,-66,23,29,-2,3,19,
+-19,-37,32,15,6,30,-34,13,
+11,-5,40,31,10,-42,4,-9,
+26,-9,-70,17,-2,-23,20,-22,
+-55,51,-24,-31,22,-22,15,-13,
+3,-10,-28,-16,56,4,-63,11,
+-18,-15,-18,-38,-35,16,-7,34,
+-1,-21,-49,-47,9,-37,7,8,
+69,55,20,6,-33,-45,-10,-9,
+6,-9,12,71,15,-3,-42,-7,
+-24,32,-35,-2,-42,-17,-5,0,
+-2,-33,-54,13,-12,-34,47,23,
+19,55,7,-8,74,31,14,16,
+-23,-26,19,12,-18,-49,-28,-31,
+-20,2,-14,-20,-47,78,40,13,
+-23,-11,21,-6,18,1,47,5,
+38,35,32,46,22,8,13,16,
+-14,18,51,19,40,39,11,-26,
+-1,-17,47,2,-53,-15,31,-22,
+38,21,-15,-16,5,-33,53,15,
+-38,86,11,-3,-24,49,13,-4,
+-11,-18,28,20,-12,-27,-26,35,
+-25,-35,-3,-20,-61,30,10,-55,
+-12,-22,-52,-54,-14,19,-32,-12,
+45,15,-8,-48,-9,11,-32,8,
+-16,-34,-13,51,18,38,-2,-32,
+-17,22,-2,-18,-28,-70,59,27,
+-28,-19,-10,-20,-9,-9,-8,-21,
+21,-8,35,-2,45,-3,-9,12,
+0,30,7,-39,43,27,-38,-91,
+30,26,19,-55,-4,63,14,-17,
+13,9,13,2,7,4,6,61,
+72,-1,-17,29,-1,-22,-17,8,
+-28,-37,63,44,41,3,2,14,
+9,-6,75,-8,-7,-12,-15,-12,
+13,9,-4,30,-22,-65,15,0,
+-45,4,-4,1,5,22,11,23};
--- /dev/null
+/*---------------------------------------------------------------------------*\
+Original copyright
+ FILE........: lsp.c
+ AUTHOR......: David Rowe
+ DATE CREATED: 24/2/93
+
+Heavily modified by Jean-Marc Valin (c) 2002-2006 (fixed-point,
+ optimizations, additional functions, ...)
+
+ This file contains functions for converting Linear Prediction
+ Coefficients (LPC) to Line Spectral Pair (LSP) and back. Note that the
+ LSP coefficients are not in radians format but in the x domain of the
+ unit circle.
+
+ Speex License:
+
+ 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.
+*/
+
+/*---------------------------------------------------------------------------*\
+
+ Introduction to Line Spectrum Pairs (LSPs)
+ ------------------------------------------
+
+ LSPs are used to encode the LPC filter coefficients {ak} for
+ transmission over the channel. LSPs have several properties (like
+ less sensitivity to quantisation noise) that make them superior to
+ direct quantisation of {ak}.
+
+ A(z) is a polynomial of order lpcrdr with {ak} as the coefficients.
+
+ A(z) is transformed to P(z) and Q(z) (using a substitution and some
+ algebra), to obtain something like:
+
+ A(z) = 0.5[P(z)(z+z^-1) + Q(z)(z-z^-1)] (1)
+
+ As you can imagine A(z) has complex zeros all over the z-plane. P(z)
+ and Q(z) have the very neat property of only having zeros _on_ the
+ unit circle. So to find them we take a test point z=exp(jw) and
+ evaluate P (exp(jw)) and Q(exp(jw)) using a grid of points between 0
+ and pi.
+
+ The zeros (roots) of P(z) also happen to alternate, which is why we
+ swap coefficients as we find roots. So the process of finding the
+ LSP frequencies is basically finding the roots of 5th order
+ polynomials.
+
+ The root so P(z) and Q(z) occur in symmetrical pairs at +/-w, hence
+ the name Line Spectrum Pairs (LSPs).
+
+ To convert back to ak we just evaluate (1), "clocking" an impulse
+ thru it lpcrdr times gives us the impulse response of A(z) which is
+ {ak}.
+
+\*---------------------------------------------------------------------------*/
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <math.h>
+#include "lsp.h"
+#include "stack_alloc.h"
+#include "math_approx.h"
+
+#ifndef M_PI
+#define M_PI 3.14159265358979323846 /* pi */
+#endif
+
+#ifndef NULL
+#define NULL 0
+#endif
+
+#ifdef FIXED_POINT
+
+#define FREQ_SCALE 16384
+
+/*#define ANGLE2X(a) (32768*cos(((a)/8192.)))*/
+#define ANGLE2X(a) (SHL16(spx_cos(a),2))
+
+/*#define X2ANGLE(x) (acos(.00006103515625*(x))*LSP_SCALING)*/
+#define X2ANGLE(x) (spx_acos(x))
+
+#ifdef BFIN_ASM
+#include "lsp_bfin.h"
+#endif
+
+#else
+
+/*#define C1 0.99940307
+#define C2 -0.49558072
+#define C3 0.03679168*/
+
+#define FREQ_SCALE 1.
+#define ANGLE2X(a) (spx_cos(a))
+#define X2ANGLE(x) (acos(x))
+
+#endif
+
+
+/*---------------------------------------------------------------------------*\
+
+ FUNCTION....: cheb_poly_eva()
+
+ AUTHOR......: David Rowe
+ DATE CREATED: 24/2/93
+
+ This function evaluates a series of Chebyshev polynomials
+
+\*---------------------------------------------------------------------------*/
+
+#ifdef FIXED_POINT
+
+#ifndef OVERRIDE_CHEB_POLY_EVA
+static inline spx_word32_t cheb_poly_eva(
+ spx_word16_t *coef, /* P or Q coefs in Q13 format */
+ spx_word16_t x, /* cos of freq (-1.0 to 1.0) in Q14 format */
+ int m, /* LPC order/2 */
+ char *stack
+)
+{
+ int i;
+ spx_word16_t b0, b1;
+ spx_word32_t sum;
+
+ /*Prevents overflows*/
+ if (x>16383)
+ x = 16383;
+ if (x<-16383)
+ x = -16383;
+
+ /* Initialise values */
+ b1=16384;
+ b0=x;
+
+ /* Evaluate Chebyshev series formulation usin g iterative approach */
+ sum = ADD32(EXTEND32(coef[m]), EXTEND32(MULT16_16_P14(coef[m-1],x)));
+ for(i=2;i<=m;i++)
+ {
+ spx_word16_t tmp=b0;
+ b0 = SUB16(MULT16_16_Q13(x,b0), b1);
+ b1 = tmp;
+ sum = ADD32(sum, EXTEND32(MULT16_16_P14(coef[m-i],b0)));
+ }
+
+ return sum;
+}
+#endif
+
+#else
+
+static float cheb_poly_eva(spx_word32_t *coef, spx_word16_t x, int m, char *stack)
+{
+ int k;
+ float b0, b1, tmp;
+
+ /* Initial conditions */
+ b0=0; /* b_(m+1) */
+ b1=0; /* b_(m+2) */
+
+ x*=2;
+
+ /* Calculate the b_(k) */
+ for(k=m;k>0;k--)
+ {
+ tmp=b0; /* tmp holds the previous value of b0 */
+ b0=x*b0-b1+coef[m-k]; /* b0 holds its new value based on b0 and b1 */
+ b1=tmp; /* b1 holds the previous value of b0 */
+ }
+
+ return(-b1+.5*x*b0+coef[m]);
+}
+#endif
+
+/*---------------------------------------------------------------------------*\
+
+ FUNCTION....: lpc_to_lsp()
+
+ AUTHOR......: David Rowe
+ DATE CREATED: 24/2/93
+
+ This function converts LPC coefficients to LSP
+ coefficients.
+
+\*---------------------------------------------------------------------------*/
+
+#ifdef FIXED_POINT
+#define SIGN_CHANGE(a,b) (((a)&0x70000000)^((b)&0x70000000)||(b==0))
+#else
+#define SIGN_CHANGE(a,b) (((a)*(b))<0.0)
+#endif
+
+
+int lpc_to_lsp (spx_coef_t *a,int lpcrdr,spx_lsp_t *freq,int nb,spx_word16_t delta, char *stack)
+/* float *a lpc coefficients */
+/* int lpcrdr order of LPC coefficients (10) */
+/* float *freq LSP frequencies in the x domain */
+/* int nb number of sub-intervals (4) */
+/* float delta grid spacing interval (0.02) */
+
+
+{
+ spx_word16_t temp_xr,xl,xr,xm=0;
+ spx_word32_t psuml,psumr,psumm,temp_psumr/*,temp_qsumr*/;
+ int i,j,m,flag,k;
+ VARDECL(spx_word32_t *Q); /* ptrs for memory allocation */
+ VARDECL(spx_word32_t *P);
+ VARDECL(spx_word16_t *Q16); /* ptrs for memory allocation */
+ VARDECL(spx_word16_t *P16);
+ spx_word32_t *px; /* ptrs of respective P'(z) & Q'(z) */
+ spx_word32_t *qx;
+ spx_word32_t *p;
+ spx_word32_t *q;
+ spx_word16_t *pt; /* ptr used for cheb_poly_eval()
+ whether P' or Q' */
+ int roots=0; /* DR 8/2/94: number of roots found */
+ flag = 1; /* program is searching for a root when,
+ 1 else has found one */
+ m = lpcrdr/2; /* order of P'(z) & Q'(z) polynomials */
+
+ /* Allocate memory space for polynomials */
+ ALLOC(Q, (m+1), spx_word32_t);
+ ALLOC(P, (m+1), spx_word32_t);
+
+ /* determine P'(z)'s and Q'(z)'s coefficients where
+ P'(z) = P(z)/(1 + z^(-1)) and Q'(z) = Q(z)/(1-z^(-1)) */
+
+ px = P; /* initialise ptrs */
+ qx = Q;
+ p = px;
+ q = qx;
+
+#ifdef FIXED_POINT
+ *px++ = LPC_SCALING;
+ *qx++ = LPC_SCALING;
+ for(i=0;i<m;i++){
+ *px++ = SUB32(ADD32(EXTEND32(a[i]),EXTEND32(a[lpcrdr-i-1])), *p++);
+ *qx++ = ADD32(SUB32(EXTEND32(a[i]),EXTEND32(a[lpcrdr-i-1])), *q++);
+ }
+ px = P;
+ qx = Q;
+ for(i=0;i<m;i++)
+ {
+ /*if (fabs(*px)>=32768)
+ speex_warning_int("px", *px);
+ if (fabs(*qx)>=32768)
+ speex_warning_int("qx", *qx);*/
+ *px = PSHR32(*px,2);
+ *qx = PSHR32(*qx,2);
+ px++;
+ qx++;
+ }
+ /* The reason for this lies in the way cheb_poly_eva() is implemented for fixed-point */
+ P[m] = PSHR32(P[m],3);
+ Q[m] = PSHR32(Q[m],3);
+#else
+ *px++ = LPC_SCALING;
+ *qx++ = LPC_SCALING;
+ for(i=0;i<m;i++){
+ *px++ = (a[i]+a[lpcrdr-1-i]) - *p++;
+ *qx++ = (a[i]-a[lpcrdr-1-i]) + *q++;
+ }
+ px = P;
+ qx = Q;
+ for(i=0;i<m;i++){
+ *px = 2**px;
+ *qx = 2**qx;
+ px++;
+ qx++;
+ }
+#endif
+
+ px = P; /* re-initialise ptrs */
+ qx = Q;
+
+ /* now that we have computed P and Q convert to 16 bits to
+ speed up cheb_poly_eval */
+
+ ALLOC(P16, m+1, spx_word16_t);
+ ALLOC(Q16, m+1, spx_word16_t);
+
+ for (i=0;i<m+1;i++)
+ {
+ P16[i] = P[i];
+ Q16[i] = Q[i];
+ }
+
+ /* Search for a zero in P'(z) polynomial first and then alternate to Q'(z).
+ Keep alternating between the two polynomials as each zero is found */
+
+ xr = 0; /* initialise xr to zero */
+ xl = FREQ_SCALE; /* start at point xl = 1 */
+
+ for(j=0;j<lpcrdr;j++){
+ if(j&1) /* determines whether P' or Q' is eval. */
+ pt = Q16;
+ else
+ pt = P16;
+
+ psuml = cheb_poly_eva(pt,xl,m,stack); /* evals poly. at xl */
+ flag = 1;
+ while(flag && (xr >= -FREQ_SCALE)){
+ spx_word16_t dd;
+ /* Modified by JMV to provide smaller steps around x=+-1 */
+#ifdef FIXED_POINT
+ dd = MULT16_16_Q15(delta,SUB16(FREQ_SCALE, MULT16_16_Q14(MULT16_16_Q14(xl,xl),14000)));
+ if (psuml<512 && psuml>-512)
+ dd = PSHR16(dd,1);
+#else
+ dd=delta*(1-.9*xl*xl);
+ if (fabs(psuml)<.2)
+ dd *= .5;
+#endif
+ xr = SUB16(xl, dd); /* interval spacing */
+ psumr = cheb_poly_eva(pt,xr,m,stack);/* poly(xl-delta_x) */
+ temp_psumr = psumr;
+ temp_xr = xr;
+
+ /* if no sign change increment xr and re-evaluate poly(xr). Repeat til
+ sign change.
+ if a sign change has occurred the interval is bisected and then
+ checked again for a sign change which determines in which
+ interval the zero lies in.
+ If there is no sign change between poly(xm) and poly(xl) set interval
+ between xm and xr else set interval between xl and xr and repeat till
+ root is located within the specified limits */
+
+ if(SIGN_CHANGE(psumr,psuml))
+ {
+ roots++;
+
+ psumm=psuml;
+ for(k=0;k<=nb;k++){
+#ifdef FIXED_POINT
+ xm = ADD16(PSHR16(xl,1),PSHR16(xr,1)); /* bisect the interval */
+#else
+ xm = .5*(xl+xr); /* bisect the interval */
+#endif
+ psumm=cheb_poly_eva(pt,xm,m,stack);
+ /*if(psumm*psuml>0.)*/
+ if(!SIGN_CHANGE(psumm,psuml))
+ {
+ psuml=psumm;
+ xl=xm;
+ } else {
+ psumr=psumm;
+ xr=xm;
+ }
+ }
+
+ /* once zero is found, reset initial interval to xr */
+ freq[j] = X2ANGLE(xm);
+ xl = xm;
+ flag = 0; /* reset flag for next search */
+ }
+ else{
+ psuml=temp_psumr;
+ xl=temp_xr;
+ }
+ }
+ }
+ return(roots);
+}
+
+/*---------------------------------------------------------------------------*\
+
+ FUNCTION....: lsp_to_lpc()
+
+ AUTHOR......: David Rowe
+ DATE CREATED: 24/2/93
+
+ Converts LSP coefficients to LPC coefficients.
+
+\*---------------------------------------------------------------------------*/
+
+#ifdef FIXED_POINT
+
+void lsp_to_lpc(spx_lsp_t *freq,spx_coef_t *ak,int lpcrdr, char *stack)
+/* float *freq array of LSP frequencies in the x domain */
+/* float *ak array of LPC coefficients */
+/* int lpcrdr order of LPC coefficients */
+{
+ int i,j;
+ spx_word32_t xout1,xout2,xin;
+ spx_word32_t mult, a;
+ VARDECL(spx_word16_t *freqn);
+ VARDECL(spx_word32_t **xp);
+ VARDECL(spx_word32_t *xpmem);
+ VARDECL(spx_word32_t **xq);
+ VARDECL(spx_word32_t *xqmem);
+ int m = lpcrdr>>1;
+
+ /*
+
+ Reconstruct P(z) and Q(z) by cascading second order polynomials
+ in form 1 - 2cos(w)z(-1) + z(-2), where w is the LSP frequency.
+ In the time domain this is:
+
+ y(n) = x(n) - 2cos(w)x(n-1) + x(n-2)
+
+ This is what the ALLOCS below are trying to do:
+
+ int xp[m+1][lpcrdr+1+2]; // P matrix in QIMP
+ int xq[m+1][lpcrdr+1+2]; // Q matrix in QIMP
+
+ These matrices store the output of each stage on each row. The
+ final (m-th) row has the output of the final (m-th) cascaded
+ 2nd order filter. The first row is the impulse input to the
+ system (not written as it is known).
+
+ The version below takes advantage of the fact that a lot of the
+ outputs are zero or known, for example if we put an inpulse
+ into the first section the "clock" it 10 times only the first 3
+ outputs samples are non-zero (it's an FIR filter).
+ */
+
+ ALLOC(xp, (m+1), spx_word32_t*);
+ ALLOC(xpmem, (m+1)*(lpcrdr+1+2), spx_word32_t);
+
+ ALLOC(xq, (m+1), spx_word32_t*);
+ ALLOC(xqmem, (m+1)*(lpcrdr+1+2), spx_word32_t);
+
+ for(i=0; i<=m; i++) {
+ xp[i] = xpmem + i*(lpcrdr+1+2);
+ xq[i] = xqmem + i*(lpcrdr+1+2);
+ }
+
+ /* work out 2cos terms in Q14 */
+
+ ALLOC(freqn, lpcrdr, spx_word16_t);
+ for (i=0;i<lpcrdr;i++)
+ freqn[i] = ANGLE2X(freq[i]);
+
+ #define QIMP 21 /* scaling for impulse */
+
+ xin = SHL32(EXTEND32(1), (QIMP-1)); /* 0.5 in QIMP format */
+
+ /* first col and last non-zero values of each row are trivial */
+
+ for(i=0;i<=m;i++) {
+ xp[i][1] = 0;
+ xp[i][2] = xin;
+ xp[i][2+2*i] = xin;
+ xq[i][1] = 0;
+ xq[i][2] = xin;
+ xq[i][2+2*i] = xin;
+ }
+
+ /* 2nd row (first output row) is trivial */
+
+ xp[1][3] = -MULT16_32_Q14(freqn[0],xp[0][2]);
+ xq[1][3] = -MULT16_32_Q14(freqn[1],xq[0][2]);
+
+ xout1 = xout2 = 0;
+
+ /* now generate remaining rows */
+
+ for(i=1;i<m;i++) {
+
+ for(j=1;j<2*(i+1)-1;j++) {
+ mult = MULT16_32_Q14(freqn[2*i],xp[i][j+1]);
+ xp[i+1][j+2] = ADD32(SUB32(xp[i][j+2], mult), xp[i][j]);
+ mult = MULT16_32_Q14(freqn[2*i+1],xq[i][j+1]);
+ xq[i+1][j+2] = ADD32(SUB32(xq[i][j+2], mult), xq[i][j]);
+ }
+
+ /* for last col xp[i][j+2] = xq[i][j+2] = 0 */
+
+ mult = MULT16_32_Q14(freqn[2*i],xp[i][j+1]);
+ xp[i+1][j+2] = SUB32(xp[i][j], mult);
+ mult = MULT16_32_Q14(freqn[2*i+1],xq[i][j+1]);
+ xq[i+1][j+2] = SUB32(xq[i][j], mult);
+ }
+
+ /* process last row to extra a{k} */
+
+ for(j=1;j<=lpcrdr;j++) {
+ int shift = QIMP-13;
+
+ /* final filter sections */
+ a = PSHR32(xp[m][j+2] + xout1 + xq[m][j+2] - xout2, shift);
+ xout1 = xp[m][j+2];
+ xout2 = xq[m][j+2];
+
+ /* hard limit ak's to +/- 32767 */
+
+ if (a < -32767) a = -32767;
+ if (a > 32767) a = 32767;
+ ak[j-1] = (short)a;
+
+ }
+
+}
+
+#else
+
+void lsp_to_lpc(spx_lsp_t *freq,spx_coef_t *ak,int lpcrdr, char *stack)
+/* float *freq array of LSP frequencies in the x domain */
+/* float *ak array of LPC coefficients */
+/* int lpcrdr order of LPC coefficients */
+
+
+{
+ int i,j;
+ float xout1,xout2,xin1,xin2;
+ VARDECL(float *Wp);
+ float *pw,*n1,*n2,*n3,*n4=NULL;
+ VARDECL(float *x_freq);
+ int m = lpcrdr>>1;
+
+ ALLOC(Wp, 4*m+2, float);
+ pw = Wp;
+
+ /* initialise contents of array */
+
+ for(i=0;i<=4*m+1;i++){ /* set contents of buffer to 0 */
+ *pw++ = 0.0;
+ }
+
+ /* Set pointers up */
+
+ pw = Wp;
+ xin1 = 1.0;
+ xin2 = 1.0;
+
+ ALLOC(x_freq, lpcrdr, float);
+ for (i=0;i<lpcrdr;i++)
+ x_freq[i] = ANGLE2X(freq[i]);
+
+ /* reconstruct P(z) and Q(z) by cascading second order
+ polynomials in form 1 - 2xz(-1) +z(-2), where x is the
+ LSP coefficient */
+
+ for(j=0;j<=lpcrdr;j++){
+ int i2=0;
+ for(i=0;i<m;i++,i2+=2){
+ n1 = pw+(i*4);
+ n2 = n1 + 1;
+ n3 = n2 + 1;
+ n4 = n3 + 1;
+ xout1 = xin1 - 2.f*x_freq[i2] * *n1 + *n2;
+ xout2 = xin2 - 2.f*x_freq[i2+1] * *n3 + *n4;
+ *n2 = *n1;
+ *n4 = *n3;
+ *n1 = xin1;
+ *n3 = xin2;
+ xin1 = xout1;
+ xin2 = xout2;
+ }
+ xout1 = xin1 + *(n4+1);
+ xout2 = xin2 - *(n4+2);
+ if (j>0)
+ ak[j-1] = (xout1 + xout2)*0.5f;
+ *(n4+1) = xin1;
+ *(n4+2) = xin2;
+
+ xin1 = 0.0;
+ xin2 = 0.0;
+ }
+
+}
+#endif
+
+
+#ifdef FIXED_POINT
+
+/*Makes sure the LSPs are stable*/
+void lsp_enforce_margin(spx_lsp_t *lsp, int len, spx_word16_t margin)
+{
+ int i;
+ spx_word16_t m = margin;
+ spx_word16_t m2 = 25736-margin;
+
+ if (lsp[0]<m)
+ lsp[0]=m;
+ if (lsp[len-1]>m2)
+ lsp[len-1]=m2;
+ for (i=1;i<len-1;i++)
+ {
+ if (lsp[i]<lsp[i-1]+m)
+ lsp[i]=lsp[i-1]+m;
+
+ if (lsp[i]>lsp[i+1]-m)
+ lsp[i]= SHR16(lsp[i],1) + SHR16(lsp[i+1]-m,1);
+ }
+}
+
+
+void lsp_interpolate(spx_lsp_t *old_lsp, spx_lsp_t *new_lsp, spx_lsp_t *interp_lsp, int len, int subframe, int nb_subframes)
+{
+ int i;
+ spx_word16_t tmp = DIV32_16(SHL32(EXTEND32(1 + subframe),14),nb_subframes);
+ spx_word16_t tmp2 = 16384-tmp;
+ for (i=0;i<len;i++)
+ {
+ interp_lsp[i] = MULT16_16_P14(tmp2,old_lsp[i]) + MULT16_16_P14(tmp,new_lsp[i]);
+ }
+}
+
+#else
+
+/*Makes sure the LSPs are stable*/
+void lsp_enforce_margin(spx_lsp_t *lsp, int len, spx_word16_t margin)
+{
+ int i;
+ if (lsp[0]<LSP_SCALING*margin)
+ lsp[0]=LSP_SCALING*margin;
+ if (lsp[len-1]>LSP_SCALING*(M_PI-margin))
+ lsp[len-1]=LSP_SCALING*(M_PI-margin);
+ for (i=1;i<len-1;i++)
+ {
+ if (lsp[i]<lsp[i-1]+LSP_SCALING*margin)
+ lsp[i]=lsp[i-1]+LSP_SCALING*margin;
+
+ if (lsp[i]>lsp[i+1]-LSP_SCALING*margin)
+ lsp[i]= .5f* (lsp[i] + lsp[i+1]-LSP_SCALING*margin);
+ }
+}
+
+
+void lsp_interpolate(spx_lsp_t *old_lsp, spx_lsp_t *new_lsp, spx_lsp_t *interp_lsp, int len, int subframe, int nb_subframes)
+{
+ int i;
+ float tmp = (1.0f + subframe)/nb_subframes;
+ for (i=0;i<len;i++)
+ {
+ interp_lsp[i] = (1-tmp)*old_lsp[i] + tmp*new_lsp[i];
+ }
+}
+
+#endif
--- /dev/null
+/*---------------------------------------------------------------------------*\\r
+Original Copyright\r
+ FILE........: AK2LSPD.H\r
+\r
+Modified by Jean-Marc Valin\r
+\r
+ This file contains functions for converting Linear Prediction\r
+ Coefficients (LPC) to Line Spectral Pair (LSP) and back. Note that the\r
+ LSP coefficients are not in radians format but in the x domain of the\r
+ unit circle.\r
+\r
+\*---------------------------------------------------------------------------*/\r
+/**\r
+ @file lsp.h\r
+ @brief Line Spectral Pair (LSP) functions.\r
+*/\r
+/* Speex License:\r
+\r
+ Redistribution and use in source and binary forms, with or without\r
+ modification, are permitted provided that the following conditions\r
+ are met:\r
+ \r
+ - Redistributions of source code must retain the above copyright\r
+ notice, this list of conditions and the following disclaimer.\r
+ \r
+ - Redistributions in binary form must reproduce the above copyright\r
+ notice, this list of conditions and the following disclaimer in the\r
+ documentation and/or other materials provided with the distribution.\r
+ \r
+ - Neither the name of the Xiph.org Foundation nor the names of its\r
+ contributors may be used to endorse or promote products derived from\r
+ this software without specific prior written permission.\r
+ \r
+ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS\r
+ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT\r
+ LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR\r
+ A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR\r
+ CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,\r
+ EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,\r
+ PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR\r
+ PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF\r
+ LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING\r
+ NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS\r
+ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.\r
+*/\r
+\r
+#ifndef __AK2LSPD__\r
+#define __AK2LSPD__\r
+\r
+#include "arch.h"\r
+\r
+int lpc_to_lsp (spx_coef_t *a, int lpcrdr, spx_lsp_t *freq, int nb, spx_word16_t delta, char *stack);\r
+void lsp_to_lpc(spx_lsp_t *freq, spx_coef_t *ak, int lpcrdr, char *stack);\r
+\r
+/*Added by JMV*/\r
+void lsp_enforce_margin(spx_lsp_t *lsp, int len, spx_word16_t margin);\r
+\r
+void lsp_interpolate(spx_lsp_t *old_lsp, spx_lsp_t *new_lsp, spx_lsp_t *interp_lsp, int len, int subframe, int nb_subframes);\r
+\r
+#endif /* __AK2LSPD__ */\r
--- /dev/null
+/* Copyright (C) 2002 Jean-Marc Valin
+ File: lsp_tables_nb.c
+ Codebooks for LSPs in narrowband CELP mode
+
+ Redistribution and use in source and binary forms, with or without
+ modification, are permitted provided that the following conditions are
+ met:
+
+ 1. Redistributions of source code must retain the above copyright notice,
+ this list of conditions and the following disclaimer.
+
+ 2. 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.
+
+ 3. The name of the author may not be used to endorse or promote products
+ derived from this software without specific prior written permission.
+
+ THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``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 AUTHOR 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.
+*/
+
+const signed char cdbk_nb[640]={
+30,19,38,34,40,32,46,43,58,43,
+5,-18,-25,-40,-33,-55,-52,20,34,28,
+-20,-63,-97,-92,61,53,47,49,53,75,
+-14,-53,-77,-79,0,-3,-5,19,22,26,
+-9,-53,-55,66,90,72,85,68,74,52,
+-4,-41,-58,-31,-18,-31,27,32,30,18,
+24,3,8,5,-12,-3,26,28,74,63,
+-2,-39,-67,-77,-106,-74,59,59,73,65,
+44,40,71,72,82,83,98,88,89,60,
+-6,-31,-47,-48,-13,-39,-9,7,2,79,
+-1,-39,-60,-17,87,81,65,50,45,19,
+-21,-67,-91,-87,-41,-50,7,18,39,74,
+10,-31,-28,39,24,13,23,5,56,45,
+29,10,-5,-13,-11,-35,-18,-8,-10,-8,
+-25,-71,-77,-21,2,16,50,63,87,87,
+5,-32,-40,-51,-68,0,12,6,54,34,
+5,-12,32,52,68,64,69,59,65,45,
+14,-16,-31,-40,-65,-67,41,49,47,37,
+-11,-52,-75,-84,-4,57,48,42,42,33,
+-11,-51,-68,-6,13,0,8,-8,26,32,
+-23,-53,0,36,56,76,97,105,111,97,
+-1,-28,-39,-40,-43,-54,-44,-40,-18,35,
+16,-20,-19,-28,-42,29,47,38,74,45,
+3,-29,-48,-62,-80,-104,-33,56,59,59,
+10,17,46,72,84,101,117,123,123,106,
+-7,-33,-49,-51,-70,-67,-27,-31,70,67,
+-16,-62,-85,-20,82,71,86,80,85,74,
+-19,-58,-75,-45,-29,-33,-18,-25,45,57,
+-12,-42,-5,12,28,36,52,64,81,82,
+13,-9,-27,-28,22,3,2,22,26,6,
+-6,-44,-51,2,15,10,48,43,49,34,
+-19,-62,-84,-89,-102,-24,8,17,61,68,
+39,24,23,19,16,-5,12,15,27,15,
+-8,-44,-49,-60,-18,-32,-28,52,54,62,
+-8,-48,-77,-70,66,101,83,63,61,37,
+-12,-50,-75,-64,33,17,13,25,15,77,
+1,-42,-29,72,64,46,49,31,61,44,
+-8,-47,-54,-46,-30,19,20,-1,-16,0,
+16,-12,-18,-9,-26,-27,-10,-22,53,45,
+-10,-47,-75,-82,-105,-109,8,25,49,77,
+50,65,114,117,124,118,115,96,90,61,
+-9,-45,-63,-60,-75,-57,8,11,20,29,
+0,-35,-49,-43,40,47,35,40,55,38,
+-24,-76,-103,-112,-27,3,23,34,52,75,
+8,-29,-43,12,63,38,35,29,24,8,
+25,11,1,-15,-18,-43,-7,37,40,21,
+-20,-56,-19,-19,-4,-2,11,29,51,63,
+-2,-44,-62,-75,-89,30,57,51,74,51,
+50,46,68,64,65,52,63,55,65,43,
+18,-9,-26,-35,-55,-69,3,6,8,17,
+-15,-61,-86,-97,1,86,93,74,78,67,
+-1,-38,-66,-48,48,39,29,25,17,-1,
+13,13,29,39,50,51,69,82,97,98,
+-2,-36,-46,-27,-16,-30,-13,-4,-7,-4,
+25,-5,-11,-6,-25,-21,33,12,31,29,
+-8,-38,-52,-63,-68,-89,-33,-1,10,74,
+-2,-15,59,91,105,105,101,87,84,62,
+-7,-33,-50,-35,-54,-47,25,17,82,81,
+-13,-56,-83,21,58,31,42,25,72,65,
+-24,-66,-91,-56,9,-2,21,10,69,75,
+2,-24,11,22,25,28,38,34,48,33,
+7,-29,-26,17,15,-1,14,0,-2,0,
+-6,-41,-67,6,-2,-9,19,2,85,74,
+-22,-67,-84,-71,-50,3,11,-9,2,62};
+
+const signed char cdbk_nb_low1[320]={
+-34,-52,-15,45,2,
+23,21,52,24,-33,
+-9,-1,9,-44,-41,
+-13,-17,44,22,-17,
+-6,-4,-1,22,38,
+26,16,2,50,27,
+-35,-34,-9,-41,6,
+0,-16,-34,51,8,
+-14,-31,-49,15,-33,
+45,49,33,-11,-37,
+-62,-54,45,11,-5,
+-72,11,-1,-12,-11,
+24,27,-11,-43,46,
+43,33,-12,-9,-1,
+1,-4,-23,-57,-71,
+11,8,16,17,-8,
+-20,-31,-41,53,48,
+-16,3,65,-24,-8,
+-23,-32,-37,-32,-49,
+-10,-17,6,38,5,
+-9,-17,-46,8,52,
+3,6,45,40,39,
+-7,-6,-34,-74,31,
+8,1,-16,43,68,
+-11,-19,-31,4,6,
+0,-6,-17,-16,-38,
+-16,-30,2,9,-39,
+-16,-1,43,-10,48,
+3,3,-16,-31,-3,
+62,68,43,13,3,
+-10,8,20,-56,12,
+12,-2,-18,22,-15,
+-40,-36,1,7,41,
+0,1,46,-6,-62,
+-4,-12,-2,-11,-83,
+-13,-2,91,33,-10,
+0,4,-11,-16,79,
+32,37,14,9,51,
+-21,-28,-56,-34,0,
+21,9,-26,11,28,
+-42,-54,-23,-2,-15,
+31,30,8,-39,-66,
+-39,-36,31,-28,-40,
+-46,35,40,22,24,
+33,48,23,-34,14,
+40,32,17,27,-3,
+25,26,-13,-61,-17,
+11,4,31,60,-6,
+-26,-41,-64,13,16,
+-26,54,31,-11,-23,
+-9,-11,-34,-71,-21,
+-34,-35,55,50,29,
+-22,-27,-50,-38,57,
+33,42,57,48,26,
+11,0,-49,-31,26,
+-4,-14,5,78,37,
+17,0,-49,-12,-23,
+26,14,2,2,-43,
+-17,-12,10,-8,-4,
+8,18,12,-6,20,
+-12,-6,-13,-25,34,
+15,40,49,7,8,
+13,20,20,-19,-22,
+-2,-8,2,51,-51};
+
+const signed char cdbk_nb_low2[320]={
+-6,53,-21,-24,4,
+26,17,-4,-37,25,
+17,-36,-13,31,3,
+-6,27,15,-10,31,
+28,26,-10,-10,-40,
+16,-7,15,13,41,
+-9,0,-4,50,-6,
+-7,14,38,22,0,
+-48,2,1,-13,-19,
+32,-3,-60,11,-17,
+-1,-24,-34,-1,35,
+-5,-27,28,44,13,
+25,15,42,-11,15,
+51,35,-36,20,8,
+-4,-12,-29,19,-47,
+49,-15,-4,16,-29,
+-39,14,-30,4,25,
+-9,-5,-51,-14,-3,
+-40,-32,38,5,-9,
+-8,-4,-1,-22,71,
+-3,14,26,-18,-22,
+24,-41,-25,-24,6,
+23,19,-10,39,-26,
+-27,65,45,2,-7,
+-26,-8,22,-12,16,
+15,16,-35,-5,33,
+-21,-8,0,23,33,
+34,6,21,36,6,
+-7,-22,8,-37,-14,
+31,38,11,-4,-3,
+-39,-32,-8,32,-23,
+-6,-12,16,20,-28,
+-4,23,13,-52,-1,
+22,6,-33,-40,-6,
+4,-62,13,5,-26,
+35,39,11,2,57,
+-11,9,-20,-28,-33,
+52,-5,-6,-2,22,
+-14,-16,-48,35,1,
+-58,20,13,33,-1,
+-74,56,-18,-22,-31,
+12,6,-14,4,-2,
+-9,-47,10,-3,29,
+-17,-5,61,14,47,
+-12,2,72,-39,-17,
+92,64,-53,-51,-15,
+-30,-38,-41,-29,-28,
+27,9,36,9,-35,
+-42,81,-21,20,25,
+-16,-5,-17,-35,21,
+15,-28,48,2,-2,
+9,-19,29,-40,30,
+-18,-18,18,-16,-57,
+15,-20,-12,-15,-37,
+-15,33,-39,21,-22,
+-13,35,11,13,-38,
+-63,29,23,-27,32,
+18,3,-26,42,33,
+-64,-66,-17,16,56,
+2,36,3,31,21,
+-41,-39,8,-57,14,
+37,-2,19,-36,-19,
+-23,-29,-16,1,-3,
+-8,-10,31,64,-65};
+
+const signed char cdbk_nb_high1[320]={
+-26,-8,29,21,4,
+19,-39,33,-7,-36,
+56,54,48,40,29,
+-4,-24,-42,-66,-43,
+-60,19,-2,37,41,
+-10,-37,-60,-64,18,
+-22,77,73,40,25,
+4,19,-19,-66,-2,
+11,5,21,14,26,
+-25,-86,-4,18,1,
+26,-37,10,37,-1,
+24,-12,-59,-11,20,
+-6,34,-16,-16,42,
+19,-28,-51,53,32,
+4,10,62,21,-12,
+-34,27,4,-48,-48,
+-50,-49,31,-7,-21,
+-42,-25,-4,-43,-22,
+59,2,27,12,-9,
+-6,-16,-8,-32,-58,
+-16,-29,-5,41,23,
+-30,-33,-46,-13,-10,
+-38,52,52,1,-17,
+-9,10,26,-25,-6,
+33,-20,53,55,25,
+-32,-5,-42,23,21,
+66,5,-28,20,9,
+75,29,-7,-42,-39,
+15,3,-23,21,6,
+11,1,-29,14,63,
+10,54,26,-24,-51,
+-49,7,-23,-51,15,
+-66,1,60,25,10,
+0,-30,-4,-15,17,
+19,59,40,4,-5,
+33,6,-22,-58,-70,
+-5,23,-6,60,44,
+-29,-16,-47,-29,52,
+-19,50,28,16,35,
+31,36,0,-21,6,
+21,27,22,42,7,
+-66,-40,-8,7,19,
+46,0,-4,60,36,
+45,-7,-29,-6,-32,
+-39,2,6,-9,33,
+20,-51,-34,18,-6,
+19,6,11,5,-19,
+-29,-2,42,-11,-45,
+-21,-55,57,37,2,
+-14,-67,-16,-27,-38,
+69,48,19,2,-17,
+20,-20,-16,-34,-17,
+-25,-61,10,73,45,
+16,-40,-64,-17,-29,
+-22,56,17,-39,8,
+-11,8,-25,-18,-13,
+-19,8,54,57,36,
+-17,-26,-4,6,-21,
+40,42,-4,20,31,
+53,10,-34,-53,31,
+-17,35,0,15,-6,
+-20,-63,-73,22,25,
+29,17,8,-29,-39,
+-69,18,15,-15,-5};
+
+const signed char cdbk_nb_high2[320]={
+11,47,16,-9,-46,
+-32,26,-64,34,-5,
+38,-7,47,20,2,
+-73,-99,-3,-45,20,
+70,-52,15,-6,-7,
+-82,31,21,47,51,
+39,-3,9,0,-41,
+-7,-15,-54,2,0,
+27,-31,9,-45,-22,
+-38,-24,-24,8,-33,
+23,5,50,-36,-17,
+-18,-51,-2,13,19,
+43,12,-15,-12,61,
+38,38,7,13,0,
+6,-1,3,62,9,
+27,22,-33,38,-35,
+-9,30,-43,-9,-32,
+-1,4,-4,1,-5,
+-11,-8,38,31,11,
+-10,-42,-21,-37,1,
+43,15,-13,-35,-19,
+-18,15,23,-26,59,
+1,-21,53,8,-41,
+-50,-14,-28,4,21,
+25,-28,-40,5,-40,
+-41,4,51,-33,-8,
+-8,1,17,-60,12,
+25,-41,17,34,43,
+19,45,7,-37,24,
+-15,56,-2,35,-10,
+48,4,-47,-2,5,
+-5,-54,5,-3,-33,
+-10,30,-2,-44,-24,
+-38,9,-9,42,4,
+6,-56,44,-16,9,
+-40,-26,18,-20,10,
+28,-41,-21,-4,13,
+-18,32,-30,-3,37,
+15,22,28,50,-40,
+3,-29,-64,7,51,
+-19,-11,17,-27,-40,
+-64,24,-12,-7,-27,
+3,37,48,-1,2,
+-9,-38,-34,46,1,
+27,-6,19,-13,26,
+10,34,20,25,40,
+50,-6,-7,30,9,
+-24,0,-23,71,-61,
+22,58,-34,-4,2,
+-49,-33,25,30,-8,
+-6,-16,77,2,38,
+-8,-35,-6,-30,56,
+78,31,33,-20,13,
+-39,20,22,4,21,
+-8,4,-6,10,-83,
+-41,9,-25,-43,15,
+-7,-12,-34,-39,-37,
+-33,19,30,16,-33,
+42,-25,25,-68,44,
+-15,-11,-4,23,50,
+14,4,-39,-43,20,
+-30,60,9,-20,7,
+16,19,-33,37,29,
+16,-35,7,38,-27};
--- /dev/null
+/* Copyright (C) 2002 Jean-Marc Valin */
+/**
+ @file math_approx.h
+ @brief Various math approximation functions for Speex
+*/
+/*
+ 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 MATH_APPROX_H
+#define MATH_APPROX_H
+
+#include "arch.h"
+
+#ifndef FIXED_POINT
+
+#define spx_sqrt sqrt
+#define spx_acos acos
+#define spx_exp exp
+#define spx_cos_norm(x) (cos((.5f*M_PI)*(x)))
+#define spx_atan atan
+
+/** Generate a pseudo-random number */
+static inline spx_word16_t speex_rand(spx_word16_t std, spx_int32_t *seed)
+{
+ const unsigned int jflone = 0x3f800000;
+ const unsigned int jflmsk = 0x007fffff;
+ union {int i; float f;} ran;
+ *seed = 1664525 * *seed + 1013904223;
+ ran.i = jflone | (jflmsk & *seed);
+ ran.f -= 1.5;
+ return 3.4642*std*ran.f;
+}
+
+
+#endif
+
+
+static inline spx_int16_t spx_ilog2(spx_uint32_t x)
+{
+ int r=0;
+ if (x>=(spx_int32_t)65536)
+ {
+ x >>= 16;
+ r += 16;
+ }
+ if (x>=256)
+ {
+ x >>= 8;
+ r += 8;
+ }
+ if (x>=16)
+ {
+ x >>= 4;
+ r += 4;
+ }
+ if (x>=4)
+ {
+ x >>= 2;
+ r += 2;
+ }
+ if (x>=2)
+ {
+ r += 1;
+ }
+ return r;
+}
+
+static inline spx_int16_t spx_ilog4(spx_uint32_t x)
+{
+ int r=0;
+ if (x>=(spx_int32_t)65536)
+ {
+ x >>= 16;
+ r += 8;
+ }
+ if (x>=256)
+ {
+ x >>= 8;
+ r += 4;
+ }
+ if (x>=16)
+ {
+ x >>= 4;
+ r += 2;
+ }
+ if (x>=4)
+ {
+ r += 1;
+ }
+ return r;
+}
+
+#ifdef FIXED_POINT
+
+/** Generate a pseudo-random number */
+static inline spx_word16_t speex_rand(spx_word16_t std, spx_int32_t *seed)
+{
+ spx_word32_t res;
+ *seed = 1664525 * *seed + 1013904223;
+ res = MULT16_16(EXTRACT16(SHR32(*seed,16)),std);
+ return EXTRACT16(PSHR32(SUB32(res, SHR32(res, 3)),14));
+}
+
+/* sqrt(x) ~= 0.22178 + 1.29227*x - 0.77070*x^2 + 0.25723*x^3 (for .25 < x < 1) */
+/*#define C0 3634
+#define C1 21173
+#define C2 -12627
+#define C3 4215*/
+
+/* sqrt(x) ~= 0.22178 + 1.29227*x - 0.77070*x^2 + 0.25659*x^3 (for .25 < x < 1) */
+#define C0 3634
+#define C1 21173
+#define C2 -12627
+#define C3 4204
+
+static inline spx_word16_t spx_sqrt(spx_word32_t x)
+{
+ int k;
+ spx_word32_t rt;
+ k = spx_ilog4(x)-6;
+ x = VSHR32(x, (k<<1));
+ rt = ADD16(C0, MULT16_16_Q14(x, ADD16(C1, MULT16_16_Q14(x, ADD16(C2, MULT16_16_Q14(x, (C3)))))));
+ rt = VSHR32(rt,7-k);
+ return rt;
+}
+
+/* log(x) ~= -2.18151 + 4.20592*x - 2.88938*x^2 + 0.86535*x^3 (for .5 < x < 1) */
+
+
+#define A1 16469
+#define A2 2242
+#define A3 1486
+
+static inline spx_word16_t spx_acos(spx_word16_t x)
+{
+ int s=0;
+ spx_word16_t ret;
+ spx_word16_t sq;
+ if (x<0)
+ {
+ s=1;
+ x = NEG16(x);
+ }
+ x = SUB16(16384,x);
+
+ x = x >> 1;
+ sq = MULT16_16_Q13(x, ADD16(A1, MULT16_16_Q13(x, ADD16(A2, MULT16_16_Q13(x, (A3))))));
+ ret = spx_sqrt(SHL32(EXTEND32(sq),13));
+
+ /*ret = spx_sqrt(67108864*(-1.6129e-04 + 2.0104e+00*f + 2.7373e-01*f*f + 1.8136e-01*f*f*f));*/
+ if (s)
+ ret = SUB16(25736,ret);
+ return ret;
+}
+
+
+#define K1 8192
+#define K2 -4096
+#define K3 340
+#define K4 -10
+
+static inline spx_word16_t spx_cos(spx_word16_t x)
+{
+ spx_word16_t x2;
+
+ if (x<12868)
+ {
+ x2 = MULT16_16_P13(x,x);
+ return ADD32(K1, MULT16_16_P13(x2, ADD32(K2, MULT16_16_P13(x2, ADD32(K3, MULT16_16_P13(K4, x2))))));
+ } else {
+ x = SUB16(25736,x);
+ x2 = MULT16_16_P13(x,x);
+ return SUB32(-K1, MULT16_16_P13(x2, ADD32(K2, MULT16_16_P13(x2, ADD32(K3, MULT16_16_P13(K4, x2))))));
+ }
+}
+
+#define L1 32767
+#define L2 -7651
+#define L3 8277
+#define L4 -626
+
+static inline spx_word16_t _spx_cos_pi_2(spx_word16_t x)
+{
+ spx_word16_t x2;
+
+ x2 = MULT16_16_P15(x,x);
+ return ADD16(1,MIN16(32766,ADD32(SUB16(L1,x2), MULT16_16_P15(x2, ADD32(L2, MULT16_16_P15(x2, ADD32(L3, MULT16_16_P15(L4, x2))))))));
+}
+
+static inline spx_word16_t spx_cos_norm(spx_word32_t x)
+{
+ x = x&0x0001ffff;
+ if (x>SHL32(EXTEND32(1), 16))
+ x = SUB32(SHL32(EXTEND32(1), 17),x);
+ if (x&0x00007fff)
+ {
+ if (x<SHL32(EXTEND32(1), 15))
+ {
+ return _spx_cos_pi_2(EXTRACT16(x));
+ } else {
+ return NEG32(_spx_cos_pi_2(EXTRACT16(65536-x)));
+ }
+ } else {
+ if (x&0x0000ffff)
+ return 0;
+ else if (x&0x0001ffff)
+ return -32767;
+ else
+ return 32767;
+ }
+}
+
+/*
+ K0 = 1
+ K1 = log(2)
+ K2 = 3-4*log(2)
+ K3 = 3*log(2) - 2
+*/
+#define D0 16384
+#define D1 11356
+#define D2 3726
+#define D3 1301
+/* Input in Q11 format, output in Q16 */
+static inline spx_word32_t spx_exp2(spx_word16_t x)
+{
+ int integer;
+ spx_word16_t frac;
+ integer = SHR16(x,11);
+ if (integer>14)
+ return 0x7fffffff;
+ else if (integer < -15)
+ return 0;
+ frac = SHL16(x-SHL16(integer,11),3);
+ frac = ADD16(D0, MULT16_16_Q14(frac, ADD16(D1, MULT16_16_Q14(frac, ADD16(D2 , MULT16_16_Q14(D3,frac))))));
+ return VSHR32(EXTEND32(frac), -integer-2);
+}
+
+/* Input in Q11 format, output in Q16 */
+static inline spx_word32_t spx_exp(spx_word16_t x)
+{
+ if (x>21290)
+ return 0x7fffffff;
+ else if (x<-21290)
+ return 0;
+ else
+ return spx_exp2(MULT16_16_P14(23637,x));
+}
+#define M1 32767
+#define M2 -21
+#define M3 -11943
+#define M4 4936
+
+static inline spx_word16_t spx_atan01(spx_word16_t x)
+{
+ return MULT16_16_P15(x, ADD32(M1, MULT16_16_P15(x, ADD32(M2, MULT16_16_P15(x, ADD32(M3, MULT16_16_P15(M4, x)))))));
+}
+
+#undef M1
+#undef M2
+#undef M3
+#undef M4
+
+/* Input in Q15, output in Q14 */
+static inline spx_word16_t spx_atan(spx_word32_t x)
+{
+ if (x <= 32767)
+ {
+ return SHR16(spx_atan01(x),1);
+ } else {
+ int e = spx_ilog2(x);
+ if (e>=29)
+ return 25736;
+ x = DIV32_16(SHL32(EXTEND32(32767),29-e), EXTRACT16(SHR32(x, e-14)));
+ return SUB16(25736, SHR16(spx_atan01(x),1));
+ }
+}
+#else
+
+#ifndef M_PI
+#define M_PI 3.14159265358979323846 /* pi */
+#endif
+
+#define C1 0.9999932946f
+#define C2 -0.4999124376f
+#define C3 0.0414877472f
+#define C4 -0.0012712095f
+
+
+#define SPX_PI_2 1.5707963268
+static inline spx_word16_t spx_cos(spx_word16_t x)
+{
+ if (x<SPX_PI_2)
+ {
+ x *= x;
+ return C1 + x*(C2+x*(C3+C4*x));
+ } else {
+ x = M_PI-x;
+ x *= x;
+ return NEG16(C1 + x*(C2+x*(C3+C4*x)));
+ }
+}
+
+#endif
+
+
+#endif
--- /dev/null
+/* Copyright (C) 2007 Jean-Marc Valin
+
+ File: os_support.h
+ This is the (tiny) OS abstraction layer. Aside from math.h, this is the
+ only place where system headers are allowed.
+
+ Redistribution and use in source and binary forms, with or without
+ modification, are permitted provided that the following conditions are
+ met:
+
+ 1. Redistributions of source code must retain the above copyright notice,
+ this list of conditions and the following disclaimer.
+
+ 2. 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.
+
+ 3. The name of the author may not be used to endorse or promote products
+ derived from this software without specific prior written permission.
+
+ THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``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 AUTHOR 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 OS_SUPPORT_H
+#define OS_SUPPORT_H
+
+#include <string.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+#ifdef OS_SUPPORT_CUSTOM
+#include "os_support_custom.h"
+#endif
+
+/** Speex wrapper for calloc. To do your own dynamic allocation, all you need to do is replace this function, speex_realloc and speex_free
+ NOTE: speex_alloc needs to CLEAR THE MEMORY */
+#ifndef OVERRIDE_SPEEX_ALLOC
+static inline void *speex_alloc (int size)
+{
+ /* WARNING: this is not equivalent to malloc(). If you want to use malloc()
+ or your own allocator, YOU NEED TO CLEAR THE MEMORY ALLOCATED. Otherwise
+ you will experience strange bugs */
+ return calloc(size,1);
+}
+#endif
+
+/** Same as speex_alloc, except that the area is only needed inside a Speex call (might cause problem with wideband though) */
+#ifndef OVERRIDE_SPEEX_ALLOC_SCRATCH
+static inline void *speex_alloc_scratch (int size)
+{
+ /* Scratch space doesn't need to be cleared */
+ return calloc(size,1);
+}
+#endif
+
+/** Speex wrapper for realloc. To do your own dynamic allocation, all you need to do is replace this function, speex_alloc and speex_free */
+#ifndef OVERRIDE_SPEEX_REALLOC
+static inline void *speex_realloc (void *ptr, int size)
+{
+ return realloc(ptr, size);
+}
+#endif
+
+/** Speex wrapper for calloc. To do your own dynamic allocation, all you need to do is replace this function, speex_realloc and speex_alloc */
+#ifndef OVERRIDE_SPEEX_FREE
+static inline void speex_free (void *ptr)
+{
+ free(ptr);
+}
+#endif
+
+/** Same as speex_free, except that the area is only needed inside a Speex call (might cause problem with wideband though) */
+#ifndef OVERRIDE_SPEEX_FREE_SCRATCH
+static inline void speex_free_scratch (void *ptr)
+{
+ free(ptr);
+}
+#endif
+
+/** Copy n bytes of memory from src to dst. The 0* term provides compile-time type checking */
+#ifndef OVERRIDE_SPEEX_COPY
+#define SPEEX_COPY(dst, src, n) (memcpy((dst), (src), (n)*sizeof(*(dst)) + 0*((dst)-(src)) ))
+#endif
+
+/** Copy n bytes of memory from src to dst, allowing overlapping regions. The 0* term
+ provides compile-time type checking */
+#ifndef OVERRIDE_SPEEX_MOVE
+#define SPEEX_MOVE(dst, src, n) (memmove((dst), (src), (n)*sizeof(*(dst)) + 0*((dst)-(src)) ))
+#endif
+
+/** Set n bytes of memory to value of c, starting at address s */
+#ifndef OVERRIDE_SPEEX_MEMSET
+#define SPEEX_MEMSET(dst, c, n) (memset((dst), (c), (n)*sizeof(*(dst))))
+#endif
+
+
+#ifndef OVERRIDE_SPEEX_FATAL
+static inline void _speex_fatal(const char *str, const char *file, int line)
+{
+ fprintf (stderr, "Fatal (internal) error in %s, line %d: %s\n", file, line, str);
+ exit(1);
+}
+#endif
+
+#ifndef OVERRIDE_SPEEX_WARNING
+static inline void speex_warning(const char *str)
+{
+#ifndef DISABLE_WARNINGS
+ fprintf (stderr, "warning: %s\n", str);
+#endif
+}
+#endif
+
+#ifndef OVERRIDE_SPEEX_WARNING_INT
+static inline void speex_warning_int(const char *str, int val)
+{
+#ifndef DISABLE_WARNINGS
+ fprintf (stderr, "warning: %s %d\n", str, val);
+#endif
+}
+#endif
+
+#ifndef OVERRIDE_SPEEX_NOTIFY
+static inline void speex_notify(const char *str)
+{
+#ifndef DISABLE_NOTIFICATIONS
+ fprintf (stderr, "notification: %s\n", str);
+#endif
+}
+#endif
+
+#ifndef OVERRIDE_SPEEX_PUTC
+/** Speex wrapper for putc */
+static inline void _speex_putc(int ch, void *file)
+{
+ FILE *f = (FILE *)file;
+ fprintf(f, "%c", ch);
+}
+#endif
+
+#define speex_fatal(str) _speex_fatal(str, __FILE__, __LINE__);
+#define speex_assert(cond) {if (!(cond)) {speex_fatal("assertion failed: " #cond);}}
+
+#ifndef RELEASE
+static inline void print_vec(float *vec, int len, char *name)
+{
+ int i;
+ printf ("%s ", name);
+ for (i=0;i<len;i++)
+ printf (" %f", vec[i]);
+ printf ("\n");
+}
+#endif
+
+#endif
+
--- /dev/null
+/* Copyright (C) 2002 Jean-Marc Valin
+ File: quant_lsp.c
+ LSP vector quantization
+
+ 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.
+*/
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "quant_lsp.h"
+#include "os_support.h"
+#include <math.h>
+#ifndef M_PI
+#define M_PI 3.14159265358979323846
+#endif
+
+#include "arch.h"
+
+#ifdef BFIN_ASM
+#include "quant_lsp_bfin.h"
+#endif
+
+#ifdef FIXED_POINT
+
+#define LSP_LINEAR(i) (SHL16(i+1,11))
+#define LSP_LINEAR_HIGH(i) (ADD16(MULT16_16_16(i,2560),6144))
+#define LSP_DIV_256(x) (SHL16((spx_word16_t)x, 5))
+#define LSP_DIV_512(x) (SHL16((spx_word16_t)x, 4))
+#define LSP_DIV_1024(x) (SHL16((spx_word16_t)x, 3))
+#define LSP_PI 25736
+
+#else
+
+#define LSP_LINEAR(i) (.25*(i)+.25)
+#define LSP_LINEAR_HIGH(i) (.3125*(i)+.75)
+#define LSP_SCALE 256.
+#define LSP_DIV_256(x) (0.0039062*(x))
+#define LSP_DIV_512(x) (0.0019531*(x))
+#define LSP_DIV_1024(x) (0.00097656*(x))
+#define LSP_PI M_PI
+
+#endif
+
+static void compute_quant_weights(spx_lsp_t *qlsp, spx_word16_t *quant_weight, int order)
+{
+ int i;
+ spx_word16_t tmp1, tmp2;
+ for (i=0;i<order;i++)
+ {
+ if (i==0)
+ tmp1 = qlsp[i];
+ else
+ tmp1 = qlsp[i]-qlsp[i-1];
+ if (i==order-1)
+ tmp2 = LSP_PI-qlsp[i];
+ else
+ tmp2 = qlsp[i+1]-qlsp[i];
+ if (tmp2<tmp1)
+ tmp1 = tmp2;
+#ifdef FIXED_POINT
+ quant_weight[i] = DIV32_16(81920,ADD16(300,tmp1));
+#else
+ quant_weight[i] = 10/(.04+tmp1);
+#endif
+ }
+
+}
+
+/* Note: x is modified*/
+#ifndef OVERRIDE_LSP_QUANT
+static int lsp_quant(spx_word16_t *x, const signed char *cdbk, int nbVec, int nbDim)
+{
+ int i,j;
+ spx_word32_t dist;
+ spx_word16_t tmp;
+ spx_word32_t best_dist=VERY_LARGE32;
+ int best_id=0;
+ const signed char *ptr=cdbk;
+ for (i=0;i<nbVec;i++)
+ {
+ dist=0;
+ for (j=0;j<nbDim;j++)
+ {
+ tmp=SUB16(x[j],SHL16((spx_word16_t)*ptr++,5));
+ dist=MAC16_16(dist,tmp,tmp);
+ }
+ if (dist<best_dist)
+ {
+ best_dist=dist;
+ best_id=i;
+ }
+ }
+
+ for (j=0;j<nbDim;j++)
+ x[j] = SUB16(x[j],SHL16((spx_word16_t)cdbk[best_id*nbDim+j],5));
+
+ return best_id;
+}
+#endif
+
+/* Note: x is modified*/
+#ifndef OVERRIDE_LSP_WEIGHT_QUANT
+static int lsp_weight_quant(spx_word16_t *x, spx_word16_t *weight, const signed char *cdbk, int nbVec, int nbDim)
+{
+ int i,j;
+ spx_word32_t dist;
+ spx_word16_t tmp;
+ spx_word32_t best_dist=VERY_LARGE32;
+ int best_id=0;
+ const signed char *ptr=cdbk;
+ for (i=0;i<nbVec;i++)
+ {
+ dist=0;
+ for (j=0;j<nbDim;j++)
+ {
+ tmp=SUB16(x[j],SHL16((spx_word16_t)*ptr++,5));
+ dist=MAC16_32_Q15(dist,weight[j],MULT16_16(tmp,tmp));
+ }
+ if (dist<best_dist)
+ {
+ best_dist=dist;
+ best_id=i;
+ }
+ }
+
+ for (j=0;j<nbDim;j++)
+ x[j] = SUB16(x[j],SHL16((spx_word16_t)cdbk[best_id*nbDim+j],5));
+ return best_id;
+}
+#endif
+
+void lsp_quant_nb(spx_lsp_t *lsp, spx_lsp_t *qlsp, int order, SpeexBits *bits)
+{
+ int i;
+ int id;
+ spx_word16_t quant_weight[10];
+
+ for (i=0;i<order;i++)
+ qlsp[i]=lsp[i];
+
+ compute_quant_weights(qlsp, quant_weight, order);
+
+ for (i=0;i<order;i++)
+ qlsp[i]=SUB16(qlsp[i],LSP_LINEAR(i));
+
+#ifndef FIXED_POINT
+ for (i=0;i<order;i++)
+ qlsp[i] = LSP_SCALE*qlsp[i];
+#endif
+ id = lsp_quant(qlsp, cdbk_nb, NB_CDBK_SIZE, order);
+ speex_bits_pack(bits, id, 6);
+
+ for (i=0;i<order;i++)
+ qlsp[i]*=2;
+
+ id = lsp_weight_quant(qlsp, quant_weight, cdbk_nb_low1, NB_CDBK_SIZE_LOW1, 5);
+ speex_bits_pack(bits, id, 6);
+
+ for (i=0;i<5;i++)
+ qlsp[i]*=2;
+
+ id = lsp_weight_quant(qlsp, quant_weight, cdbk_nb_low2, NB_CDBK_SIZE_LOW2, 5);
+ speex_bits_pack(bits, id, 6);
+
+ id = lsp_weight_quant(qlsp+5, quant_weight+5, cdbk_nb_high1, NB_CDBK_SIZE_HIGH1, 5);
+ speex_bits_pack(bits, id, 6);
+
+ for (i=5;i<10;i++)
+ qlsp[i]*=2;
+
+ id = lsp_weight_quant(qlsp+5, quant_weight+5, cdbk_nb_high2, NB_CDBK_SIZE_HIGH2, 5);
+ speex_bits_pack(bits, id, 6);
+
+#ifdef FIXED_POINT
+ for (i=0;i<order;i++)
+ qlsp[i]=PSHR16(qlsp[i],2);
+#else
+ for (i=0;i<order;i++)
+ qlsp[i]=qlsp[i] * .00097656;
+#endif
+
+ for (i=0;i<order;i++)
+ qlsp[i]=lsp[i]-qlsp[i];
+}
+
+void lsp_unquant_nb(spx_lsp_t *lsp, int order, SpeexBits *bits)
+{
+ int i, id;
+ for (i=0;i<order;i++)
+ lsp[i]=LSP_LINEAR(i);
+
+
+ id=speex_bits_unpack_unsigned(bits, 6);
+ for (i=0;i<10;i++)
+ lsp[i] = ADD32(lsp[i], LSP_DIV_256(cdbk_nb[id*10+i]));
+
+ id=speex_bits_unpack_unsigned(bits, 6);
+ for (i=0;i<5;i++)
+ lsp[i] = ADD16(lsp[i], LSP_DIV_512(cdbk_nb_low1[id*5+i]));
+
+ id=speex_bits_unpack_unsigned(bits, 6);
+ for (i=0;i<5;i++)
+ lsp[i] = ADD32(lsp[i], LSP_DIV_1024(cdbk_nb_low2[id*5+i]));
+
+ id=speex_bits_unpack_unsigned(bits, 6);
+ for (i=0;i<5;i++)
+ lsp[i+5] = ADD32(lsp[i+5], LSP_DIV_512(cdbk_nb_high1[id*5+i]));
+
+ id=speex_bits_unpack_unsigned(bits, 6);
+ for (i=0;i<5;i++)
+ lsp[i+5] = ADD32(lsp[i+5], LSP_DIV_1024(cdbk_nb_high2[id*5+i]));
+}
+
+
+void lsp_quant_lbr(spx_lsp_t *lsp, spx_lsp_t *qlsp, int order, SpeexBits *bits)
+{
+ int i;
+ int id;
+ spx_word16_t quant_weight[10];
+
+ for (i=0;i<order;i++)
+ qlsp[i]=lsp[i];
+
+ compute_quant_weights(qlsp, quant_weight, order);
+
+ for (i=0;i<order;i++)
+ qlsp[i]=SUB16(qlsp[i],LSP_LINEAR(i));
+#ifndef FIXED_POINT
+ for (i=0;i<order;i++)
+ qlsp[i]=qlsp[i]*LSP_SCALE;
+#endif
+ id = lsp_quant(qlsp, cdbk_nb, NB_CDBK_SIZE, order);
+ speex_bits_pack(bits, id, 6);
+
+ for (i=0;i<order;i++)
+ qlsp[i]*=2;
+
+ id = lsp_weight_quant(qlsp, quant_weight, cdbk_nb_low1, NB_CDBK_SIZE_LOW1, 5);
+ speex_bits_pack(bits, id, 6);
+
+ id = lsp_weight_quant(qlsp+5, quant_weight+5, cdbk_nb_high1, NB_CDBK_SIZE_HIGH1, 5);
+ speex_bits_pack(bits, id, 6);
+
+#ifdef FIXED_POINT
+ for (i=0;i<order;i++)
+ qlsp[i] = PSHR16(qlsp[i],1);
+#else
+ for (i=0;i<order;i++)
+ qlsp[i] = qlsp[i]*0.0019531;
+#endif
+
+ for (i=0;i<order;i++)
+ qlsp[i]=lsp[i]-qlsp[i];
+}
+
+void lsp_unquant_lbr(spx_lsp_t *lsp, int order, SpeexBits *bits)
+{
+ int i, id;
+ for (i=0;i<order;i++)
+ lsp[i]=LSP_LINEAR(i);
+
+
+ id=speex_bits_unpack_unsigned(bits, 6);
+ for (i=0;i<10;i++)
+ lsp[i] += LSP_DIV_256(cdbk_nb[id*10+i]);
+
+ id=speex_bits_unpack_unsigned(bits, 6);
+ for (i=0;i<5;i++)
+ lsp[i] += LSP_DIV_512(cdbk_nb_low1[id*5+i]);
+
+ id=speex_bits_unpack_unsigned(bits, 6);
+ for (i=0;i<5;i++)
+ lsp[i+5] += LSP_DIV_512(cdbk_nb_high1[id*5+i]);
+
+}
+
+
+#ifdef DISABLE_WIDEBAND
+void lsp_quant_high(spx_lsp_t *lsp, spx_lsp_t *qlsp, int order, SpeexBits *bits)
+{
+ speex_fatal("Wideband and Ultra-wideband are disabled");
+}
+void lsp_unquant_high(spx_lsp_t *lsp, int order, SpeexBits *bits)
+{
+ speex_fatal("Wideband and Ultra-wideband are disabled");
+}
+#else
+extern const signed char high_lsp_cdbk[];
+extern const signed char high_lsp_cdbk2[];
+
+
+void lsp_quant_high(spx_lsp_t *lsp, spx_lsp_t *qlsp, int order, SpeexBits *bits)
+{
+ int i;
+ int id;
+ spx_word16_t quant_weight[10];
+
+ for (i=0;i<order;i++)
+ qlsp[i]=lsp[i];
+
+ compute_quant_weights(qlsp, quant_weight, order);
+
+ /* quant_weight[0] = 10/(qlsp[1]-qlsp[0]);
+ quant_weight[order-1] = 10/(qlsp[order-1]-qlsp[order-2]);
+ for (i=1;i<order-1;i++)
+ {
+ tmp1 = 10/(qlsp[i]-qlsp[i-1]);
+ tmp2 = 10/(qlsp[i+1]-qlsp[i]);
+ quant_weight[i] = tmp1 > tmp2 ? tmp1 : tmp2;
+ }*/
+
+ for (i=0;i<order;i++)
+ qlsp[i]=SUB16(qlsp[i],LSP_LINEAR_HIGH(i));
+#ifndef FIXED_POINT
+ for (i=0;i<order;i++)
+ qlsp[i] = qlsp[i]*LSP_SCALE;
+#endif
+ id = lsp_quant(qlsp, high_lsp_cdbk, 64, order);
+ speex_bits_pack(bits, id, 6);
+
+ for (i=0;i<order;i++)
+ qlsp[i]*=2;
+
+ id = lsp_weight_quant(qlsp, quant_weight, high_lsp_cdbk2, 64, order);
+ speex_bits_pack(bits, id, 6);
+
+#ifdef FIXED_POINT
+ for (i=0;i<order;i++)
+ qlsp[i] = PSHR16(qlsp[i],1);
+#else
+ for (i=0;i<order;i++)
+ qlsp[i] = qlsp[i]*0.0019531;
+#endif
+
+ for (i=0;i<order;i++)
+ qlsp[i]=lsp[i]-qlsp[i];
+}
+
+void lsp_unquant_high(spx_lsp_t *lsp, int order, SpeexBits *bits)
+{
+
+ int i, id;
+ for (i=0;i<order;i++)
+ lsp[i]=LSP_LINEAR_HIGH(i);
+
+
+ id=speex_bits_unpack_unsigned(bits, 6);
+ for (i=0;i<order;i++)
+ lsp[i] += LSP_DIV_256(high_lsp_cdbk[id*order+i]);
+
+
+ id=speex_bits_unpack_unsigned(bits, 6);
+ for (i=0;i<order;i++)
+ lsp[i] += LSP_DIV_512(high_lsp_cdbk2[id*order+i]);
+}
+
+#endif
+
--- /dev/null
+/* Copyright (C) 2002 Jean-Marc Valin */
+/**
+ @file quant_lsp.h
+ @brief LSP vector quantization
+*/
+/*
+ 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 QUANT_LSP_H
+#define QUANT_LSP_H
+
+#include <speex/speex_bits.h>
+#include "arch.h"
+
+#define MAX_LSP_SIZE 20
+
+#define NB_CDBK_SIZE 64
+#define NB_CDBK_SIZE_LOW1 64
+#define NB_CDBK_SIZE_LOW2 64
+#define NB_CDBK_SIZE_HIGH1 64
+#define NB_CDBK_SIZE_HIGH2 64
+
+/*Narrowband codebooks*/
+extern const signed char cdbk_nb[];
+extern const signed char cdbk_nb_low1[];
+extern const signed char cdbk_nb_low2[];
+extern const signed char cdbk_nb_high1[];
+extern const signed char cdbk_nb_high2[];
+
+/* Quantizes narrowband LSPs with 30 bits */
+void lsp_quant_nb(spx_lsp_t *lsp, spx_lsp_t *qlsp, int order, SpeexBits *bits);
+
+/* Decodes quantized narrowband LSPs */
+void lsp_unquant_nb(spx_lsp_t *lsp, int order, SpeexBits *bits);
+
+/* Quantizes low bit-rate narrowband LSPs with 18 bits */
+void lsp_quant_lbr(spx_lsp_t *lsp, spx_lsp_t *qlsp, int order, SpeexBits *bits);
+
+/* Decodes quantized low bit-rate narrowband LSPs */
+void lsp_unquant_lbr(spx_lsp_t *lsp, int order, SpeexBits *bits);
+
+/* Quantizes high-band LSPs with 12 bits */
+void lsp_quant_high(spx_lsp_t *lsp, spx_lsp_t *qlsp, int order, SpeexBits *bits);
+
+/* Decodes high-band LSPs */
+void lsp_unquant_high(spx_lsp_t *lsp, int order, SpeexBits *bits);
+
+#endif
--- /dev/null
+/* Copyright (C) 2002 Jean-Marc Valin */
+/**
+ @file speex_bits.h
+ @brief Handles bit packing/unpacking
+*/
+/*
+ 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 BITS_H
+#define BITS_H
+/** @defgroup SpeexBits SpeexBits: Bit-stream manipulations
+ * This is the structure that holds the bit-stream when encoding or decoding
+ * with Speex. It allows some manipulations as well.
+ * @{
+ */
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/** Bit-packing data structure representing (part of) a bit-stream. */
+typedef struct SpeexBits {
+ char *chars; /**< "raw" data */
+ int nbBits; /**< Total number of bits stored in the stream*/
+ int charPtr; /**< Position of the byte "cursor" */
+ int bitPtr; /**< Position of the bit "cursor" within the current char */
+ int owner; /**< Does the struct "own" the "raw" buffer (member "chars") */
+ int overflow;/**< Set to one if we try to read past the valid data */
+ int buf_size;/**< Allocated size for buffer */
+ int reserved1; /**< Reserved for future use */
+ void *reserved2; /**< Reserved for future use */
+} SpeexBits;
+
+/** Initializes and allocates resources for a SpeexBits struct */
+void speex_bits_init(SpeexBits *bits);
+
+/** Initializes SpeexBits struct using a pre-allocated buffer*/
+void speex_bits_init_buffer(SpeexBits *bits, void *buff, int buf_size);
+
+/** Sets the bits in a SpeexBits struct to use data from an existing buffer (for decoding without copying data) */
+void speex_bits_set_bit_buffer(SpeexBits *bits, void *buff, int buf_size);
+
+/** Frees all resources associated to a SpeexBits struct. Right now this does nothing since no resources are allocated, but this could change in the future.*/
+void speex_bits_destroy(SpeexBits *bits);
+
+/** Resets bits to initial value (just after initialization, erasing content)*/
+void speex_bits_reset(SpeexBits *bits);
+
+/** Rewind the bit-stream to the beginning (ready for read) without erasing the content */
+void speex_bits_rewind(SpeexBits *bits);
+
+/** Initializes the bit-stream from the data in an area of memory */
+void speex_bits_read_from(SpeexBits *bits, char *bytes, int len);
+
+/** Append bytes to the bit-stream
+ *
+ * @param bits Bit-stream to operate on
+ * @param bytes pointer to the bytes what will be appended
+ * @param len Number of bytes of append
+ */
+void speex_bits_read_whole_bytes(SpeexBits *bits, char *bytes, int len);
+
+/** Write the content of a bit-stream to an area of memory
+ *
+ * @param bits Bit-stream to operate on
+ * @param bytes Memory location where to write the bits
+ * @param max_len Maximum number of bytes to write (i.e. size of the "bytes" buffer)
+ * @return Number of bytes written to the "bytes" buffer
+*/
+int speex_bits_write(SpeexBits *bits, char *bytes, int max_len);
+
+/** Like speex_bits_write, but writes only the complete bytes in the stream. Also removes the written bytes from the stream */
+int speex_bits_write_whole_bytes(SpeexBits *bits, char *bytes, int max_len);
+
+/** Append bits to the bit-stream
+ * @param bits Bit-stream to operate on
+ * @param data Value to append as integer
+ * @param nbBits number of bits to consider in "data"
+ */
+void speex_bits_pack(SpeexBits *bits, int data, int nbBits);
+
+/** Interpret the next bits in the bit-stream as a signed integer
+ *
+ * @param bits Bit-stream to operate on
+ * @param nbBits Number of bits to interpret
+ * @return A signed integer represented by the bits read
+ */
+int speex_bits_unpack_signed(SpeexBits *bits, int nbBits);
+
+/** Interpret the next bits in the bit-stream as an unsigned integer
+ *
+ * @param bits Bit-stream to operate on
+ * @param nbBits Number of bits to interpret
+ * @return An unsigned integer represented by the bits read
+ */
+unsigned int speex_bits_unpack_unsigned(SpeexBits *bits, int nbBits);
+
+/** Returns the number of bytes in the bit-stream, including the last one even if it is not "full"
+ *
+ * @param bits Bit-stream to operate on
+ * @return Number of bytes in the stream
+ */
+int speex_bits_nbytes(SpeexBits *bits);
+
+/** Same as speex_bits_unpack_unsigned, but without modifying the cursor position
+ *
+ * @param bits Bit-stream to operate on
+ * @param nbBits Number of bits to look for
+ * @return Value of the bits peeked, interpreted as unsigned
+ */
+unsigned int speex_bits_peek_unsigned(SpeexBits *bits, int nbBits);
+
+/** Get the value of the next bit in the stream, without modifying the
+ * "cursor" position
+ *
+ * @param bits Bit-stream to operate on
+ * @return Value of the bit peeked (one bit only)
+ */
+int speex_bits_peek(SpeexBits *bits);
+
+/** Advances the position of the "bit cursor" in the stream
+ *
+ * @param bits Bit-stream to operate on
+ * @param n Number of bits to advance
+ */
+void speex_bits_advance(SpeexBits *bits, int n);
+
+/** Returns the number of bits remaining to be read in a stream
+ *
+ * @param bits Bit-stream to operate on
+ * @return Number of bits that can still be read from the stream
+ */
+int speex_bits_remaining(SpeexBits *bits);
+
+/** Insert a terminator so that the data can be sent as a packet while auto-detecting
+ * the number of frames in each packet
+ *
+ * @param bits Bit-stream to operate on
+ */
+void speex_bits_insert_terminator(SpeexBits *bits);
+
+#ifdef __cplusplus
+}
+#endif
+
+/* @} */
+#endif
--- /dev/null
+/* speex_types.h taken from libogg */
+/********************************************************************
+ * *
+ * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. *
+ * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS *
+ * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
+ * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. *
+ * *
+ * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2002 *
+ * by the Xiph.Org Foundation http://www.xiph.org/ *
+ * *
+ ********************************************************************
+
+ function: #ifdef jail to whip a few platforms into the UNIX ideal.
+ last mod: $Id: os_types.h 7524 2004-08-11 04:20:36Z conrad $
+
+ ********************************************************************/
+/**
+ @file speex_types.h
+ @brief Speex types
+*/
+#ifndef _SPEEX_TYPES_H
+#define _SPEEX_TYPES_H
+
+#if defined(_WIN32)
+
+# if defined(__CYGWIN__)
+# include <_G_config.h>
+ typedef _G_int32_t spx_int32_t;
+ typedef _G_uint32_t spx_uint32_t;
+ typedef _G_int16_t spx_int16_t;
+ typedef _G_uint16_t spx_uint16_t;
+# elif defined(__MINGW32__)
+ typedef short spx_int16_t;
+ typedef unsigned short spx_uint16_t;
+ typedef int spx_int32_t;
+ typedef unsigned int spx_uint32_t;
+# elif defined(__MWERKS__)
+ typedef int spx_int32_t;
+ typedef unsigned int spx_uint32_t;
+ typedef short spx_int16_t;
+ typedef unsigned short spx_uint16_t;
+# else
+ /* MSVC/Borland */
+ typedef __int32 spx_int32_t;
+ typedef unsigned __int32 spx_uint32_t;
+ typedef __int16 spx_int16_t;
+ typedef unsigned __int16 spx_uint16_t;
+# endif
+
+#elif defined(__MACOS__)
+
+# include <sys/types.h>
+ typedef SInt16 spx_int16_t;
+ typedef UInt16 spx_uint16_t;
+ typedef SInt32 spx_int32_t;
+ typedef UInt32 spx_uint32_t;
+
+#elif (defined(__APPLE__) && defined(__MACH__)) /* MacOS X Framework build */
+
+# include <sys/types.h>
+ typedef int16_t spx_int16_t;
+ typedef u_int16_t spx_uint16_t;
+ typedef int32_t spx_int32_t;
+ typedef u_int32_t spx_uint32_t;
+
+#elif defined(__BEOS__)
+
+ /* Be */
+# include <inttypes.h>
+ typedef int16_t spx_int16_t;
+ typedef u_int16_t spx_uint16_t;
+ typedef int32_t spx_int32_t;
+ typedef u_int32_t spx_uint32_t;
+
+#elif defined (__EMX__)
+
+ /* OS/2 GCC */
+ typedef short spx_int16_t;
+ typedef unsigned short spx_uint16_t;
+ typedef int spx_int32_t;
+ typedef unsigned int spx_uint32_t;
+
+#elif defined (DJGPP)
+
+ /* DJGPP */
+ typedef short spx_int16_t;
+ typedef int spx_int32_t;
+ typedef unsigned int spx_uint32_t;
+
+#elif defined(R5900)
+
+ /* PS2 EE */
+ typedef int spx_int32_t;
+ typedef unsigned spx_uint32_t;
+ typedef short spx_int16_t;
+
+#elif defined(__SYMBIAN32__)
+
+ /* Symbian GCC */
+ typedef signed short spx_int16_t;
+ typedef unsigned short spx_uint16_t;
+ typedef signed int spx_int32_t;
+ typedef unsigned int spx_uint32_t;
+
+#elif defined(CONFIG_TI_C54X) || defined (CONFIG_TI_C55X)
+
+ typedef short spx_int16_t;
+ typedef unsigned short spx_uint16_t;
+ typedef long spx_int32_t;
+ typedef unsigned long spx_uint32_t;
+
+#elif defined(CONFIG_TI_C6X)
+
+ typedef short spx_int16_t;
+ typedef unsigned short spx_uint16_t;
+ typedef int spx_int32_t;
+ typedef unsigned int spx_uint32_t;
+
+#else
+
+# include <speex/speex_config_types.h>
+
+#endif
+
+#endif /* _SPEEX_TYPES_H */
--- /dev/null
+/* Copyright (C) 2002 Jean-Marc Valin */
+/**
+ @file stack_alloc.h
+ @brief Temporary memory allocation on stack
+*/
+/*
+ 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 STACK_ALLOC_H
+#define STACK_ALLOC_H
+
+#ifdef USE_ALLOCA
+# ifdef WIN32
+# include <malloc.h>
+# else
+# ifdef HAVE_ALLOCA_H
+# include <alloca.h>
+# else
+# include <stdlib.h>
+# endif
+# endif
+#endif
+
+/**
+ * @def ALIGN(stack, size)
+ *
+ * Aligns the stack to a 'size' boundary
+ *
+ * @param stack Stack
+ * @param size New size boundary
+ */
+
+/**
+ * @def PUSH(stack, size, type)
+ *
+ * Allocates 'size' elements of type 'type' on the stack
+ *
+ * @param stack Stack
+ * @param size Number of elements
+ * @param type Type of element
+ */
+
+/**
+ * @def VARDECL(var)
+ *
+ * Declare variable on stack
+ *
+ * @param var Variable to declare
+ */
+
+/**
+ * @def ALLOC(var, size, type)
+ *
+ * Allocate 'size' elements of 'type' on stack
+ *
+ * @param var Name of variable to allocate
+ * @param size Number of elements
+ * @param type Type of element
+ */
+
+#ifdef ENABLE_VALGRIND
+
+#include <valgrind/memcheck.h>
+
+#define ALIGN(stack, size) ((stack) += ((size) - (long)(stack)) & ((size) - 1))
+
+#define PUSH(stack, size, type) (VALGRIND_MAKE_NOACCESS(stack, 1000),ALIGN((stack),sizeof(type)),VALGRIND_MAKE_WRITABLE(stack, ((size)*sizeof(type))),(stack)+=((size)*sizeof(type)),(type*)((stack)-((size)*sizeof(type))))
+
+#else
+
+#define ALIGN(stack, size) ((stack) += ((size) - (long)(stack)) & ((size) - 1))
+
+#define PUSH(stack, size, type) (ALIGN((stack),sizeof(type)),(stack)+=((size)*sizeof(type)),(type*)((stack)-((size)*sizeof(type))))
+
+#endif
+
+#if defined(VAR_ARRAYS)
+#define VARDECL(var)
+#define ALLOC(var, size, type) type var[size]
+#elif defined(USE_ALLOCA)
+#define VARDECL(var) var
+#define ALLOC(var, size, type) var = alloca(sizeof(type)*(size))
+#else
+#define VARDECL(var) var
+#define ALLOC(var, size, type) var = PUSH(stack, size, type)
+#endif
+
+
+#endif
--- /dev/null
+/*---------------------------------------------------------------------------*\
+
+ FILE........: comp.h
+ AUTHOR......: David Rowe
+ DATE CREATED: 24/08/09
+
+ Complex number definition.
+
+\*---------------------------------------------------------------------------*/
+
+/*
+ Copyright (C) 2009 David Rowe
+
+ All rights reserved.
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License version 2, as
+ published by the Free Software Foundation. This program is
+ distributed in the hope that it will be useful, but WITHOUT ANY
+ WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
+ License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program; if not, write to the Free Software
+ Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+*/
+
+#ifndef __COMP__
+#define __COMP__
+
+/* Complex number */
+
+typedef struct {
+ float real;
+ float imag;
+} COMP;
+
+#endif
-CFLAGS = -I../src -Wall -g
+CFLAGS = -I. -I../src -I../speex -Wall -g -DFLOATING_POINT -DVAR_ARRAYS
-all: genres
+all: genres lsptest
genres: genres.o ../src/lpc.o
gcc $(CFLAGS) -o genres genres.o ../src/lpc.o -lm
+LSP_TEST_OBJ = lsptest.o ../src/lpc.o ../speex/lsp.o sd.o ../src/four1.o \
+ ../speex/quant_lsp.o ../speex/bits.o ../speex/lsp_tables_nb.o \
+ ../speex/high_lsp_tables.c
+
+lsptest: $(LSP_TEST_OBJ)
+ gcc $(CFLAGS) -o lsptest $(LSP_TEST_OBJ) -lm
+
%.o : %.c
$(CC) -c $(CFLAGS) $< -o $@
--- /dev/null
+/*---------------------------------------------------------------------------*\
+
+ FILE........: lsptest.c
+ AUTHOR......: David Rowe
+ DATE CREATED: 24/8/09
+
+ Test Speech LPC to LSP conversion and quantisation.
+
+\*---------------------------------------------------------------------------*/
+
+/*
+ Copyright (C) 2009 David Rowe
+
+ All rights reserved.
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License version 2, as
+ published by the Free Software Foundation. This program is
+ distributed in the hope that it will be useful, but WITHOUT ANY
+ WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
+ License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program; if not, write to the Free Software
+ Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+*/
+
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <lpc.h>
+#include <lsp.h>
+#include <sd.h>
+#include <speex_bits.h>
+#include <quant_lsp.h>
+
+#define N 160
+#define P 10
+
+#define LPC_FLOOR 0.0002 /* autocorrelation floor */
+#define LSP_DELTA1 0.2 /* grid spacing for LSP root searches */
+#define NDFT 256 /* DFT size for SD calculation */
+
+/* Speex lag window */
+
+const float lag_window[11] = {
+ 1.00000, 0.99716, 0.98869, 0.97474, 0.95554, 0.93140, 0.90273, 0.86998,
+ 0.83367, 0.79434, 0.75258
+};
+
+/*---------------------------------------------------------------------------*\
+
+ find_aks_for_lsp()
+
+ This function takes a frame of samples, and determines the linear
+ prediction coefficients for that frame of samples. Modified version of
+ find_aks from lpc.c to include autocorrelation noise floor and lag window
+ to match Speex processing steps prior to LSP conversion.
+
+\*---------------------------------------------------------------------------*/
+
+void find_aks_for_lsp(
+ float Sn[], /* Nsam samples with order sample memory */
+ float a[], /* order+1 LPCs with first coeff 1.0 */
+ int Nsam, /* number of input speech samples */
+ int order, /* order of the LPC analysis */
+ float *E /* residual energy */
+)
+{
+ float Wn[N]; /* windowed frame of Nsam speech samples */
+ float R[P+1]; /* order+1 autocorrelation values of Sn[] */
+ int i;
+
+ hanning_window(Sn,Wn,Nsam);
+
+ autocorrelate(Wn,R,Nsam,order);
+ R[0] += LPC_FLOOR;
+ assert(order == 10); /* lag window only defined for order == 10 */
+ for(i=0; i<=order; i++)
+ R[i] *= lag_window[i];
+ levinson_durbin(R,a,order);
+
+ *E = 0.0;
+ for(i=0; i<=order; i++)
+ *E += a[i]*R[i];
+ if (*E < 0.0)
+ *E = 1E-12;
+}
+
+/*---------------------------------------------------------------------------*\
+
+ MAIN
+
+\*---------------------------------------------------------------------------*/
+
+int main(int argc, char *argv[])
+{
+ FILE *fin; /* input speech files */
+ short buf[N]; /* buffer of 16 bit speech samples */
+ float Sn[P+N]; /* input speech samples */
+ float E;
+ float ak[P+1]; /* LP coeffs */
+ float ak_[P+1]; /* quantised LP coeffs */
+ float lsp[P];
+ float lsp_[P]; /* quantised LSPs */
+ int roots; /* number of LSP roots found */
+ int frames; /* frames processed so far */
+ int i; /* loop variables */
+
+ SpeexBits bits;
+
+ float sd; /* SD for this frame */
+ float totsd; /* accumulated SD so far */
+ int gt2,gt4; /* number of frames > 2 and 4 dB SD */
+ int unstables; /* number of unstable LSP frames */
+
+ if (argc < 2) {
+ printf("usage: %s InputFile\n", argv[0]);
+ exit(0);
+ }
+
+ /* Open files */
+
+ if ((fin = fopen(argv[1],"rb")) == NULL) {
+ printf("Error opening input file: %s\n",argv[1]);
+ exit(0);
+ }
+
+ /* Initialise */
+
+ frames = 0;
+ for(i=0; i<P; i++) {
+ Sn[i] = 0.0;
+ }
+ ak_[0] = 1.0;
+
+ speex_bits_init(&bits);
+
+ totsd = 0.0;
+ unstables = 0;
+ gt2 = 0; gt4 = 0;
+
+ /* Main loop */
+
+ while( (fread(buf,sizeof(short),N,fin)) == N) {
+ frames++;
+ for(i=0; i<N; i++)
+ Sn[P+i] = (float)buf[i];
+
+ /* convert to LSP domain and back */
+
+ find_aks(&Sn[P], ak, N, P, &E);
+ roots = lpc_to_lsp(&ak[1], P , lsp, 10, LSP_DELTA1, NULL);
+ if (roots == P) {
+
+ speex_bits_reset(&bits);
+ lsp_quant_lbr(lsp, lsp_, P, &bits);
+ lsp_to_lpc(lsp_, &ak_[1], P, NULL);
+
+ /* measure spectral distortion */
+ sd = spectral_dist(ak, ak_, P, NDFT);
+ if (sd > 2.0) gt2++;
+ if (sd > 4.0) gt4++;
+ totsd += sd;
+ }
+ else
+ unstables++;
+ }
+
+ fclose(fin);
+
+ printf("frames = %d Av sd = %3.2f dB", frames, totsd/frames);
+ printf(" >2 dB %3.2f%% >4 dB %3.2f%% unstables: %d\n",gt2*100.0/frames,
+ gt4*100.0/frames, unstables);
+
+ return 0;
+}
+
--- /dev/null
+/*--------------------------------------------------------------------------*\
+
+ FILE........: sd.c
+ AUTHOR......: David Rowe
+ DATE CREATED: 20/7/93
+
+ Function to determine spectral distortion between two sets of LPCs.
+
+\*--------------------------------------------------------------------------*/
+
+/*
+ Copyright (C) 2009 David Rowe
+
+ All rights reserved.
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License version 2, as
+ published by the Free Software Foundation. This program is
+ distributed in the hope that it will be useful, but WITHOUT ANY
+ WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
+ License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program; if not, write to the Free Software
+ Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+*/
+
+#define MAX_N 2048 /* maximum DFT size */
+
+#include <math.h>
+#include "four1.h"
+#include "comp.h"
+#include "sd.h"
+
+/*---------------------------------------------------------------------------*\
+
+ FUNCTION....: spectral_dist()
+
+ AUTHOR......: David Rowe
+ DATE CREATED: 20/7/93
+
+ This function returns the soectral distoertion between two
+ sets of LPCs.
+
+\*---------------------------------------------------------------------------*/
+
+float spectral_dist(float ak1[], float ak2[], int p, int n)
+/* float ak1[]; unquantised set of p+1 LPCs */
+/* float ak2[]; quantised set of p+1 LPCs */
+/* int p; LP order */
+/* int n; DFT size to use for SD calculations (power of 2) */
+{
+ COMP A1[MAX_N]; /* DFT of ak1[] */
+ COMP A2[MAX_N]; /* DFT of ak2[] */
+ float P1,P2; /* power of current bin */
+ float sd;
+ int i;
+
+ for(i=0; i<n; i++) {
+ A1[i].real = 0.0;
+ A1[i].imag = 0.0;
+ A2[i].real = 0.0;
+ A2[i].imag = 0.0;
+ }
+
+ for(i=0; i<p+1; i++) {
+ A1[i].real = ak1[i];
+ A2[i].real = ak2[i];
+ }
+
+ four1(&A1[-1].imag,n,-1);
+ four1(&A2[-1].imag,n,-1);
+
+ sd = 0.0;
+ for(i=0; i<n; i++) {
+ P1 = A1[i].real*A1[i].real + A1[i].imag*A1[i].imag;
+ P2 = A2[i].real*A2[i].real + A2[i].imag*A2[i].imag;
+ sd += pow(log10(P2/P1),2.0);
+ }
+ sd = 10.0*sqrt(sd/n); /* sd in dB */
+
+ return(sd);
+}
--- /dev/null
+/*--------------------------------------------------------------------------*\
+
+ FILE........: sd.h
+ AUTHOR......: David Rowe
+ DATE CREATED: 22/7/93
+
+ Function to determine spectral distortion between two sets of LPCs.
+
+\*--------------------------------------------------------------------------*/
+
+/*
+ Copyright (C) 2009 David Rowe
+
+ All rights reserved.
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License version 2, as
+ published by the Free Software Foundation. This program is
+ distributed in the hope that it will be useful, but WITHOUT ANY
+ WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
+ License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program; if not, write to the Free Software
+ Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+*/
+
+#ifndef __SD__
+#define __SD__
+
+float spectral_dist(float ak1[], float ak2[], int p, int n);
+
+#endif /* __SD__ */
projects / freetel-svn-tracking.git / commitdiff