From 1f330c1e2a3703da9cfb5cd6543cd607c8e2f882 Mon Sep 17 00:00:00 2001 From: drowe67 Date: Mon, 24 Aug 2009 09:34:39 +0000 Subject: [PATCH] Speex LSP quantiser working and LSP unit test running git-svn-id: https://svn.code.sf.net/p/freetel/code@30 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2/README.txt | 2 + codec2/speex/arch.h | 239 ++++++++++++ codec2/speex/bits.c | 372 +++++++++++++++++++ codec2/speex/fixed_generic.h | 106 ++++++ codec2/speex/high_lsp_tables.c | 163 ++++++++ codec2/speex/lsp.c | 656 +++++++++++++++++++++++++++++++++ codec2/speex/lsp.h | 60 +++ codec2/speex/lsp_tables_nb.c | 360 ++++++++++++++++++ codec2/speex/math_approx.h | 332 +++++++++++++++++ codec2/speex/os_support.h | 169 +++++++++ codec2/speex/quant_lsp.c | 385 +++++++++++++++++++ codec2/speex/quant_lsp.h | 74 ++++ codec2/speex/speex_bits.h | 174 +++++++++ codec2/speex/speex_types.h | 126 +++++++ codec2/speex/stack_alloc.h | 115 ++++++ codec2/src/comp.h | 39 ++ codec2/unittest/Makefile | 11 +- codec2/unittest/lsptest.c | 179 +++++++++ codec2/unittest/sd.c | 84 +++++ codec2/unittest/sd.h | 34 ++ 20 files changed, 3678 insertions(+), 2 deletions(-) create mode 100644 codec2/speex/arch.h create mode 100644 codec2/speex/bits.c create mode 100644 codec2/speex/fixed_generic.h create mode 100644 codec2/speex/high_lsp_tables.c create mode 100644 codec2/speex/lsp.c create mode 100644 codec2/speex/lsp.h create mode 100644 codec2/speex/lsp_tables_nb.c create mode 100644 codec2/speex/math_approx.h create mode 100644 codec2/speex/os_support.h create mode 100644 codec2/speex/quant_lsp.c create mode 100644 codec2/speex/quant_lsp.h create mode 100644 codec2/speex/speex_bits.h create mode 100644 codec2/speex/speex_types.h create mode 100644 codec2/speex/stack_alloc.h create mode 100644 codec2/src/comp.h create mode 100644 codec2/unittest/lsptest.c create mode 100644 codec2/unittest/sd.c create mode 100644 codec2/unittest/sd.h diff --git a/codec2/README.txt b/codec2/README.txt index 8ace7817..27b19bec 100644 --- a/codec2/README.txt +++ b/codec2/README.txt @@ -39,6 +39,8 @@ Plan 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 diff --git a/codec2/speex/arch.h b/codec2/speex/arch.h new file mode 100644 index 00000000..d38c36ce --- /dev/null +++ b/codec2/speex/arch.h @@ -0,0 +1,239 @@ +/* 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 diff --git a/codec2/speex/bits.c b/codec2/speex/bits.c new file mode 100644 index 00000000..f61c9329 --- /dev/null +++ b/codec2/speex/bits.c @@ -0,0 +1,372 @@ +/* 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 +#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<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;ichars[i]=HTOLS(chars[i]); + + bits->nbBits=nchars<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<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;ichars[pos+i]=HTOLS(chars[i]); + bits->nbBits+=nchars<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;ichars[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;ichars[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)<charPtr<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<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<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<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<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); +} diff --git a/codec2/speex/fixed_generic.h b/codec2/speex/fixed_generic.h new file mode 100644 index 00000000..3fb096ed --- /dev/null +++ b/codec2/speex/fixed_generic.h @@ -0,0 +1,106 @@ +/* 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 diff --git a/codec2/speex/high_lsp_tables.c b/codec2/speex/high_lsp_tables.c new file mode 100644 index 00000000..e82e8755 --- /dev/null +++ b/codec2/speex/high_lsp_tables.c @@ -0,0 +1,163 @@ +/* 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}; diff --git a/codec2/speex/lsp.c b/codec2/speex/lsp.c new file mode 100644 index 00000000..a73d8835 --- /dev/null +++ b/codec2/speex/lsp.c @@ -0,0 +1,656 @@ +/*---------------------------------------------------------------------------*\ +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 +#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=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= -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 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;i0) + 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]m2) + lsp[len-1]=m2; + for (i=1;ilsp[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;iLSP_SCALING*(M_PI-margin)) + lsp[len-1]=LSP_SCALING*(M_PI-margin); + for (i=1;ilsp[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=(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 (x14) + 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 +#include +#include + +#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 +#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 tmp2 ? tmp1 : tmp2; + }*/ + + for (i=0;i +#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 diff --git a/codec2/speex/speex_bits.h b/codec2/speex/speex_bits.h new file mode 100644 index 00000000..a26fb4ce --- /dev/null +++ b/codec2/speex/speex_bits.h @@ -0,0 +1,174 @@ +/* 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 diff --git a/codec2/speex/speex_types.h b/codec2/speex/speex_types.h new file mode 100644 index 00000000..852fed80 --- /dev/null +++ b/codec2/speex/speex_types.h @@ -0,0 +1,126 @@ +/* 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 + 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 + 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 + 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 + +#endif + +#endif /* _SPEEX_TYPES_H */ diff --git a/codec2/speex/stack_alloc.h b/codec2/speex/stack_alloc.h new file mode 100644 index 00000000..5264e666 --- /dev/null +++ b/codec2/speex/stack_alloc.h @@ -0,0 +1,115 @@ +/* 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 +# else +# ifdef HAVE_ALLOCA_H +# include +# else +# include +# 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 + +#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 diff --git a/codec2/src/comp.h b/codec2/src/comp.h new file mode 100644 index 00000000..4e2eb3a2 --- /dev/null +++ b/codec2/src/comp.h @@ -0,0 +1,39 @@ +/*---------------------------------------------------------------------------*\ + + 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 diff --git a/codec2/unittest/Makefile b/codec2/unittest/Makefile index 898c13cc..eb9d67e4 100644 --- a/codec2/unittest/Makefile +++ b/codec2/unittest/Makefile @@ -1,9 +1,16 @@ -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 $@ diff --git a/codec2/unittest/lsptest.c b/codec2/unittest/lsptest.c new file mode 100644 index 00000000..81d23670 --- /dev/null +++ b/codec2/unittest/lsptest.c @@ -0,0 +1,179 @@ +/*---------------------------------------------------------------------------*\ + + 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 +#include +#include +#include +#include +#include +#include +#include + +#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 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; +} + diff --git a/codec2/unittest/sd.c b/codec2/unittest/sd.c new file mode 100644 index 00000000..2cdc6d82 --- /dev/null +++ b/codec2/unittest/sd.c @@ -0,0 +1,84 @@ +/*--------------------------------------------------------------------------*\ + + 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 +#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