diff options
Diffstat (limited to 'gr-codec2-vocoder/src/lib')
72 files changed, 8500 insertions, 0 deletions
diff --git a/gr-codec2-vocoder/src/lib/Makefile.am b/gr-codec2-vocoder/src/lib/Makefile.am new file mode 100644 index 000000000..cef4b20b0 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/Makefile.am @@ -0,0 +1,61 @@ +# +# Copyright 2004,2005,2008,2009,2010 Free Software Foundation, Inc. +# +# This file is part of GNU Radio +# +# GNU Radio is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 3, or (at your option) +# any later version. +# +# GNU Radio 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 GNU Radio; see the file COPYING. If not, write to +# the Free Software Foundation, Inc., 51 Franklin Street, +# Boston, MA 02110-1301, USA. +# + +include $(top_srcdir)/Makefile.common +include $(top_srcdir)/Makefile.swig + +SUBDIRS = codec2 . + +AM_CPPFLAGS = $(STD_DEFINES_AND_INCLUDES) $(PYTHON_CPPFLAGS) $(WITH_INCLUDES) + +# C/C++ headers get installed in ${prefix}/include/gnuradio +grinclude_HEADERS = \ + codec2_decode_ps.h \ + codec2_encode_sp.h + +lib_LTLIBRARIES = libgnuradio-codec2-vocoder.la + +libgnuradio_codec2_vocoder_la_SOURCES = \ + codec2_decode_ps.cc \ + codec2_encode_sp.cc + +libgnuradio_codec2_vocoder_la_LIBADD = \ + $(GNURADIO_CORE_LA) \ + codec2/libcodec2.la + +libgnuradio_codec2_vocoder_la_LDFLAGS = $(NO_UNDEFINED) $(LTVERSIONFLAGS) + + +# SWIG interface and library +TOP_SWIG_IFILES = \ + codec2.i + +# Install so that they end up available as: +# import gnuradio.vocoder.codec2 +# This ends up at: +# ${prefix}/lib/python${python_version}/site-packages/gnuradio/vocoder +codec2_pythondir_category = \ + gnuradio/vocoder + +# additional libraries for linking with the SWIG-generated library +codec2_la_swig_libadd = \ + libgnuradio-codec2-vocoder.la + diff --git a/gr-codec2-vocoder/src/lib/Makefile.swig.gen b/gr-codec2-vocoder/src/lib/Makefile.swig.gen new file mode 100644 index 000000000..bb81f107d --- /dev/null +++ b/gr-codec2-vocoder/src/lib/Makefile.swig.gen @@ -0,0 +1,145 @@ +# -*- Makefile -*- +# +# Copyright 2009 Free Software Foundation, Inc. +# +# This file is part of GNU Radio +# +# GNU Radio is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 3, or (at your option) +# any later version. +# +# GNU Radio 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 GNU Radio; see the file COPYING. If not, write to +# the Free Software Foundation, Inc., 51 Franklin Street, +# Boston, MA 02110-1301, USA. +# + +# Makefile.swig.gen for codec2.i + +## Default install locations for these files: +## +## Default location for the Python directory is: +## ${prefix}/lib/python${python_version}/site-packages/[category]/codec2 +## Default location for the Python exec directory is: +## ${exec_prefix}/lib/python${python_version}/site-packages/[category]/codec2 +## +## The following can be overloaded to change the install location, but +## this has to be done in the including Makefile.am -before- +## Makefile.swig is included. + +codec2_pythondir_category ?= gnuradio/codec2 +codec2_pylibdir_category ?= $(codec2_pythondir_category) +codec2_pythondir = $(pythondir)/$(codec2_pythondir_category) +codec2_pylibdir = $(pyexecdir)/$(codec2_pylibdir_category) + +# The .so libraries for the guile modules get installed whereever guile +# is installed, usually /usr/lib/guile/gnuradio/ +# FIXME: determince whether these should be installed with gnuradio. +codec2_scmlibdir = $(libdir) + +# The scm files for the guile modules get installed where ever guile +# is installed, usually /usr/share/guile/site/codec2 +# FIXME: determince whether these should be installed with gnuradio. +codec2_scmdir = $(guiledir) + +## SWIG headers are always installed into the same directory. + +codec2_swigincludedir = $(swigincludedir) + +## This is a template file for a "generated" Makefile addition (in +## this case, "Makefile.swig.gen"). By including the top-level +## Makefile.swig, this file will be used to generate the SWIG +## dependencies. Assign the variable TOP_SWIG_FILES to be the list of +## SWIG .i files to generated wrappings for; there can be more than 1 +## so long as the names are unique (no sorting is done on the +## TOP_SWIG_FILES list). This file explicitly assumes that a SWIG .i +## file will generate .cc, .py, and possibly .h files -- meaning that +## all of these files will have the same base name (that provided for +## the SWIG .i file). +## +## This code is setup to ensure parallel MAKE ("-j" or "-jN") does the +## right thing. For more info, see < +## http://sources.redhat.com/automake/automake.html#Multiple-Outputs > + +## Other cleaned files: dependency files generated by SWIG or this Makefile + +MOSTLYCLEANFILES += $(DEPDIR)/*.S* + +## Various SWIG variables. These can be overloaded in the including +## Makefile.am by setting the variable value there, then including +## Makefile.swig . + +codec2_swiginclude_HEADERS = \ + codec2.i \ + $(codec2_swiginclude_headers) + +if PYTHON +codec2_pylib_LTLIBRARIES = \ + _codec2.la + +_codec2_la_SOURCES = \ + python/codec2.cc \ + $(codec2_la_swig_sources) + +codec2_python_PYTHON = \ + codec2.py \ + $(codec2_python) + +_codec2_la_LIBADD = \ + $(STD_SWIG_LA_LIB_ADD) \ + $(codec2_la_swig_libadd) + +_codec2_la_LDFLAGS = \ + $(STD_SWIG_LA_LD_FLAGS) \ + $(codec2_la_swig_ldflags) + +_codec2_la_CXXFLAGS = \ + $(STD_SWIG_CXX_FLAGS) \ + -I$(top_builddir) \ + $(codec2_la_swig_cxxflags) + +python/codec2.cc: codec2.py +codec2.py: codec2.i + +# Include the python dependencies for this file +-include python/codec2.d + +endif # end of if python + +if GUILE + +codec2_scmlib_LTLIBRARIES = \ + libguile-gnuradio-codec2.la +libguile_gnuradio_codec2_la_SOURCES = \ + guile/codec2.cc \ + $(codec2_la_swig_sources) +nobase_codec2_scm_DATA = \ + gnuradio/codec2.scm \ + gnuradio/codec2-primitive.scm +libguile_gnuradio_codec2_la_LIBADD = \ + $(STD_SWIG_LA_LIB_ADD) \ + $(codec2_la_swig_libadd) +libguile_gnuradio_codec2_la_LDFLAGS = \ + $(STD_SWIG_LA_LD_FLAGS) \ + $(codec2_la_swig_ldflags) +libguile_gnuradio_codec2_la_CXXFLAGS = \ + $(STD_SWIG_CXX_FLAGS) \ + -I$(top_builddir) \ + $(codec2_la_swig_cxxflags) + +guile/codec2.cc: gnuradio/codec2.scm +gnuradio/codec2.scm: codec2.i +gnuradio/codec2-primitive.scm: gnuradio/codec2.scm + +# Include the guile dependencies for this file +-include guile/codec2.d + +endif # end of GUILE + + diff --git a/gr-codec2-vocoder/src/lib/codec2.i b/gr-codec2-vocoder/src/lib/codec2.i new file mode 100644 index 000000000..d29eba918 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2.i @@ -0,0 +1,58 @@ +/* -*- c++ -*- */ +/* + * Copyright 2005,2009 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * GNU Radio is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3, or (at your option) + * any later version. + * + * GNU Radio 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 GNU Radio; see the file COPYING. If not, write to + * the Free Software Foundation, Inc., 51 Franklin Street, + * Boston, MA 02110-1301, USA. + */ + +%include "gnuradio.i" + +%{ +#include "codec2_encode_sp.h" +#include "codec2_decode_ps.h" +%} + +GR_SWIG_BLOCK_MAGIC(codec2,encode_sp); + +codec2_encode_sp_sptr codec2_make_encode_sp (); + +class codec2_encode_sp : public gr_sync_decimator { +public: + ~codec2_encode_sp (); +}; + +// ---------------------------------------------------------------- + +GR_SWIG_BLOCK_MAGIC(codec2,decode_ps); + +codec2_decode_ps_sptr codec2_make_decode_ps (); + +class codec2_decode_ps : public gr_sync_interpolator { +public: + ~codec2_decode_ps (); +}; + +#if SWIGGUILE +%scheme %{ +(load-extension-global "libguile-gnuradio-codec2" "scm_init_gnuradio_codec2_module") +%} + +%goops %{ +(use-modules (gnuradio gnuradio_core_runtime)) +%} +#endif diff --git a/gr-codec2-vocoder/src/lib/codec2/Makefile.am b/gr-codec2-vocoder/src/lib/codec2/Makefile.am new file mode 100644 index 000000000..1c7c16971 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/Makefile.am @@ -0,0 +1,103 @@ +AM_CFLAGS = -I../src -fPIC -Wall -O3 +AUTOMAKE_OPTS = gnu +NAME = codec2 +AM_CPPFLAGS = $(AM_CFLAGS) + +D=codebook + +# lsp quantisers + +CODEBOOKS= \ + $D/lsp1.txt \ + $D/lsp2.txt \ + $D/lsp3.txt \ + $D/lsp4.txt \ + $D/lsp5.txt \ + $D/lsp6.txt \ + $D/lsp7.txt \ + $D/lsp8.txt \ + $D/lsp9.txt \ + $D/lsp10.txt + +# lspd quantisers + +CODEBOOKSD= \ + $D/dlsp1.txt \ + $D/dlsp2.txt \ + $D/dlsp3.txt \ + $D/dlsp4.txt \ + $D/dlsp5.txt \ + $D/dlsp6.txt \ + $D/dlsp7.txt \ + $D/dlsp8.txt \ + $D/dlsp9.txt \ + $D/dlsp10.txt + +# lspd VQ quantisers + +CODEBOOKSDVQ= \ + $D/dlsp1.txt \ + $D/dlsp2.txt \ + $D/dlsp3.txt \ + $D/dlsp4.txt \ + $D/dlsp5.txt + +noinst_PROGRAMS = generate_codebook + +codebook.$(OBJEXT): codebook.c +codebookd.$(OBJEXT): codebookd.c +codebookdvq.$(OBJEXT): codebookdvq.c + +codebook.lo: codebook.c + +codebook.c: generate_codebook $(CODEBOOKS) + ./generate_codebook lsp_cb $(CODEBOOKS) > codebook.c + +codebookd.c: generate_codebook $(CODEBOOKSD) + ./generate_codebook lsp_cbd $(CODEBOOKSD) > codebookd.c + +codebookdvq.c: generate_codebook $(CODEBOOKSDVQ) + ./generate_codebook lsp_cbdvq $(CODEBOOKSDVQ) > codebookdvq.c + +clean-local: + -rm -f codebook.c codebookd.c codebookdvq.c + +lib_LTLIBRARIES = libcodec2.la +libcodec2_la_SOURCES = dump.c \ +lpc.c \ +nlp.c \ +postfilter.c \ +sine.c \ +codec2.c \ +fft.c \ +kiss_fft.c \ +interp.c \ +lsp.c \ +phase.c \ +quantise.c \ +pack.c \ +codebook.c \ +codebookd.c \ +codebookdvq.c + +libcodec2_la_CFLAGS = $(AM_CFLAGS) +libcodec2_la_LDFLAGS = $(LIBS) + +generate_codebook_LDFLAGS = $(LIBS) -lm + +library_includedir = $(prefix) +library_include_HEADERS = codec2.h \ +defines.h \ +kiss_fft.h\ +_kiss_fft_guts.h \ +fft.h \ +interp.h \ +lsp.h \ +phase.h \ +quantise.h \ +comp.h \ +dump.h \ +lpc.h \ +nlp.h \ +postfilter.h \ +sine.h diff --git a/gr-codec2-vocoder/src/lib/codec2/_kiss_fft_guts.h b/gr-codec2-vocoder/src/lib/codec2/_kiss_fft_guts.h new file mode 100644 index 000000000..ba6614440 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/_kiss_fft_guts.h @@ -0,0 +1,164 @@ +/* +Copyright (c) 2003-2010, Mark Borgerding + +All rights reserved. + +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 author nor the names of any 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 COPYRIGHT OWNER 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. +*/ + +/* kiss_fft.h + defines kiss_fft_scalar as either short or a float type + and defines + typedef struct { kiss_fft_scalar r; kiss_fft_scalar i; }kiss_fft_cpx; */ +#include "kiss_fft.h" +#include <limits.h> + +#define MAXFACTORS 32 +/* e.g. an fft of length 128 has 4 factors + as far as kissfft is concerned + 4*4*4*2 + */ + +struct kiss_fft_state{ + int nfft; + int inverse; + int factors[2*MAXFACTORS]; + kiss_fft_cpx twiddles[1]; +}; + +/* + Explanation of macros dealing with complex math: + + C_MUL(m,a,b) : m = a*b + C_FIXDIV( c , div ) : if a fixed point impl., c /= div. noop otherwise + C_SUB( res, a,b) : res = a - b + C_SUBFROM( res , a) : res -= a + C_ADDTO( res , a) : res += a + * */ +#ifdef FIXED_POINT +#if (FIXED_POINT==32) +# define FRACBITS 31 +# define SAMPPROD int64_t +#define SAMP_MAX 2147483647 +#else +# define FRACBITS 15 +# define SAMPPROD int32_t +#define SAMP_MAX 32767 +#endif + +#define SAMP_MIN -SAMP_MAX + +#if defined(CHECK_OVERFLOW) +# define CHECK_OVERFLOW_OP(a,op,b) \ + if ( (SAMPPROD)(a) op (SAMPPROD)(b) > SAMP_MAX || (SAMPPROD)(a) op (SAMPPROD)(b) < SAMP_MIN ) { \ + fprintf(stderr,"WARNING:overflow @ " __FILE__ "(%d): (%d " #op" %d) = %ld\n",__LINE__,(a),(b),(SAMPPROD)(a) op (SAMPPROD)(b) ); } +#endif + + +# define smul(a,b) ( (SAMPPROD)(a)*(b) ) +# define sround( x ) (kiss_fft_scalar)( ( (x) + (1<<(FRACBITS-1)) ) >> FRACBITS ) + +# define S_MUL(a,b) sround( smul(a,b) ) + +# define C_MUL(m,a,b) \ + do{ (m).r = sround( smul((a).r,(b).r) - smul((a).i,(b).i) ); \ + (m).i = sround( smul((a).r,(b).i) + smul((a).i,(b).r) ); }while(0) + +# define DIVSCALAR(x,k) \ + (x) = sround( smul( x, SAMP_MAX/k ) ) + +# define C_FIXDIV(c,div) \ + do { DIVSCALAR( (c).r , div); \ + DIVSCALAR( (c).i , div); }while (0) + +# define C_MULBYSCALAR( c, s ) \ + do{ (c).r = sround( smul( (c).r , s ) ) ;\ + (c).i = sround( smul( (c).i , s ) ) ; }while(0) + +#else /* not FIXED_POINT*/ + +# define S_MUL(a,b) ( (a)*(b) ) +#define C_MUL(m,a,b) \ + do{ (m).r = (a).r*(b).r - (a).i*(b).i;\ + (m).i = (a).r*(b).i + (a).i*(b).r; }while(0) +# define C_FIXDIV(c,div) /* NOOP */ +# define C_MULBYSCALAR( c, s ) \ + do{ (c).r *= (s);\ + (c).i *= (s); }while(0) +#endif + +#ifndef CHECK_OVERFLOW_OP +# define CHECK_OVERFLOW_OP(a,op,b) /* noop */ +#endif + +#define C_ADD( res, a,b)\ + do { \ + CHECK_OVERFLOW_OP((a).r,+,(b).r)\ + CHECK_OVERFLOW_OP((a).i,+,(b).i)\ + (res).r=(a).r+(b).r; (res).i=(a).i+(b).i; \ + }while(0) +#define C_SUB( res, a,b)\ + do { \ + CHECK_OVERFLOW_OP((a).r,-,(b).r)\ + CHECK_OVERFLOW_OP((a).i,-,(b).i)\ + (res).r=(a).r-(b).r; (res).i=(a).i-(b).i; \ + }while(0) +#define C_ADDTO( res , a)\ + do { \ + CHECK_OVERFLOW_OP((res).r,+,(a).r)\ + CHECK_OVERFLOW_OP((res).i,+,(a).i)\ + (res).r += (a).r; (res).i += (a).i;\ + }while(0) + +#define C_SUBFROM( res , a)\ + do {\ + CHECK_OVERFLOW_OP((res).r,-,(a).r)\ + CHECK_OVERFLOW_OP((res).i,-,(a).i)\ + (res).r -= (a).r; (res).i -= (a).i; \ + }while(0) + + +#ifdef FIXED_POINT +# define KISS_FFT_COS(phase) floor(.5+SAMP_MAX * cos (phase)) +# define KISS_FFT_SIN(phase) floor(.5+SAMP_MAX * sin (phase)) +# define HALF_OF(x) ((x)>>1) +#elif defined(USE_SIMD) +# define KISS_FFT_COS(phase) _mm_set1_ps( cos(phase) ) +# define KISS_FFT_SIN(phase) _mm_set1_ps( sin(phase) ) +# define HALF_OF(x) ((x)*_mm_set1_ps(.5)) +#else +# define KISS_FFT_COS(phase) (kiss_fft_scalar) cos(phase) +# define KISS_FFT_SIN(phase) (kiss_fft_scalar) sin(phase) +# define HALF_OF(x) ((x)*.5) +#endif + +#define kf_cexp(x,phase) \ + do{ \ + (x)->r = KISS_FFT_COS(phase);\ + (x)->i = KISS_FFT_SIN(phase);\ + }while(0) + + +/* a debugging function */ +#define pcpx(c)\ + fprintf(stderr,"%g + %gi\n",(double)((c)->r),(double)((c)->i) ) + + +#ifdef KISS_FFT_USE_ALLOCA +// define this to allow use of alloca instead of malloc for temporary buffers +// Temporary buffers are used in two case: +// 1. FFT sizes that have "bad" factors. i.e. not 2,3 and 5 +// 2. "in-place" FFTs. Notice the quotes, since kissfft does not really do an in-place transform. +#include <alloca.h> +#define KISS_FFT_TMP_ALLOC(nbytes) alloca(nbytes) +#define KISS_FFT_TMP_FREE(ptr) +#else +#define KISS_FFT_TMP_ALLOC(nbytes) KISS_FFT_MALLOC(nbytes) +#define KISS_FFT_TMP_FREE(ptr) KISS_FFT_FREE(ptr) +#endif diff --git a/gr-codec2-vocoder/src/lib/codec2/c2dec.c b/gr-codec2-vocoder/src/lib/codec2/c2dec.c new file mode 100644 index 000000000..b866d04d6 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/c2dec.c @@ -0,0 +1,82 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: c2dec.c + AUTHOR......: David Rowe + DATE CREATED: 23/8/2010 + + Decodes a file of bits to a file of raw speech samples using codec2. + +\*---------------------------------------------------------------------------*/ + +/* + Copyright (C) 2010 David Rowe + + All rights reserved. + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +#include "codec2.h" + +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <errno.h> + +#define BITS_SIZE ((CODEC2_BITS_PER_FRAME + 7) / 8) + +int main(int argc, char *argv[]) +{ + void *codec2; + FILE *fin; + FILE *fout; + short buf[CODEC2_SAMPLES_PER_FRAME]; + unsigned char bits[BITS_SIZE]; + + if (argc != 3) { + printf("usage: %s InputBitFile OutputRawSpeechFile\n", argv[0]); + exit(1); + } + + if (strcmp(argv[1], "-") == 0) fin = stdin; + else if ( (fin = fopen(argv[1],"rb")) == NULL ) { + fprintf(stderr, "Error opening input bit file: %s: %s.\n", + argv[1], strerror(errno)); + exit(1); + } + + if (strcmp(argv[2], "-") == 0) fout = stdout; + else if ( (fout = fopen(argv[2],"wb")) == NULL ) { + fprintf(stderr, "Error opening output speech file: %s: %s.\n", + argv[2], strerror(errno)); + exit(1); + } + + codec2 = codec2_create(); + + while(fread(bits, sizeof(char), BITS_SIZE, fin) == BITS_SIZE) { + codec2_decode(codec2, buf, bits); + fwrite(buf, sizeof(short), CODEC2_SAMPLES_PER_FRAME, fout); + //if this is in a pipeline, we probably don't want the usual + //buffering to occur + if (fout == stdout) fflush(stdout); + if (fin == stdin) fflush(stdin); + + } + + codec2_destroy(codec2); + + fclose(fin); + fclose(fout); + + return 0; +} diff --git a/gr-codec2-vocoder/src/lib/codec2/c2demo.c b/gr-codec2-vocoder/src/lib/codec2/c2demo.c new file mode 100644 index 000000000..efa8d6449 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/c2demo.c @@ -0,0 +1,86 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: c2demo.c + AUTHOR......: David Rowe + DATE CREATED: 15/11/2010 + + Encodes and decodes a file of raw speech samples using Codec 2. + Demonstrates use of Codec 2 function API. + + Note to convert a wave file to raw and vice-versa: + + $ sox file.wav -r 8000 -s -2 file.raw + $ sox -r 8000 -s -2 file.raw file.wav + +\*---------------------------------------------------------------------------*/ + +/* + Copyright (C) 2010 David Rowe + + All rights reserved. + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +#include "codec2.h" + +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <errno.h> + +#define BITS_SIZE ((CODEC2_BITS_PER_FRAME + 7) / 8) + +int main(int argc, char *argv[]) +{ + void *codec2; + FILE *fin; + FILE *fout; + short buf[CODEC2_SAMPLES_PER_FRAME]; + unsigned char bits[BITS_SIZE]; + + if (argc != 3) { + printf("usage: %s InputRawSpeechFile OutputRawSpeechFile\n", argv[0]); + exit(1); + } + + if ( (fin = fopen(argv[1],"rb")) == NULL ) { + fprintf(stderr, "Error opening input speech file: %s: %s.\n", + argv[1], strerror(errno)); + exit(1); + } + + if ( (fout = fopen(argv[2],"wb")) == NULL ) { + fprintf(stderr, "Error opening output speech file: %s: %s.\n", + argv[2], strerror(errno)); + exit(1); + } + + /* Note only one set of Codec 2 states is required for an encoder + and decoder pair. */ + + codec2 = codec2_create(); + + while(fread(buf, sizeof(short), CODEC2_SAMPLES_PER_FRAME, fin) == + CODEC2_SAMPLES_PER_FRAME) { + codec2_encode(codec2, bits, buf); + codec2_decode(codec2, buf, bits); + fwrite(buf, sizeof(short), CODEC2_SAMPLES_PER_FRAME, fout); + } + + codec2_destroy(codec2); + + fclose(fin); + fclose(fout); + + return 0; +} diff --git a/gr-codec2-vocoder/src/lib/codec2/c2enc.c b/gr-codec2-vocoder/src/lib/codec2/c2enc.c new file mode 100644 index 000000000..4d1d019df --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/c2enc.c @@ -0,0 +1,83 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: c2enc.c + AUTHOR......: David Rowe + DATE CREATED: 23/8/2010 + + Encodes a file of raw speech samples using codec2 and outputs a file + of bits. + +\*---------------------------------------------------------------------------*/ + +/* + Copyright (C) 2010 David Rowe + + All rights reserved. + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +#include "codec2.h" + +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <errno.h> + +#define BITS_SIZE ((CODEC2_BITS_PER_FRAME + 7) / 8) + +int main(int argc, char *argv[]) +{ + void *codec2; + FILE *fin; + FILE *fout; + short buf[CODEC2_SAMPLES_PER_FRAME]; + unsigned char bits[BITS_SIZE]; + + if (argc != 3) { + printf("usage: %s InputRawspeechFile OutputBitFile\n", argv[0]); + exit(1); + } + + if (strcmp(argv[1], "-") == 0) fin = stdin; + else if ( (fin = fopen(argv[1],"rb")) == NULL ) { + fprintf(stderr, "Error opening input bit file: %s: %s.\n", + argv[1], strerror(errno)); + exit(1); + } + + if (strcmp(argv[2], "-") == 0) fout = stdout; + else if ( (fout = fopen(argv[2],"wb")) == NULL ) { + fprintf(stderr, "Error opening output speech file: %s: %s.\n", + argv[2], strerror(errno)); + exit(1); + } + + codec2 = codec2_create(); + + while(fread(buf, sizeof(short), CODEC2_SAMPLES_PER_FRAME, fin) == + CODEC2_SAMPLES_PER_FRAME) { + codec2_encode(codec2, bits, buf); + fwrite(bits, sizeof(char), BITS_SIZE, fout); + //if this is in a pipeline, we probably don't want the usual + //buffering to occur + if (fout == stdout) fflush(stdout); + if (fin == stdin) fflush(stdin); + } + + codec2_destroy(codec2); + + fclose(fin); + fclose(fout); + + return 0; +} diff --git a/gr-codec2-vocoder/src/lib/codec2/c2sim.c b/gr-codec2-vocoder/src/lib/codec2/c2sim.c new file mode 100644 index 000000000..bb49c7899 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/c2sim.c @@ -0,0 +1,469 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: c2sim.c + AUTHOR......: David Rowe + DATE CREATED: 20/8/2010 + + Codec2 simulation. Combines encoder and decoder and allows switching in + out various algorithms and quantisation steps. + +\*---------------------------------------------------------------------------*/ + +/* + 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 Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +#include <assert.h> +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <errno.h> +#include <math.h> + +#include "defines.h" +#include "sine.h" +#include "nlp.h" +#include "dump.h" +#include "lpc.h" +#include "lsp.h" +#include "quantise.h" +#include "phase.h" +#include "postfilter.h" +#include "interp.h" + +/*---------------------------------------------------------------------------*\ + + switch_present() + + Searches the command line arguments for a "switch". If the switch is + found, returns the command line argument where it ws found, else returns + NULL. + +\*---------------------------------------------------------------------------*/ + +int switch_present(sw,argc,argv) +register char sw[]; /* switch in string form */ +register int argc; /* number of command line arguments */ +register char *argv[]; /* array of command line arguments in string form */ +{ + register int i; /* loop variable */ + + for(i=1; i<argc; i++) + if (!strcmp(sw,argv[i])) + return(i); + + return 0; +} + +void synth_one_frame(short buf[], MODEL *model, float Sn_[], float Pn[]); + +/*---------------------------------------------------------------------------*\ + + MAIN + +\*---------------------------------------------------------------------------*/ + +int main(int argc, char *argv[]) +{ + FILE *fout; /* output speech file */ + FILE *fin; /* input speech file */ + short buf[N]; /* input/output buffer */ + float Sn[M]; /* float input speech samples */ + COMP Sw[FFT_ENC]; /* DFT of Sn[] */ + float w[M]; /* time domain hamming window */ + COMP W[FFT_ENC]; /* DFT of w[] */ + MODEL model; + float Pn[2*N]; /* trapezoidal synthesis window */ + float Sn_[2*N]; /* synthesised speech */ + int i; /* loop variable */ + int frames; + float prev_Wo; + float pitch; + int voiced1 = 0; + + char out_file[MAX_STR]; + int arg; + float snr; + float sum_snr; + + int lpc_model, order = LPC_ORD; + int lsp, lspd, lspdvq, lsp_quantiser; + float ak[LPC_MAX]; + COMP Sw_[FFT_ENC]; + COMP Ew[FFT_ENC]; + + int dump; + + int phase0; + float ex_phase[MAX_AMP+1]; + + int postfilt; + float bg_est; + + int hand_voicing; + FILE *fvoicing = 0; + + MODEL prev_model, interp_model; + int decimate; + float lsps[LPC_ORD]; + float prev_lsps[LPC_ORD]; + float e, prev_e; + float ak_interp[LPC_MAX]; + + void *nlp_states; + float hpf_states[2]; + int resample; + float AresdB_prev[MAX_AMP]; + + for(i=0; i<MAX_AMP; i++) + AresdB_prev[i] = 0.0; + + for(i=0; i<M; i++) + Sn[i] = 1.0; + for(i=0; i<2*N; i++) + Sn_[i] = 0; + + prev_Wo = TWO_PI/P_MAX; + + prev_model.Wo = TWO_PI/P_MIN; + prev_model.L = floor(PI/prev_model.Wo); + for(i=1; i<=prev_model.L; i++) { + prev_model.A[i] = 0.0; + prev_model.phi[i] = 0.0; + } + for(i=1; i<=MAX_AMP; i++) { + ex_phase[i] = 0.0; + } + for(i=0; i<LPC_ORD; i++) { + prev_lsps[i] = i*PI/(LPC_ORD+1); + } + e = prev_e = 1; + hpf_states[0] = hpf_states[1] = 0.0; + + nlp_states = nlp_create(); + + if (argc < 2) { + fprintf(stderr, "\nCodec2 - 2400 bit/s speech codec - Simulation Program\n" + "\thttp://rowetel.com/codec2.html\n\n" + "usage: %s InputFile [-o OutputFile]\n" + "\t[--lpc Order]\n" + "\t[--lsp]\n" + "\t[--lspd]\n" + "\t[--lspdvq]\n" + "\t[--phase0]\n" + "\t[--postfilter]\n" + "\t[--hand_voicing]\n" + "\t[--dec]\n" + "\t[--dump DumpFilePrefix]\n", argv[0]); + exit(1); + } + + /* Interpret command line arguments -------------------------------------*/ + + /* Input file */ + + if ((fin = fopen(argv[1],"rb")) == NULL) { + fprintf(stderr, "Error opening input speech file: %s: %s.\n", + argv[1], strerror(errno)); + exit(1); + } + + /* Output file */ + + if ((arg = switch_present("-o",argc,argv))) { + if ((fout = fopen(argv[arg+1],"wb")) == NULL) { + fprintf(stderr, "Error opening output speech file: %s: %s.\n", + argv[arg+1], strerror(errno)); + exit(1); + } + strcpy(out_file,argv[arg+1]); + } + else + fout = NULL; + + lpc_model = 0; + if ((arg = switch_present("--lpc",argc,argv))) { + lpc_model = 1; + order = atoi(argv[arg+1]); + if ((order < 4) || (order > 20)) { + fprintf(stderr, "Error in lpc order: %d\n", order); + exit(1); + } + } + + dump = switch_present("--dump",argc,argv); +#ifdef DUMP + if (dump) + dump_on(argv[dump+1]); +#endif + + lsp = switch_present("--lsp",argc,argv); + lsp_quantiser = 0; + if (lsp) + assert(order == LPC_ORD); + + lspd = switch_present("--lspd",argc,argv); + if (lspd) + assert(order == LPC_ORD); + + lspdvq = switch_present("--lspdvq",argc,argv); + if (lspdvq) + assert(order == LPC_ORD); + + phase0 = switch_present("--phase0",argc,argv); + if (phase0) { + ex_phase[0] = 0; + } + + hand_voicing = switch_present("--hand_voicing",argc,argv); + if (hand_voicing) { + fvoicing = fopen(argv[hand_voicing+1],"rt"); + assert(fvoicing != NULL); + } + + bg_est = 0.0; + postfilt = switch_present("--postfilter",argc,argv); + + decimate = switch_present("--dec",argc,argv); + + arg = switch_present("--resample",argc,argv); + resample = atoi(argv[arg+1]); + + /* Initialise ------------------------------------------------------------*/ + + make_analysis_window(w,W); + make_synthesis_window(Pn); + quantise_init(); + + /* Main loop ------------------------------------------------------------*/ + + frames = 0; + sum_snr = 0; + while(fread(buf,sizeof(short),N,fin)) { + frames++; + //printf("frame: %d", frames); + + /* Read input speech */ + + for(i=0; i<M-N; i++) + Sn[i] = Sn[i+N]; + for(i=0; i<N; i++) { + //Sn[i+M-N] = hpf((float)buf[i], hpf_states); + Sn[i+M-N] = (float)buf[i]; + } + + /* Estimate pitch */ + + nlp(nlp_states,Sn,N,M,P_MIN,P_MAX,&pitch,Sw,&prev_Wo); + model.Wo = TWO_PI/pitch; + + /* estimate model parameters */ + + dft_speech(Sw, Sn, w); + two_stage_pitch_refinement(&model, Sw); + estimate_amplitudes(&model, Sw, W); +#ifdef DUMP + dump_Sn(Sn); dump_Sw(Sw); dump_model(&model); +#endif + + /* optional zero-phase modelling */ + + if (phase0) { + float Wn[M]; /* windowed speech samples */ + float Rk[LPC_MAX+1]; /* autocorrelation coeffs */ + +#ifdef DUMP + dump_phase(&model.phi[0], model.L); +#endif + + /* find aks here, these are overwritten if LPC modelling is enabled */ + + for(i=0; i<M; i++) + Wn[i] = Sn[i]*w[i]; + autocorrelate(Wn,Rk,M,order); + levinson_durbin(Rk,ak,order); + +#ifdef DUMP + dump_ak(ak, LPC_ORD); +#endif + + /* determine voicing */ + + snr = est_voicing_mbe(&model, Sw, W, Sw_, Ew, prev_Wo); +#ifdef DUMP + dump_Sw_(Sw_); + dump_Ew(Ew); + dump_snr(snr); +#endif + + /* just to make sure we are not cheating - kill all phases */ + + for(i=0; i<MAX_AMP; i++) + model.phi[i] = 0; + + if (hand_voicing) { + fscanf(fvoicing,"%d\n",&model.voiced); + } + } + + /* optional LPC model amplitudes */ + + if (lpc_model) { + int lsp_indexes[LPC_MAX]; + + e = speech_to_uq_lsps(lsps, ak, Sn, w, order); + + if (lsp) { + encode_lsps(lsp_indexes, lsps, LPC_ORD); + decode_lsps(lsps, lsp_indexes, LPC_ORD); + bw_expand_lsps(lsps, LPC_ORD); + lsp_to_lpc(lsps, ak, LPC_ORD); + } + + if (lspd) { + float lsps_[LPC_ORD]; + + lspd_quantise(lsps, lsps_, LPC_ORD); + lsp_to_lpc(lsps_, ak, LPC_ORD); + } + + if (lspdvq) { + float lsps_[LPC_ORD]; + + lspdvq_quantise(lsps, lsps_, LPC_ORD); + lsp_to_lpc(lsps_, ak, LPC_ORD); + } + + e = decode_energy(encode_energy(e)); + model.Wo = decode_Wo(encode_Wo(model.Wo)); + + aks_to_M2(ak, order, &model, e, &snr, 1); + apply_lpc_correction(&model); + sum_snr += snr; +#ifdef DUMP + dump_quantised_model(&model); +#endif + } + + /* optional resampling of model amplitudes */ + + printf("frames=%d\n", frames); + if (resample) { + snr = resample_amp_nl(&model, resample, AresdB_prev); + sum_snr += snr; +#ifdef DUMP + dump_quantised_model(&model); +#endif + } + + /* option decimation to 20ms rate, which enables interpolation + routine to synthesise in between frame */ + + if (decimate) { + if (!phase0) { + printf("needs --phase0 to resample phase for interpolated Wo\n"); + exit(0); + } + if (!lpc_model) { + printf("needs --lpc 10 to resample amplitudes\n"); + exit(0); + } + + /* odd frame - interpolate */ + + if (frames%2) { + + interp_model.voiced = voiced1; + + #ifdef LOG_LIN_INTERP + interpolate(&interp_model, &prev_model, &model); + #else + interpolate_lsp(&interp_model, &prev_model, &model, + prev_lsps, prev_e, lsps, e, ak_interp); + apply_lpc_correction(&interp_model); + #endif + + if (phase0) + phase_synth_zero_order(&interp_model, ak_interp, ex_phase, + order); + if (postfilt) + postfilter(&interp_model, &bg_est); + synth_one_frame(buf, &interp_model, Sn_, Pn); + if (fout != NULL) fwrite(buf,sizeof(short),N,fout); + + if (phase0) + phase_synth_zero_order(&model, ak, ex_phase, order); + if (postfilt) + postfilter(&model, &bg_est); + synth_one_frame(buf, &model, Sn_, Pn); + if (fout != NULL) fwrite(buf,sizeof(short),N,fout); + + prev_model = model; + for(i=0; i<LPC_ORD; i++) + prev_lsps[i] = lsps[i]; + prev_e = e; + } + else { + voiced1 = model.voiced; + } + } + else { + if (phase0) + phase_synth_zero_order(&model, ak, ex_phase, order); + if (postfilt) + postfilter(&model, &bg_est); + synth_one_frame(buf, &model, Sn_, Pn); + if (fout != NULL) fwrite(buf,sizeof(short),N,fout); + } + prev_Wo = TWO_PI/pitch; + } + fclose(fin); + + if (fout != NULL) + fclose(fout); + + if (lpc_model || resample) + printf("SNR av = %5.2f dB\n", sum_snr/frames); + +#ifdef DUMP + if (dump) + dump_off(); +#endif + + if (hand_voicing) + fclose(fvoicing); + + nlp_destroy(nlp_states); + + return 0; +} + +void synth_one_frame(short buf[], MODEL *model, float Sn_[], float Pn[]) +{ + int i; + + synthesise(Sn_, model, Pn, 1); + + for(i=0; i<N; i++) { + if (Sn_[i] > 32767.0) + buf[i] = 32767; + else if (Sn_[i] < -32767.0) + buf[i] = -32767; + else + buf[i] = Sn_[i]; + } + +} diff --git a/gr-codec2-vocoder/src/lib/codec2/codebook.c b/gr-codec2-vocoder/src/lib/codec2/codebook.c new file mode 100644 index 000000000..70e1d2abc --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/codebook.c @@ -0,0 +1,245 @@ +/* THIS IS A GENERATED FILE. Edit generate_codebook.c and its input */ + +/* + * This intermediary file and the files that used to create it are under + * The LGPL. See the file COPYING. + */ + +#include "defines.h" + + /* codebook/lsp1.txt */ +static const float codes0[] = { + 225, + 250, + 275, + 300, + 325, + 350, + 375, + 400, + 425, + 450, + 475, + 500, + 525, + 550, + 575, + 600 +}; + /* codebook/lsp2.txt */ +static const float codes1[] = { + 325, + 350, + 375, + 400, + 425, + 450, + 475, + 500, + 525, + 550, + 575, + 600, + 625, + 650, + 675, + 700 +}; + /* codebook/lsp3.txt */ +static const float codes2[] = { + 500, + 550, + 600, + 650, + 700, + 750, + 800, + 850, + 900, + 950, + 1000, + 1050, + 1100, + 1150, + 1200, + 1250 +}; + /* codebook/lsp4.txt */ +static const float codes3[] = { + 700, + 800, + 900, + 1000, + 1100, + 1200, + 1300, + 1400, + 1500, + 1600, + 1700, + 1800, + 1900, + 2000, + 2100, + 2200 +}; + /* codebook/lsp5.txt */ +static const float codes4[] = { + 950, + 1050, + 1150, + 1250, + 1350, + 1450, + 1550, + 1650, + 1750, + 1850, + 1950, + 2050, + 2150, + 2250, + 2350, + 2450 +}; + /* codebook/lsp6.txt */ +static const float codes5[] = { + 1100, + 1200, + 1300, + 1400, + 1500, + 1600, + 1700, + 1800, + 1900, + 2000, + 2100, + 2200, + 2300, + 2400, + 2500, + 2600 +}; + /* codebook/lsp7.txt */ +static const float codes6[] = { + 1500, + 1600, + 1700, + 1800, + 1900, + 2000, + 2100, + 2200, + 2300, + 2400, + 2500, + 2600, + 2700, + 2800, + 2900, + 3000 +}; + /* codebook/lsp8.txt */ +static const float codes7[] = { + 2300, + 2400, + 2500, + 2600, + 2700, + 2800, + 2900, + 3000 +}; + /* codebook/lsp9.txt */ +static const float codes8[] = { + 2500, + 2600, + 2700, + 2800, + 2900, + 3000, + 3100, + 3200 +}; + /* codebook/lsp10.txt */ +static const float codes9[] = { + 2900, + 3100, + 3300, + 3500 +}; + +const struct lsp_codebook lsp_cb[] = { + /* codebook/lsp1.txt */ + { + 1, + 4, + 16, + codes0 + }, + /* codebook/lsp2.txt */ + { + 1, + 4, + 16, + codes1 + }, + /* codebook/lsp3.txt */ + { + 1, + 4, + 16, + codes2 + }, + /* codebook/lsp4.txt */ + { + 1, + 4, + 16, + codes3 + }, + /* codebook/lsp5.txt */ + { + 1, + 4, + 16, + codes4 + }, + /* codebook/lsp6.txt */ + { + 1, + 4, + 16, + codes5 + }, + /* codebook/lsp7.txt */ + { + 1, + 4, + 16, + codes6 + }, + /* codebook/lsp8.txt */ + { + 1, + 3, + 8, + codes7 + }, + /* codebook/lsp9.txt */ + { + 1, + 3, + 8, + codes8 + }, + /* codebook/lsp10.txt */ + { + 1, + 2, + 4, + codes9 + }, + { 0, 0, 0, 0 } +}; diff --git a/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp1.txt b/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp1.txt new file mode 100644 index 000000000..d126be771 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp1.txt @@ -0,0 +1,17 @@ +1 16 +225 +250 +275 +300 +325 +350 +375 +400 +425 +450 +475 +500 +525 +550 +575 +600 diff --git a/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp10.txt b/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp10.txt new file mode 100644 index 000000000..dea9dd9d8 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp10.txt @@ -0,0 +1,9 @@ +1 8 +50 +100 +200 +300 +425 +550 +675 +800 diff --git a/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp2.txt b/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp2.txt new file mode 100644 index 000000000..234bf2067 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp2.txt @@ -0,0 +1,17 @@ +1 16 +25 +50 +75 +100 +125 +150 +175 +200 +225 +250 +275 +300 +325 +350 +375 +400 diff --git a/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp3.txt b/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp3.txt new file mode 100644 index 000000000..b2ee06da4 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp3.txt @@ -0,0 +1,9 @@ +1 8 +50 +75 +100 +120 +150 +250 +350 +450 diff --git a/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp4.txt b/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp4.txt new file mode 100644 index 000000000..dea9dd9d8 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp4.txt @@ -0,0 +1,9 @@ +1 8 +50 +100 +200 +300 +425 +550 +675 +800 diff --git a/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp5.txt b/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp5.txt new file mode 100644 index 000000000..dea9dd9d8 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp5.txt @@ -0,0 +1,9 @@ +1 8 +50 +100 +200 +300 +425 +550 +675 +800 diff --git a/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp6.txt b/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp6.txt new file mode 100644 index 000000000..dea9dd9d8 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp6.txt @@ -0,0 +1,9 @@ +1 8 +50 +100 +200 +300 +425 +550 +675 +800 diff --git a/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp7.txt b/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp7.txt new file mode 100644 index 000000000..dea9dd9d8 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp7.txt @@ -0,0 +1,9 @@ +1 8 +50 +100 +200 +300 +425 +550 +675 +800 diff --git a/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp8.txt b/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp8.txt new file mode 100644 index 000000000..dea9dd9d8 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp8.txt @@ -0,0 +1,9 @@ +1 8 +50 +100 +200 +300 +425 +550 +675 +800 diff --git a/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp9.txt b/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp9.txt new file mode 100644 index 000000000..dea9dd9d8 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/codebook/dlsp9.txt @@ -0,0 +1,9 @@ +1 8 +50 +100 +200 +300 +425 +550 +675 +800 diff --git a/gr-codec2-vocoder/src/lib/codec2/codebook/lsp1.txt b/gr-codec2-vocoder/src/lib/codec2/codebook/lsp1.txt new file mode 100644 index 000000000..d126be771 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/codebook/lsp1.txt @@ -0,0 +1,17 @@ +1 16 +225 +250 +275 +300 +325 +350 +375 +400 +425 +450 +475 +500 +525 +550 +575 +600 diff --git a/gr-codec2-vocoder/src/lib/codec2/codebook/lsp10.txt b/gr-codec2-vocoder/src/lib/codec2/codebook/lsp10.txt new file mode 100644 index 000000000..39aab7c56 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/codebook/lsp10.txt @@ -0,0 +1,6 @@ +1 4 +2900 +3100 +3300 +3500 + diff --git a/gr-codec2-vocoder/src/lib/codec2/codebook/lsp2.txt b/gr-codec2-vocoder/src/lib/codec2/codebook/lsp2.txt new file mode 100644 index 000000000..597f14965 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/codebook/lsp2.txt @@ -0,0 +1,17 @@ +1 16 +325 +350 +375 +400 +425 +450 +475 +500 +525 +550 +575 +600 +625 +650 +675 +700 diff --git a/gr-codec2-vocoder/src/lib/codec2/codebook/lsp3.txt b/gr-codec2-vocoder/src/lib/codec2/codebook/lsp3.txt new file mode 100644 index 000000000..36a64b158 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/codebook/lsp3.txt @@ -0,0 +1,17 @@ +1 16 +500 +550 +600 +650 +700 +750 +800 +850 +900 +950 +1000 +1050 +1100 +1150 +1200 +1250 diff --git a/gr-codec2-vocoder/src/lib/codec2/codebook/lsp4.txt b/gr-codec2-vocoder/src/lib/codec2/codebook/lsp4.txt new file mode 100644 index 000000000..53a90bd8c --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/codebook/lsp4.txt @@ -0,0 +1,17 @@ +1 16 +700 +800 +900 +1000 +1100 +1200 +1300 +1400 +1500 +1600 +1700 +1800 +1900 +2000 +2100 +2200 diff --git a/gr-codec2-vocoder/src/lib/codec2/codebook/lsp5.txt b/gr-codec2-vocoder/src/lib/codec2/codebook/lsp5.txt new file mode 100644 index 000000000..94739b56e --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/codebook/lsp5.txt @@ -0,0 +1,19 @@ +1 16 + 950 +1050 +1150 +1250 +1350 +1450 +1550 +1650 +1750 +1850 +1950 +2050 +2150 +2250 +2350 +2450 + + diff --git a/gr-codec2-vocoder/src/lib/codec2/codebook/lsp6.txt b/gr-codec2-vocoder/src/lib/codec2/codebook/lsp6.txt new file mode 100644 index 000000000..992ea25c5 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/codebook/lsp6.txt @@ -0,0 +1,19 @@ +1 16 +1100 +1200 +1300 +1400 +1500 +1600 +1700 +1800 +1900 +2000 +2100 +2200 +2300 +2400 +2500 +2600 + + diff --git a/gr-codec2-vocoder/src/lib/codec2/codebook/lsp7.txt b/gr-codec2-vocoder/src/lib/codec2/codebook/lsp7.txt new file mode 100644 index 000000000..839cbfdd5 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/codebook/lsp7.txt @@ -0,0 +1,19 @@ +1 16 +1500 +1600 +1700 +1800 +1900 +2000 +2100 +2200 +2300 +2400 +2500 +2600 +2700 +2800 +2900 +3000 + + diff --git a/gr-codec2-vocoder/src/lib/codec2/codebook/lsp8.txt b/gr-codec2-vocoder/src/lib/codec2/codebook/lsp8.txt new file mode 100644 index 000000000..d9880c94e --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/codebook/lsp8.txt @@ -0,0 +1,11 @@ +1 8 +2300 +2400 +2500 +2600 +2700 +2800 +2900 +3000 + + diff --git a/gr-codec2-vocoder/src/lib/codec2/codebook/lsp8910.txt b/gr-codec2-vocoder/src/lib/codec2/codebook/lsp8910.txt new file mode 100644 index 000000000..93cfdd81d --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/codebook/lsp8910.txt @@ -0,0 +1,65 @@ +3 64 +2.048073 2.534502 2.645915 +2.019670 2.269744 2.605462 +1.961101 2.329646 2.562857 +1.968573 2.532712 2.616918 +2.183480 2.514381 2.629582 +2.259379 2.516615 2.620410 +2.172791 2.462460 2.567064 +2.097666 2.303933 2.421685 +2.052990 2.353242 2.546992 +2.043642 2.232362 2.499262 +2.106151 2.393131 2.488401 +2.099167 2.437862 2.558655 +2.013877 2.422875 2.530071 +2.033848 2.483776 2.584598 +2.114474 2.516856 2.602372 +2.229214 2.584056 2.678855 +2.131151 2.584299 2.674845 +1.472721 2.477091 2.630241 +2.010907 2.598415 2.682989 +2.353653 2.524066 2.619773 +2.419897 2.623938 2.699605 +2.319080 2.602148 2.689044 +1.860342 2.503881 2.616576 +1.910517 2.386693 2.610126 +1.748689 2.371809 2.496542 +1.618495 2.403425 2.554956 +1.844073 2.437026 2.533443 +1.924810 2.388543 2.502698 +1.937227 2.258363 2.501697 +1.687554 2.209123 2.545239 +1.851950 2.278628 2.565632 +1.868154 2.330150 2.444883 +1.874180 2.213118 2.351940 +1.757311 2.030626 2.433836 +1.650306 2.152371 2.243421 +1.612794 1.884686 2.339313 +1.745431 2.278895 2.389449 +1.590923 2.304155 2.408510 +1.475982 2.275548 2.509897 +1.508695 2.045463 2.455520 +1.872054 2.061777 2.246202 +1.983947 2.159155 2.445535 +1.745180 2.483765 2.593698 +1.900116 2.079600 2.407479 +1.841672 2.167042 2.486827 +1.932912 2.148464 2.569850 +2.134174 2.363673 2.584252 +2.106094 2.450645 2.638417 +1.954135 2.460313 2.666512 +1.907634 2.573801 2.674025 +1.625579 2.539569 2.656363 +1.785866 2.572616 2.676082 +1.798447 2.376454 2.624298 +2.020033 2.397244 2.619868 +1.946581 2.468791 2.564185 +2.008920 2.342400 2.469132 +1.983846 2.271044 2.395408 +1.988039 2.154150 2.317920 +2.077197 2.216622 2.389101 +2.117255 2.283907 2.512242 +2.177233 2.334622 2.458268 +2.214655 2.425510 2.620013 +2.199931 2.390272 2.520731 +2.271755 2.448682 2.552649 diff --git a/gr-codec2-vocoder/src/lib/codec2/codebook/lsp9.txt b/gr-codec2-vocoder/src/lib/codec2/codebook/lsp9.txt new file mode 100644 index 000000000..7e159af2f --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/codebook/lsp9.txt @@ -0,0 +1,11 @@ +1 8 +2500 +2600 +2700 +2800 +2900 +3000 +3100 +3200 + + diff --git a/gr-codec2-vocoder/src/lib/codec2/codebookd.c b/gr-codec2-vocoder/src/lib/codec2/codebookd.c new file mode 100644 index 000000000..fbe12e5ed --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/codebookd.c @@ -0,0 +1,209 @@ +/* THIS IS A GENERATED FILE. Edit generate_codebook.c and its input */ + +/* + * This intermediary file and the files that used to create it are under + * The LGPL. See the file COPYING. + */ + +#include "defines.h" + + /* codebook/dlsp1.txt */ +static const float codes0[] = { + 225, + 250, + 275, + 300, + 325, + 350, + 375, + 400, + 425, + 450, + 475, + 500, + 525, + 550, + 575, + 600 +}; + /* codebook/dlsp2.txt */ +static const float codes1[] = { + 25, + 50, + 75, + 100, + 125, + 150, + 175, + 200, + 225, + 250, + 275, + 300, + 325, + 350, + 375, + 400 +}; + /* codebook/dlsp3.txt */ +static const float codes2[] = { + 50, + 75, + 100, + 120, + 150, + 250, + 350, + 450 +}; + /* codebook/dlsp4.txt */ +static const float codes3[] = { + 50, + 100, + 200, + 300, + 425, + 550, + 675, + 800 +}; + /* codebook/dlsp5.txt */ +static const float codes4[] = { + 50, + 100, + 200, + 300, + 425, + 550, + 675, + 800 +}; + /* codebook/dlsp6.txt */ +static const float codes5[] = { + 50, + 100, + 200, + 300, + 425, + 550, + 675, + 800 +}; + /* codebook/dlsp7.txt */ +static const float codes6[] = { + 50, + 100, + 200, + 300, + 425, + 550, + 675, + 800 +}; + /* codebook/dlsp8.txt */ +static const float codes7[] = { + 50, + 100, + 200, + 300, + 425, + 550, + 675, + 800 +}; + /* codebook/dlsp9.txt */ +static const float codes8[] = { + 50, + 100, + 200, + 300, + 425, + 550, + 675, + 800 +}; + /* codebook/dlsp10.txt */ +static const float codes9[] = { + 50, + 100, + 200, + 300, + 425, + 550, + 675, + 800 +}; + +const struct lsp_codebook lsp_cbd[] = { + /* codebook/dlsp1.txt */ + { + 1, + 4, + 16, + codes0 + }, + /* codebook/dlsp2.txt */ + { + 1, + 4, + 16, + codes1 + }, + /* codebook/dlsp3.txt */ + { + 1, + 3, + 8, + codes2 + }, + /* codebook/dlsp4.txt */ + { + 1, + 3, + 8, + codes3 + }, + /* codebook/dlsp5.txt */ + { + 1, + 3, + 8, + codes4 + }, + /* codebook/dlsp6.txt */ + { + 1, + 3, + 8, + codes5 + }, + /* codebook/dlsp7.txt */ + { + 1, + 3, + 8, + codes6 + }, + /* codebook/dlsp8.txt */ + { + 1, + 3, + 8, + codes7 + }, + /* codebook/dlsp9.txt */ + { + 1, + 3, + 8, + codes8 + }, + /* codebook/dlsp10.txt */ + { + 1, + 3, + 8, + codes9 + }, + { 0, 0, 0, 0 } +}; diff --git a/gr-codec2-vocoder/src/lib/codec2/codebookdvq.c b/gr-codec2-vocoder/src/lib/codec2/codebookdvq.c new file mode 100644 index 000000000..63cd373e8 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/codebookdvq.c @@ -0,0 +1,119 @@ +/* THIS IS A GENERATED FILE. Edit generate_codebook.c and its input */ + +/* + * This intermediary file and the files that used to create it are under + * The LGPL. See the file COPYING. + */ + +#include "defines.h" + + /* codebook/dlsp1.txt */ +static const float codes0[] = { + 225, + 250, + 275, + 300, + 325, + 350, + 375, + 400, + 425, + 450, + 475, + 500, + 525, + 550, + 575, + 600 +}; + /* codebook/dlsp2.txt */ +static const float codes1[] = { + 25, + 50, + 75, + 100, + 125, + 150, + 175, + 200, + 225, + 250, + 275, + 300, + 325, + 350, + 375, + 400 +}; + /* codebook/dlsp3.txt */ +static const float codes2[] = { + 50, + 75, + 100, + 120, + 150, + 250, + 350, + 450 +}; + /* codebook/dlsp4.txt */ +static const float codes3[] = { + 50, + 100, + 200, + 300, + 425, + 550, + 675, + 800 +}; + /* codebook/dlsp5.txt */ +static const float codes4[] = { + 50, + 100, + 200, + 300, + 425, + 550, + 675, + 800 +}; + +const struct lsp_codebook lsp_cbdvq[] = { + /* codebook/dlsp1.txt */ + { + 1, + 4, + 16, + codes0 + }, + /* codebook/dlsp2.txt */ + { + 1, + 4, + 16, + codes1 + }, + /* codebook/dlsp3.txt */ + { + 1, + 3, + 8, + codes2 + }, + /* codebook/dlsp4.txt */ + { + 1, + 3, + 8, + codes3 + }, + /* codebook/dlsp5.txt */ + { + 1, + 3, + 8, + codes4 + }, + { 0, 0, 0, 0 } +}; diff --git a/gr-codec2-vocoder/src/lib/codec2/codec2.c b/gr-codec2-vocoder/src/lib/codec2/codec2.c new file mode 100644 index 000000000..92708ee32 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/codec2.c @@ -0,0 +1,342 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: codec2.c + AUTHOR......: David Rowe + DATE CREATED: 21/8/2010 + + Codec2 fully quantised encoder and decoder functions. If you want use + codec2, the codec2_xxx functions are for you. + +\*---------------------------------------------------------------------------*/ + +/* + Copyright (C) 2010 David Rowe + + All rights reserved. + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +#include <assert.h> +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <math.h> + +#include "defines.h" +#include "sine.h" +#include "nlp.h" +#include "dump.h" +#include "lpc.h" +#include "quantise.h" +#include "phase.h" +#include "interp.h" +#include "postfilter.h" +#include "codec2.h" +#include "codec2_internal.h" + +/*---------------------------------------------------------------------------*\ + + FUNCTIONS + +\*---------------------------------------------------------------------------*/ + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: codec2_create + AUTHOR......: David Rowe + DATE CREATED: 21/8/2010 + + Create and initialise an instance of the codec. Returns a pointer + to the codec states or NULL on failure. One set of states is + sufficient for a full duuplex codec (i.e. an encoder and decoder). + You don't need separate states for encoders and decoders. See + c2enc.c and c2dec.c for examples. + +\*---------------------------------------------------------------------------*/ + +void *codec2_create() +{ + CODEC2 *c2; + int i,l; + + c2 = (CODEC2*)malloc(sizeof(CODEC2)); + if (c2 == NULL) + return NULL; + + for(i=0; i<M; i++) + c2->Sn[i] = 1.0; + c2->hpf_states[0] = c2->hpf_states[1] = 0.0; + for(i=0; i<2*N; i++) + c2->Sn_[i] = 0; + make_analysis_window(c2->w,c2->W); + make_synthesis_window(c2->Pn); + quantise_init(); + c2->prev_Wo = 0.0; + c2->bg_est = 0.0; + c2->ex_phase = 0.0; + + for(l=1; l<MAX_AMP; l++) + c2->prev_model.A[l] = 0.0; + c2->prev_model.Wo = TWO_PI/P_MAX; + c2->prev_model.L = PI/c2->prev_model.Wo; + c2->prev_model.voiced = 0; + + for(i=0; i<LPC_ORD; i++) { + c2->prev_lsps[i] = i*PI/(LPC_ORD+1); + } + c2->prev_energy = 1; + + c2->nlp = nlp_create(); + if (c2->nlp == NULL) { + free (c2); + return NULL; + } + + return (void*)c2; +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: codec2_create + AUTHOR......: David Rowe + DATE CREATED: 21/8/2010 + + Destroy an instance of the codec. + +\*---------------------------------------------------------------------------*/ + +void codec2_destroy(void *codec2_state) +{ + CODEC2 *c2; + + assert(codec2_state != NULL); + c2 = (CODEC2*)codec2_state; + nlp_destroy(c2->nlp); + free(codec2_state); +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: codec2_encode + AUTHOR......: David Rowe + DATE CREATED: 21/8/2010 + + Encodes 160 speech samples (20ms of speech) into 51 bits. + + The codec2 algorithm actually operates internally on 10ms (80 + sample) frames, so we run the encoding algorithm twice. On the + first frame we just send the voicing bit. One the second frame we + send all model parameters. + + The bit allocation is: + + Parameter bits/frame + -------------------------------------- + Harmonic magnitudes (LSPs) 36 + Low frequency LPC correction 1 + Energy 5 + Wo (fundamental frequnecy) 7 + Voicing (10ms update) 2 + TOTAL 51 + +\*---------------------------------------------------------------------------*/ + +void codec2_encode(void *codec2_state, unsigned char * bits, short speech[]) +{ + CODEC2 *c2; + MODEL model; + int voiced1, voiced2; + int lsp_indexes[LPC_ORD]; + int energy_index; + int Wo_index; + int i; + unsigned int nbit = 0; + + assert(codec2_state != NULL); + c2 = (CODEC2*)codec2_state; + + /* first 10ms analysis frame - we just want voicing */ + + analyse_one_frame(c2, &model, speech); + voiced1 = model.voiced; + + /* second 10ms analysis frame */ + + analyse_one_frame(c2, &model, &speech[N]); + voiced2 = model.voiced; + + Wo_index = encode_Wo(model.Wo); + encode_amplitudes(lsp_indexes, + &energy_index, + &model, + c2->Sn, + c2->w); + memset(bits, '\0', ((CODEC2_BITS_PER_FRAME + 7) / 8)); + pack(bits, &nbit, Wo_index, WO_BITS); + for(i=0; i<LPC_ORD; i++) { + pack(bits, &nbit, lsp_indexes[i], lsp_bits(i)); + } + pack(bits, &nbit, energy_index, E_BITS); + pack(bits, &nbit, voiced1, 1); + pack(bits, &nbit, voiced2, 1); + + assert(nbit == CODEC2_BITS_PER_FRAME); +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: codec2_decode + AUTHOR......: David Rowe + DATE CREATED: 21/8/2010 + + Decodes frames of 51 bits into 160 samples (20ms) of speech. + +\*---------------------------------------------------------------------------*/ + +void codec2_decode(void *codec2_state, short speech[], + const unsigned char * bits) +{ + CODEC2 *c2; + MODEL model; + int voiced1, voiced2; + int lsp_indexes[LPC_ORD]; + float lsps[LPC_ORD]; + int energy_index; + float energy; + int Wo_index; + float ak[LPC_ORD+1]; + float ak_interp[LPC_ORD+1]; + int i; + unsigned int nbit = 0; + MODEL model_interp; + + assert(codec2_state != NULL); + c2 = (CODEC2*)codec2_state; + + /* unpack bit stream to integer codes */ + + Wo_index = unpack(bits, &nbit, WO_BITS); + for(i=0; i<LPC_ORD; i++) { + lsp_indexes[i] = unpack(bits, &nbit, lsp_bits(i)); + } + energy_index = unpack(bits, &nbit, E_BITS); + voiced1 = unpack(bits, &nbit, 1); + voiced2 = unpack(bits, &nbit, 1); + assert(nbit == CODEC2_BITS_PER_FRAME); + + /* decode integer codes to model parameters */ + + model.Wo = decode_Wo(Wo_index); + model.L = PI/model.Wo; + memset(&model.A, 0, (model.L+1)*sizeof(model.A[0])); + decode_amplitudes(&model, + ak, + lsp_indexes, + energy_index, + lsps, + &energy); + + model.voiced = voiced2; + model_interp.voiced = voiced1; + model_interp.Wo = P_MAX/2; + memset(&model_interp.A, 0, MAX_AMP*sizeof(model_interp.A[0])); + + /* interpolate middle frame's model parameters for adjacent frames */ + + interpolate_lsp(&model_interp, &c2->prev_model, &model, + c2->prev_lsps, c2->prev_energy, lsps, energy, ak_interp); + apply_lpc_correction(&model_interp); + + /* synthesis two 10ms frames */ + + synthesise_one_frame(c2, speech, &model_interp, ak_interp); + synthesise_one_frame(c2, &speech[N], &model, ak); + + /* update memories (decode states) for next time */ + + memcpy(&c2->prev_model, &model, sizeof(MODEL)); + memcpy(c2->prev_lsps, lsps, sizeof(lsps)); + c2->prev_energy = energy; +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: synthesise_one_frame() + AUTHOR......: David Rowe + DATE CREATED: 23/8/2010 + + Synthesise 80 speech samples (10ms) from model parameters. + +\*---------------------------------------------------------------------------*/ + +void synthesise_one_frame(CODEC2 *c2, short speech[], MODEL *model, float ak[]) +{ + int i; + + phase_synth_zero_order(model, ak, &c2->ex_phase, LPC_ORD); + postfilter(model, &c2->bg_est); + synthesise(c2->Sn_, model, c2->Pn, 1); + + for(i=0; i<N; i++) { + if (c2->Sn_[i] > 32767.0) + speech[i] = 32767; + else if (c2->Sn_[i] < -32767.0) + speech[i] = -32767; + else + speech[i] = c2->Sn_[i]; + } + +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: analyse_one_frame() + AUTHOR......: David Rowe + DATE CREATED: 23/8/2010 + + Extract sinusoidal model parameters from 80 speech samples (10ms of + speech). + +\*---------------------------------------------------------------------------*/ + +void analyse_one_frame(CODEC2 *c2, MODEL *model, short speech[]) +{ + COMP Sw[FFT_ENC]; + COMP Sw_[FFT_ENC]; + COMP Ew[FFT_ENC]; + float pitch; + int i; + + /* Read input speech */ + + for(i=0; i<M-N; i++) + c2->Sn[i] = c2->Sn[i+N]; + for(i=0; i<N; i++) + c2->Sn[i+M-N] = speech[i]; + + dft_speech(Sw, c2->Sn, c2->w); + + /* Estimate pitch */ + + nlp(c2->nlp,c2->Sn,N,M,P_MIN,P_MAX,&pitch,Sw,&c2->prev_Wo); + model->Wo = TWO_PI/pitch; + model->L = PI/model->Wo; + + /* estimate model parameters */ + + two_stage_pitch_refinement(model, Sw); + estimate_amplitudes(model, Sw, c2->W); + est_voicing_mbe(model, Sw, c2->W, Sw_, Ew, c2->prev_Wo); + + c2->prev_Wo = model->Wo; +} diff --git a/gr-codec2-vocoder/src/lib/codec2/codec2.h b/gr-codec2-vocoder/src/lib/codec2/codec2.h new file mode 100644 index 000000000..946dedca5 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/codec2.h @@ -0,0 +1,41 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: codec2.h + AUTHOR......: David Rowe + DATE CREATED: 21/8/2010 + + Codec2 fully quantised encoder and decoder functions. If you want use + codec2, these are the functions you need to call. + +\*---------------------------------------------------------------------------*/ + +/* + Copyright (C) 2010 David Rowe + + All rights reserved. + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +#ifndef __CODEC2__ +#define __CODEC2__ + +#define CODEC2_SAMPLES_PER_FRAME 160 +#define CODEC2_BITS_PER_FRAME 50 + +void *codec2_create(); +void codec2_destroy(void *codec2_state); +void codec2_encode(void *codec2_state, unsigned char * bits, short speech_in[]); +void codec2_decode(void *codec2_state, short speech_out[], + const unsigned char * bits); + +#endif diff --git a/gr-codec2-vocoder/src/lib/codec2/codec2_internal.h b/gr-codec2-vocoder/src/lib/codec2/codec2_internal.h new file mode 100644 index 000000000..3943ac29d --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/codec2_internal.h @@ -0,0 +1,63 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: codec2_internal.h + AUTHOR......: David Rowe + DATE CREATED: 22 March 2011 + + Some internal structures and states broken out here as they are useful for + testing and development. + +\*---------------------------------------------------------------------------*/ + +/* + Copyright (C) 2011 David Rowe + + All rights reserved. + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +#ifndef __CODEC2_INTERNAL__ +#define __CODEC2_INTERNAL__ + +/*---------------------------------------------------------------------------*\ + + STATES + +\*---------------------------------------------------------------------------*/ + +typedef struct { + float w[M]; /* time domain hamming window */ + COMP W[FFT_ENC]; /* DFT of w[] */ + float Pn[2*N]; /* trapezoidal synthesis window */ + float Sn[M]; /* input speech */ + float hpf_states[2]; /* high pass filter states */ + void *nlp; /* pitch predictor states */ + float Sn_[2*N]; /* synthesised output speech */ + float ex_phase; /* excitation model phase track */ + float bg_est; /* background noise estimate for post filter */ + float prev_Wo; /* previous frame's pitch estimate */ + MODEL prev_model; /* previous frame's model parameters */ + float prev_lsps[LPC_ORD]; /* previous frame's LSPs */ + float prev_energy; /* previous frame's LPC energy */ +} CODEC2; + +/*---------------------------------------------------------------------------*\ + + FUNCTION HEADERS + +\*---------------------------------------------------------------------------*/ + +void analyse_one_frame(CODEC2 *c2, MODEL *model, short speech[]); +void synthesise_one_frame(CODEC2 *c2, short speech[], MODEL *model,float ak[]); + +#endif diff --git a/gr-codec2-vocoder/src/lib/codec2/comp.h b/gr-codec2-vocoder/src/lib/codec2/comp.h new file mode 100644 index 000000000..cedcab37f --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/comp.h @@ -0,0 +1,38 @@ +/*---------------------------------------------------------------------------*\ + + 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 Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +#ifndef __COMP__ +#define __COMP__ + +/* Complex number */ + +typedef struct { + float real; + float imag; +} COMP; + +#endif diff --git a/gr-codec2-vocoder/src/lib/codec2/defines.h b/gr-codec2-vocoder/src/lib/codec2/defines.h new file mode 100644 index 000000000..2dcd527d3 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/defines.h @@ -0,0 +1,88 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: defines.h + AUTHOR......: David Rowe + DATE CREATED: 23/4/93 + + Defines and structures used throughout the codec. + +\*---------------------------------------------------------------------------*/ + +/* + 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 Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +#ifndef __DEFINES__ +#define __DEFINES__ + +/*---------------------------------------------------------------------------*\ + + DEFINES + +\*---------------------------------------------------------------------------*/ + +/* General defines */ + +#define N 80 /* number of samples per frame */ +#define MAX_AMP 80 /* maximum number of harmonics */ +#define PI 3.141592654 /* mathematical constant */ +#define TWO_PI 6.283185307 /* mathematical constant */ +#define FS 8000 /* sample rate in Hz */ +#define MAX_STR 256 /* maximum string size */ + +#define NW 279 /* analysis window size */ +#define FFT_ENC 512 /* size of FFT used for encoder */ +#define FFT_DEC 512 /* size of FFT used in decoder */ +#define TW 40 /* Trapezoidal synthesis window overlap */ +#define V_THRESH 6.0 /* voicing threshold in dB */ +#define LPC_MAX 20 /* maximum LPC order */ +#define LPC_ORD 10 /* phase modelling LPC order */ + +/* Pitch estimation defines */ + +#define M 320 /* pitch analysis frame size */ +#define P_MIN 20 /* minimum pitch */ +#define P_MAX 160 /* maximum pitch */ + +/*---------------------------------------------------------------------------*\ + + TYPEDEFS + +\*---------------------------------------------------------------------------*/ + +/* Structure to hold model parameters for one frame */ + +typedef struct { + float Wo; /* fundamental frequency estimate in radians */ + int L; /* number of harmonics */ + float A[MAX_AMP]; /* amplitiude of each harmonic */ + float phi[MAX_AMP]; /* phase of each harmonic */ + int voiced; /* non-zero if this frame is voiced */ +} MODEL; + +/* describes each codebook */ + +struct lsp_codebook { + int k; /* dimension of vector */ + int log2m; /* number of bits in m */ + int m; /* elements in codebook */ + const float * cb; /* The elements */ +}; +extern const struct lsp_codebook lsp_cb[]; +extern const struct lsp_codebook lsp_cbd[]; +extern const struct lsp_codebook lsp_cbdvq[]; + +#endif diff --git a/gr-codec2-vocoder/src/lib/codec2/dump.c b/gr-codec2-vocoder/src/lib/codec2/dump.c new file mode 100644 index 000000000..73a378e23 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/dump.c @@ -0,0 +1,469 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: dump.c + AUTHOR......: David Rowe + DATE CREATED: 25/8/09 + + Routines to dump data to text files for Octave analysis. + +\*---------------------------------------------------------------------------*/ + +/* + All rights reserved. + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +#include "defines.h" +#include "comp.h" +#include "dump.h" +#include <assert.h> +#include <stdlib.h> +#include <stdio.h> +#include <string.h> +#include <math.h> + +#ifdef DUMP +static int dumpon = 0; + +static FILE *fsn = NULL; +static FILE *fsw = NULL; +static FILE *few = NULL; +static FILE *fsw_ = NULL; +static FILE *fmodel = NULL; +static FILE *fqmodel = NULL; +static FILE *fpw = NULL; +static FILE *flsp = NULL; +static FILE *fphase = NULL; +static FILE *fphase_ = NULL; +static FILE *ffw = NULL; +static FILE *fe = NULL; +static FILE *fsq = NULL; +static FILE *fdec = NULL; +static FILE *fsnr = NULL; +static FILE *fak = NULL; +static FILE *fbg = NULL; +static FILE *fE = NULL; +static FILE *frk = NULL; +static FILE *fres = NULL; + +static char prefix[MAX_STR]; + +void dump_on(char p[]) { + dumpon = 1; + strcpy(prefix, p); +} + +void dump_off(){ + if (fsn != NULL) + fclose(fsn); + if (fsw != NULL) + fclose(fsw); + if (fsw_ != NULL) + fclose(fsw_); + if (few != NULL) + fclose(few); + if (fmodel != NULL) + fclose(fmodel); + if (fqmodel != NULL) + fclose(fqmodel); + if (fpw != NULL) + fclose(fpw); + if (flsp != NULL) + fclose(flsp); + if (fphase != NULL) + fclose(fphase); + if (fphase_ != NULL) + fclose(fphase_); + if (ffw != NULL) + fclose(ffw); + if (fe != NULL) + fclose(fe); + if (fsq != NULL) + fclose(fsq); + if (fdec != NULL) + fclose(fdec); + if (fsnr != NULL) + fclose(fsnr); + if (fak != NULL) + fclose(fak); + if (fbg != NULL) + fclose(fbg); + if (fE != NULL) + fclose(fE); + if (frk != NULL) + fclose(frk); + if (fres != NULL) + fclose(fres); +} + +void dump_Sn(float Sn[]) { + int i; + char s[MAX_STR]; + + if (!dumpon) return; + + if (fsn == NULL) { + sprintf(s,"%s_sn.txt", prefix); + fsn = fopen(s, "wt"); + assert(fsn != NULL); + } + + /* split across two lines to avoid max line length problems */ + /* reconstruct in Octave */ + + for(i=0; i<M/2; i++) + fprintf(fsn,"%f\t",Sn[i]); + fprintf(fsn,"\n"); + for(i=M/2; i<M; i++) + fprintf(fsn,"%f\t",Sn[i]); + fprintf(fsn,"\n"); +} + +void dump_Sw(COMP Sw[]) { + int i; + char s[MAX_STR]; + + if (!dumpon) return; + + if (fsw == NULL) { + sprintf(s,"%s_sw.txt", prefix); + fsw = fopen(s, "wt"); + assert(fsw != NULL); + } + + for(i=0; i<FFT_ENC/2; i++) + fprintf(fsw,"%f\t", + 10.0*log10(Sw[i].real*Sw[i].real + Sw[i].imag*Sw[i].imag)); + fprintf(fsw,"\n"); +} + +void dump_Sw_(COMP Sw_[]) { + int i; + char s[MAX_STR]; + + if (!dumpon) return; + + if (fsw_ == NULL) { + sprintf(s,"%s_sw_.txt", prefix); + fsw_ = fopen(s, "wt"); + assert(fsw_ != NULL); + } + + for(i=0; i<FFT_ENC/2; i++) + fprintf(fsw_,"%f\t", + 10.0*log10(Sw_[i].real*Sw_[i].real + Sw_[i].imag*Sw_[i].imag)); + fprintf(fsw_,"\n"); +} + +void dump_Ew(COMP Ew[]) { + int i; + char s[MAX_STR]; + + if (!dumpon) return; + + if (few == NULL) { + sprintf(s,"%s_ew.txt", prefix); + few = fopen(s, "wt"); + assert(few != NULL); + } + + for(i=0; i<FFT_ENC/2; i++) + fprintf(few,"%f\t", + 10.0*log10(Ew[i].real*Ew[i].real + Ew[i].imag*Ew[i].imag)); + fprintf(few,"\n"); +} + +void dump_model(MODEL *model) { + int l; + char s[MAX_STR]; + + if (!dumpon) return; + + if (fmodel == NULL) { + sprintf(s,"%s_model.txt", prefix); + fmodel = fopen(s, "wt"); + assert(fmodel != NULL); + } + + fprintf(fmodel,"%f\t%d\t", model->Wo, model->L); + for(l=1; l<=model->L; l++) + fprintf(fmodel,"%f\t",model->A[l]); + for(l=model->L+1; l<MAX_AMP; l++) + fprintf(fmodel,"0.0\t"); + fprintf(fmodel,"%d\t",model->voiced); + fprintf(fmodel,"\n"); +} + +void dump_quantised_model(MODEL *model) { + int l; + char s[MAX_STR]; + + if (!dumpon) return; + + if (fqmodel == NULL) { + sprintf(s,"%s_qmodel.txt", prefix); + fqmodel = fopen(s, "wt"); + assert(fqmodel != NULL); + } + + fprintf(fqmodel,"%f\t%d\t", model->Wo, model->L); + for(l=1; l<=model->L; l++) + fprintf(fqmodel,"%f\t",model->A[l]); + for(l=model->L+1; l<MAX_AMP; l++) + fprintf(fqmodel,"0.0\t"); + fprintf(fqmodel,"\n"); +} + +void dump_resample(float w[], float A[], int n) { + int l; + char s[MAX_STR]; + + if (!dumpon) return; + + if (fres == NULL) { + sprintf(s,"%s_res.txt", prefix); + fres = fopen(s, "wt"); + assert(fres != NULL); + } + + fprintf(fres,"%d\t",n); + for(l=0; l<n; l++) + fprintf(fres,"%f\t",w[l]); + for(l=0; l<n; l++) + fprintf(fres,"%f\t",A[l]); + fprintf(fres,"\n"); +} + +void dump_phase(float phase[], int L) { + int l; + char s[MAX_STR]; + + if (!dumpon) return; + + if (fphase == NULL) { + sprintf(s,"%s_phase.txt", prefix); + fphase = fopen(s, "wt"); + assert(fphase != NULL); + } + + for(l=1; l<=L; l++) + fprintf(fphase,"%f\t",phase[l]); + for(l=L+1; l<MAX_AMP; l++) + fprintf(fphase,"%f\t",0.0); + fprintf(fphase,"\n"); +} + +void dump_phase_(float phase_[], int L) { + int l; + char s[MAX_STR]; + + if (!dumpon) return; + + if (fphase_ == NULL) { + sprintf(s,"%s_phase_.txt", prefix); + fphase_ = fopen(s, "wt"); + assert(fphase_ != NULL); + } + + for(l=1; l<=L; l++) + fprintf(fphase_,"%f\t",phase_[l]); + for(l=L+1; l<MAX_AMP; l++) + fprintf(fphase_,"%f\t",0.0); + fprintf(fphase_,"\n"); +} + +void dump_snr(float snr) { + char s[MAX_STR]; + + if (!dumpon) return; + + if (fsnr == NULL) { + sprintf(s,"%s_snr.txt", prefix); + fsnr = fopen(s, "wt"); + assert(fsnr != NULL); + } + + fprintf(fsnr,"%f\n",snr); +} + +void dump_Pw(COMP Pw[]) { + int i; + char s[MAX_STR]; + + if (!dumpon) return; + + if (fpw == NULL) { + sprintf(s,"%s_pw.txt", prefix); + fpw = fopen(s, "wt"); + assert(fpw != NULL); + } + + for(i=0; i<FFT_DEC/2; i++) + fprintf(fpw,"%f\t",Pw[i].real); + fprintf(fpw,"\n"); +} + +void dump_lsp(float lsp[]) { + int i; + char s[MAX_STR]; + + if (!dumpon) return; + + if (flsp == NULL) { + sprintf(s,"%s_lsp.txt", prefix); + flsp = fopen(s, "wt"); + assert(flsp != NULL); + } + + for(i=0; i<10; i++) + fprintf(flsp,"%f\t",lsp[i]); + fprintf(flsp,"\n"); +} + +void dump_ak(float ak[], int order) { + int i; + char s[MAX_STR]; + + if (!dumpon) return; + + if (fak == NULL) { + sprintf(s,"%s_ak.txt", prefix); + fak = fopen(s, "wt"); + assert(fak != NULL); + } + + for(i=0; i<=order; i++) + fprintf(fak,"%f\t",ak[i]); + fprintf(fak,"\n"); +} + +void dump_Fw(COMP Fw[]) { + int i; + char s[MAX_STR]; + + if (!dumpon) return; + + if (ffw == NULL) { + sprintf(s,"%s_fw.txt", prefix); + ffw = fopen(s, "wt"); + assert(ffw != NULL); + } + + for(i=0; i<256; i++) + fprintf(ffw,"%f\t",Fw[i].real); + fprintf(ffw,"\n"); +} + +void dump_e(float e_hz[]) { + int i; + char s[MAX_STR]; + + if (!dumpon) return; + + if (fe == NULL) { + sprintf(s,"%s_e.txt", prefix); + fe = fopen(s, "wt"); + assert(fe != NULL); + } + + for(i=0; i<500/2; i++) + fprintf(fe,"%f\t",e_hz[i]); + fprintf(fe,"\n"); + for(i=500/2; i<500; i++) + fprintf(fe,"%f\t",e_hz[i]); + fprintf(fe,"\n"); +} + +void dump_sq(float sq[]) { + int i; + char s[MAX_STR]; + + if (!dumpon) return; + + if (fsq == NULL) { + sprintf(s,"%s_sq.txt", prefix); + fsq = fopen(s, "wt"); + assert(fsq != NULL); + } + + for(i=0; i<M/2; i++) + fprintf(fsq,"%f\t",sq[i]); + fprintf(fsq,"\n"); + for(i=M/2; i<M; i++) + fprintf(fsq,"%f\t",sq[i]); + fprintf(fsq,"\n"); +} + +void dump_dec(COMP Fw[]) { + int i; + char s[MAX_STR]; + + if (!dumpon) return; + + if (fdec == NULL) { + sprintf(s,"%s_dec.txt", prefix); + fdec = fopen(s, "wt"); + assert(fdec != NULL); + } + + for(i=0; i<320/5; i++) + fprintf(fdec,"%f\t",Fw[i].real); + fprintf(fdec,"\n"); +} + +void dump_bg(float e, float bg_est, float percent_uv) { + char s[MAX_STR]; + + if (!dumpon) return; + + if (fbg == NULL) { + sprintf(s,"%s_bg.txt", prefix); + fbg = fopen(s, "wt"); + assert(fbg != NULL); + } + + fprintf(fbg,"%f\t%f\t%f\n", e, bg_est, percent_uv); +} + +void dump_E(float E) { + char s[MAX_STR]; + + if (!dumpon) return; + + if (fE == NULL) { + sprintf(s,"%s_E.txt", prefix); + fE = fopen(s, "wt"); + assert(fE != NULL); + } + + fprintf(fE,"%f\n", 10.0*log10(E)); +} + +void dump_Rk(float Rk[]) { + int i; + char s[MAX_STR]; + + if (!dumpon) return; + + if (frk == NULL) { + sprintf(s,"%s_rk.txt", prefix); + frk = fopen(s, "wt"); + assert(frk != NULL); + } + + for(i=0; i<P_MAX; i++) + fprintf(frk,"%f\t",Rk[i]); + fprintf(frk,"\n"); +} + +#endif diff --git a/gr-codec2-vocoder/src/lib/codec2/dump.h b/gr-codec2-vocoder/src/lib/codec2/dump.h new file mode 100644 index 000000000..eeddd3406 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/dump.h @@ -0,0 +1,67 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: dump.h + AUTHOR......: David Rowe + DATE CREATED: 25/8/09 + + Routines to dump data to text files for Octave analysis. + +\*---------------------------------------------------------------------------*/ + +/* + All rights reserved. + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +#ifndef __DUMP__ +#define __DUMP__ + +#include "comp.h" + +void dump_on(char filename_prefix[]); +void dump_off(); + +void dump_Sn(float Sn[]); +void dump_Sw(COMP Sw[]); +void dump_Sw_(COMP Sw_[]); +void dump_Ew(COMP Ew[]); + +/* amplitude modelling */ + +void dump_model(MODEL *m); +void dump_quantised_model(MODEL *m); +void dump_Pw(COMP Pw[]); +void dump_lsp(float lsp[]); +void dump_ak(float ak[], int order); +void dump_E(float E); +void dump_resample(float w[], float A[], int n); + +/* phase modelling */ + +void dump_snr(float snr); +void dump_phase(float phase[], int L); +void dump_phase_(float phase[], int L); + +/* NLP states */ + +void dump_sq(float sq[]); +void dump_dec(COMP Fw[]); +void dump_Fw(COMP Fw[]); +void dump_e(float e_hz[]); +void dump_Rk(float Rk[]); + +/* post filter */ + +void dump_bg(float e, float bg_est, float percent_uv); + +#endif diff --git a/gr-codec2-vocoder/src/lib/codec2/fft.c b/gr-codec2-vocoder/src/lib/codec2/fft.c new file mode 100644 index 000000000..73c46c846 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/fft.c @@ -0,0 +1,101 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: fft.c + AUTHOR......: Bruce Robertson + DATE CREATED: 20/11/2010 + + Bridging function to the kiss_fft package. + +\*---------------------------------------------------------------------------*/ + +/* + Copyright (C) 2010 Bruce Robertson + + All rights reserved. + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +#include <assert.h> +#include "kiss_fft.h" + +/*---------------------------------------------------------------------------*\ + + GLOBALS + +\*---------------------------------------------------------------------------*/ + +kiss_fft_cpx *fin; +kiss_fft_cpx *fout; +kiss_fft_cfg cfg_forward; +kiss_fft_cfg cfg_reverse; + +/*---------------------------------------------------------------------------*\ + + initialize_fft(int n) + + Initialisation function for kiss_fft. This assumes that all calls to fft() + use the same datatypes and are one arrays of the same size. + +\*---------------------------------------------------------------------------*/ + +void +initialize_fft (int n) +{ + fin = KISS_FFT_MALLOC (n * sizeof (kiss_fft_cpx)); + assert(fin != NULL); + fout = KISS_FFT_MALLOC (n * sizeof (kiss_fft_cpx)); + assert(fout != NULL); + cfg_forward = kiss_fft_alloc (n, 0, NULL, NULL); + assert(cfg_forward != NULL); + cfg_reverse = kiss_fft_alloc (n, 1, NULL, NULL); + assert(cfg_reverse != NULL); +} + +/*---------------------------------------------------------------------------*\ + + fft(float x[], int n, int isign) + Function that calls kiss_fft with the signature of four1 from NRC. + +\*---------------------------------------------------------------------------*/ + + +void +fft (float x[], int n, int isign) +{ + if (cfg_forward == NULL) + { + initialize_fft (n); + } + int isReverse = 0; + int c; + for (c = 0; c < n * 2; c += 2) + { + fin[c / 2].r = x[c]; + fin[c / 2].i = -x[c + 1]; + } + kiss_fft_cfg cfg; + if (isign == -1) + { + cfg = cfg_reverse; + } + else + { + cfg = cfg_forward; + } + kiss_fft (cfg, fin, fout); + for (c = 0; c < n * 2; c += 2) + { + x[c] = fout[(c) / 2].r; + x[c + 1] = -fout[(c) / 2].i; + } +} diff --git a/gr-codec2-vocoder/src/lib/codec2/fft.h b/gr-codec2-vocoder/src/lib/codec2/fft.h new file mode 100644 index 000000000..84c6737bd --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/fft.h @@ -0,0 +1,16 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: fft.h + AUTHOR......: Bruce Robertson + DATE CREATED: 29/11/2010 + + Bridge between existing code and kiss_fft. + +\*---------------------------------------------------------------------------*/ + +#ifndef __FFT__ +#define __FFT__ +void fft(float x[], int n, int isign); + +#endif /* __FFT__ */ + diff --git a/gr-codec2-vocoder/src/lib/codec2/fq20.sh b/gr-codec2-vocoder/src/lib/codec2/fq20.sh new file mode 100755 index 000000000..b83784b43 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/fq20.sh @@ -0,0 +1,8 @@ +#!/bin/sh +# fq20.shsh +# David Rowe 27 July 2010 +# +# Decode a file with fully quantised codec at 20ms frame rate + +../src/sinedec ../raw/$1.raw $1.mdl -o $1_phase0_lsp_20_EWo2.raw --phase 0 --lpc 10 --lsp --postfilter --dec + diff --git a/gr-codec2-vocoder/src/lib/codec2/generate_codebook.c b/gr-codec2-vocoder/src/lib/codec2/generate_codebook.c new file mode 100644 index 000000000..0bea80d85 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/generate_codebook.c @@ -0,0 +1,179 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: generate_codebook.c + AUTHOR......: Bruce Perens + DATE CREATED: 29 Sep 2010 + + Generate header files containing LSP quantisers, runs at compile time. + +\*---------------------------------------------------------------------------*/ + +/* + All rights reserved. + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +#include <stdlib.h> +#include <stdio.h> +#include <ctype.h> +#include <math.h> + +static const char usage[] = +"Usage: %s filename array_name [filename ...]\n" +"\tCreate C code for codebook tables.\n"; + +static const char format[] = +"The table format must be:\n" +"\tTwo integers describing the dimensions of the codebook.\n" +"\tThen, enough numbers to fill the specified dimensions.\n"; + +static const char header[] = +"/* THIS IS A GENERATED FILE. Edit generate_codebook.c and its input */\n\n" +"/*\n" +" * This intermediary file and the files that used to create it are under \n" +" * The LGPL. See the file COPYING.\n" +" */\n\n" +"#include \"defines.h\"\n\n"; + +struct codebook { + unsigned int k; + unsigned int log2m; + unsigned int m; + float * cb; +}; + +static void +dump_array(const struct codebook * b, int index) +{ + int limit = b->k * b->m; + int i; + + printf("static const float codes%d[] = {\n", index); + for ( i = 0; i < limit; i++ ) { + printf(" %g", b->cb[i]); + if ( i < limit - 1 ) + printf(","); + + /* organise VQs by rows, looks prettier */ + if ( ((i+1) % b->k) == 0 ) + printf("\n"); + } + printf("};\n"); +} + +static void +dump_structure(const struct codebook * b, int index) +{ + printf(" {\n"); + printf(" %d,\n", b->k); + printf(" %g,\n", log(b->m) / log(2)); + printf(" %d,\n", b->m); + printf(" codes%d\n", index); + printf(" }"); +} + +float +get_float(FILE * in, const char * name, char * * cursor, char * buffer, + int size) +{ + for ( ; ; ) { + char * s = *cursor; + char c; + + while ( (c = *s) != '\0' && !isdigit(c) && c != '-' && c != '.' ) + s++; + + /* Comments start with "#" and continue to the end of the line. */ + if ( c != '\0' && c != '#' ) { + char * end = 0; + float f = 0; + + f = strtod(s, &end); + + if ( end != s ) + *cursor = end; + return f; + } + + if ( fgets(buffer, size, in) == NULL ) { + fprintf(stderr, "%s: Format error. %s\n", name, format); + exit(1); + } + *cursor = buffer; + } +} + +static struct codebook * +load(FILE * file, const char * name) +{ + char line[1024]; + char * cursor = line; + struct codebook * b = malloc(sizeof(struct codebook)); + int i; + int size; + + *cursor = '\0'; + + b->k = (int)get_float(file, name, &cursor, line, sizeof(line)); + b->m = (int)get_float(file, name ,&cursor, line, sizeof(line)); + size = b->k * b->m; + + b->cb = (float *)malloc(size * sizeof(float)); + + for ( i = 0; i < size; i++ ) + b->cb[i] = get_float(file, name, &cursor, line, sizeof(line)); + + return b; +} + +int +main(int argc, char * * argv) +{ + struct codebook * * cb = malloc(argc * sizeof(struct codebook *)); + int i; + + if ( argc < 2 ) { + fprintf(stderr, usage, argv[0]); + fprintf(stderr, format); + exit(1); + } + + for ( i = 0; i < argc - 2; i++ ) { + FILE * in = fopen(argv[i + 2], "r"); + + if ( in == NULL ) { + perror(argv[i + 2]); + exit(1); + } + + cb[i] = load(in, argv[i + 2]); + + fclose(in); + } + + printf(header); + for ( i = 0; i < argc - 2; i++ ) { + printf(" /* %s */\n", argv[i + 2]); + dump_array(cb[i], i); + } + printf("\nconst struct lsp_codebook %s[] = {\n", argv[1]); + for ( i = 0; i < argc - 2; i++ ) { + printf(" /* %s */\n", argv[i + 2]); + dump_structure(cb[i], i); + printf(",\n"); + } + printf(" { 0, 0, 0, 0 }\n"); + printf("};\n"); + + return 0; +} diff --git a/gr-codec2-vocoder/src/lib/codec2/globals.c b/gr-codec2-vocoder/src/lib/codec2/globals.c new file mode 100644 index 000000000..f2182f79a --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/globals.c @@ -0,0 +1,49 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: globals.c + AUTHOR......: David Rowe + DATE CREATED: 11/5/94 + + Globals for sinusoidal speech coder. + +\*---------------------------------------------------------------------------*/ + +/* + 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 Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +#include "sine.h" /* global defines for coder */ + +/* Globals used in encoder and decoder */ + +int frames; /* number of frames processed so far */ +float Sn[M]; /* float input speech samples */ +MODEL model; /* model parameters for the current frame */ +int Nw; /* number of samples in analysis window */ +float sig; /* energy of current frame */ + +/* Globals used in encoder */ + +float w[M]; /* time domain hamming window */ +COMP W[FFT_ENC]; /* DFT of w[] */ +COMP Sw[FFT_ENC]; /* DFT of current frame */ + +/* Globals used in decoder */ + +COMP Sw_[FFT_ENC]; /* DFT of all voiced synthesised signal */ +float Sn_[AW_DEC]; /* synthesised speech */ +float Pn[AW_DEC]; /* time domain Parzen (trapezoidal) window */ + diff --git a/gr-codec2-vocoder/src/lib/codec2/globals.h b/gr-codec2-vocoder/src/lib/codec2/globals.h new file mode 100644 index 000000000..cef720344 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/globals.h @@ -0,0 +1,47 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: globals.h + AUTHOR......: David Rowe + DATE CREATED: 1/11/94 + + Globals for sinusoidal speech coder. + +\*---------------------------------------------------------------------------*/ + +/* + 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 Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +/* Globals used in encoder and decoder */ + +extern int frames; /* number of frames processed so far */ +extern float Sn[]; /* float input speech samples */ +extern MODEL model; /* model parameters for the current frame */ +extern int Nw; /* number of samples in analysis window */ +extern float sig; /* energy of current frame */ + +/* Globals used in encoder */ + +extern float w[]; /* time domain hamming window */ +extern COMP W[]; /* frequency domain hamming window */ +extern COMP Sw[]; /* DFT of current frame */ +extern COMP Sw_[]; /* DFT of all voiced synthesised signal */ + +/* Globals used in decoder */ + +extern float Sn_[]; /* output synthesised speech samples */ +extern float Pn[]; /* time domain Parzen (trapezoidal) window */ + diff --git a/gr-codec2-vocoder/src/lib/codec2/glottal.c b/gr-codec2-vocoder/src/lib/codec2/glottal.c new file mode 100644 index 000000000..8ac3ff4a9 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/glottal.c @@ -0,0 +1,257 @@ +const float glottal[]={ + 0.000000, + -0.057687, + -0.115338, + -0.172917, + -0.230385, + -0.287707, + -0.344845, + -0.401762, + -0.458419, + -0.514781, + -0.570809, + -0.626467, + -0.681721, + -0.736537, + -0.790884, + -0.844733, + -0.898057, + -0.950834, + -1.003044, + -1.054670, + -1.105700, + -1.156124, + -1.205936, + -1.255132, + -1.303711, + -1.351675, + -1.399026, + -1.445769, + -1.491908, + -1.537448, + -1.582393, + -1.626747, + -1.670514, + -1.713693, + -1.756285, + -1.798288, + -1.839697, + -1.880507, + -1.920712, + -1.960302, + -1.999269, + -2.037603, + -2.075295, + -2.112335, + -2.148716, + -2.184430, + -2.219472, + -2.253839, + -2.287531, + -2.320550, + -2.352900, + -2.384588, + -2.415623, + -2.446019, + -2.475788, + -2.504946, + -2.533512, + -2.561501, + -2.588934, + -2.615827, + -2.642198, + -2.668064, + -2.693439, + -2.718337, + -2.742767, + -2.766738, + -2.790256, + -2.813322, + -2.835936, + -2.858094, + -2.879790, + -2.901016, + -2.921759, + -2.942008, + -2.961747, + -2.980961, + -2.999632, + -3.017745, + -3.035282, + -3.052228, + -3.068567, + -3.084285, + -3.099371, + -3.113813, + -3.127605, + -3.140738, + 3.129975, + 3.118167, + 3.107022, + 3.096537, + 3.086709, + 3.077531, + 3.068996, + 3.061096, + 3.053821, + 3.047159, + 3.041102, + 3.035636, + 3.030753, + 3.026441, + 3.022690, + 3.019491, + 3.016836, + 3.014718, + 3.013132, + 3.012072, + 3.011535, + 3.011521, + 3.012028, + 3.013057, + 3.014612, + 3.016695, + 3.019310, + 3.022463, + 3.026160, + 3.030407, + 3.035212, + 3.040580, + 3.046520, + 3.053038, + 3.060141, + 3.067836, + 3.076128, + 3.085023, + 3.094525, + 3.104639, + 3.115367, + 3.126712, + 3.138674, + -3.131930, + -3.118731, + -3.104915, + -3.090485, + -3.075444, + -3.059795, + -3.043543, + -3.026695, + -3.009254, + -2.991229, + -2.972625, + -2.953449, + -2.933710, + -2.913414, + -2.892567, + -2.871176, + -2.849248, + -2.826787, + -2.803798, + -2.780284, + -2.756247, + -2.731689, + -2.706609, + -2.681005, + -2.654875, + -2.628213, + -2.601015, + -2.573272, + -2.544977, + -2.516121, + -2.486694, + -2.456686, + -2.426084, + -2.394879, + -2.363060, + -2.330616, + -2.297538, + -2.263816, + -2.229444, + -2.194416, + -2.158727, + -2.122375, + -2.085359, + -2.047682, + -2.009347, + -1.970361, + -1.930732, + -1.890470, + -1.849587, + -1.808098, + -1.766017, + -1.723360, + -1.680145, + -1.636388, + -1.592105, + -1.547313, + -1.502025, + -1.456256, + -1.410016, + -1.363314, + -1.316157, + -1.268547, + -1.220486, + -1.171971, + -1.122997, + -1.073555, + -1.023636, + -0.973227, + -0.922312, + -0.870875, + -0.818899, + -0.766366, + -0.713257, + -0.659554, + -0.605242, + -0.550303, + -0.494723, + -0.438492, + -0.381598, + -0.324036, + -0.265800, + -0.206889, + -0.147303, + -0.087046, + -0.026121, + 0.035463, + 0.097698, + 0.160576, + 0.224087, + 0.288221, + 0.352969, + 0.418323, + 0.484276, + 0.550822, + 0.617958, + 0.685681, + 0.753991, + 0.822889, + 0.892378, + 0.962462, + 1.033144, + 1.104430, + 1.176325, + 1.248833, + 1.321956, + 1.395696, + 1.470051, + 1.545019, + 1.620593, + 1.696763, + 1.773516, + 1.850837, + 1.928705, + 2.007097, + 2.085987, + 2.165347, + 2.245145, + 2.325347, + 2.405919, + 2.486824, + 2.568025, + 2.649485, + 2.731167, + 2.813033, + 2.895045, + 2.977167, + 3.059362}; diff --git a/gr-codec2-vocoder/src/lib/codec2/interp.c b/gr-codec2-vocoder/src/lib/codec2/interp.c new file mode 100644 index 000000000..135d8c9e7 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/interp.c @@ -0,0 +1,472 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: interp.c + AUTHOR......: David Rowe + DATE CREATED: 9/10/09 + + Interpolation of 20ms frames to 10ms frames. + +\*---------------------------------------------------------------------------*/ + +/* + 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 Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +#include <assert.h> +#include <math.h> +#include <string.h> +#include <stdio.h> +#include <stdlib.h> + +#include "defines.h" +#include "interp.h" +#include "lsp.h" +#include "quantise.h" +#include "dump.h" + +float sample_log_amp(MODEL *model, float w); + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: interp() + AUTHOR......: David Rowe + DATE CREATED: 22/8/10 + + Given two frames decribed by model parameters 20ms apart, determines + the model parameters of the 10ms frame between them. Assumes + voicing is available for middle (interpolated) frame. Outputs are + amplitudes and Wo for the interpolated frame. + + This version can interpolate the amplitudes between two frames of + different Wo and L. + + This version works by log linear interpolation, but listening tests + showed it creates problems in background noise, e.g. hts2a and mmt1. + When this function is used (--dec mode) bg noise appears to be + amplitude modulated, and gets louder. The interp_lsp() function + below seems to do a better job. + +\*---------------------------------------------------------------------------*/ + +void interpolate( + MODEL *interp, /* interpolated model params */ + MODEL *prev, /* previous frames model params */ + MODEL *next /* next frames model params */ +) +{ + int l; + float w,log_amp; + + /* Wo depends on voicing of this and adjacent frames */ + + if (interp->voiced) { + if (prev->voiced && next->voiced) + interp->Wo = (prev->Wo + next->Wo)/2.0; + if (!prev->voiced && next->voiced) + interp->Wo = next->Wo; + if (prev->voiced && !next->voiced) + interp->Wo = prev->Wo; + } + else { + interp->Wo = TWO_PI/P_MAX; + } + interp->L = PI/interp->Wo; + + /* Interpolate amplitudes using linear interpolation in log domain */ + + for(l=1; l<=interp->L; l++) { + w = l*interp->Wo; + log_amp = (sample_log_amp(prev, w) + sample_log_amp(next, w))/2.0; + interp->A[l] = pow(10.0, log_amp); + } +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: sample_log_amp() + AUTHOR......: David Rowe + DATE CREATED: 22/8/10 + + Samples the amplitude envelope at an arbitrary frequency w. Uses + linear interpolation in the log domain to sample between harmonic + amplitudes. + +\*---------------------------------------------------------------------------*/ + +float sample_log_amp(MODEL *model, float w) +{ + int m; + float f, log_amp; + + assert(w > 0.0); assert (w <= PI); + + m = 0; + while ((m+1)*model->Wo < w) m++; + f = (w - m*model->Wo)/model->Wo; + assert(f <= 1.0); + + if (m < 1) { + log_amp = f*log10(model->A[1] + 1E-6); + } + else if ((m+1) > model->L) { + log_amp = (1.0-f)*log10(model->A[model->L] + 1E-6); + } + else { + log_amp = (1.0-f)*log10(model->A[m] + 1E-6) + + f*log10(model->A[m+1] + 1E-6); + //printf("m=%d A[m] %f A[m+1] %f x %f %f %f\n", m, model->A[m], + // model->A[m+1], pow(10.0, log_amp), + // (1-f), f); + } + + return log_amp; +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: sample_log_amp_quad() + AUTHOR......: David Rowe + DATE CREATED: 9 March 2011 + + Samples the amplitude envelope at an arbitrary frequency w. Uses + quadratic interpolation in the log domain to sample between harmonic + amplitudes. + + y(x) = ax*x + bx + c + + We assume three points are x=-1, x=0, x=1, which we map to m-1,m,m+1 + + c = y(0) + b = (y(1) - y(-1))/2 + a = y(-1) + b - y(0) + +\*---------------------------------------------------------------------------*/ + +float sample_log_amp_quad(MODEL *model, float w) +{ + int m; + float a,b,c,x, log_amp; + + assert(w > 0.0); assert (w <= PI); + + m = floor(w/model->Wo + 0.5); + if (m < 2) m = 2; + if (m > (model->L-1)) m = model->L-1; + c = log10(model->A[m]+1E-6); + b = (log10(model->A[m+1]+1E-6) - log10(model->A[m-1]+1E-6))/2.0; + a = log10(model->A[m-1]+1E-6) + b - c; + x = (w - m*model->Wo)/model->Wo; + + log_amp = a*x*x + b*x + c; + //printf("m=%d A[m-1] %f A[m] %f A[m+1] %f w %f x %f log_amp %f\n", m, + // model->A[m-1], + // model->A[m], model->A[m+1], w, x, pow(10.0, log_amp)); + return log_amp; +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: sample_log_amp_quad_nl() + AUTHOR......: David Rowe + DATE CREATED: 10 March 2011 + + Samples the amplitude envelope at an arbitrary frequency w. Uses + quadratic interpolation in the log domain to sample between harmonic + amplitudes. This version can handle non-linear steps along a freq + axis defined by arbitrary steps. + + y(x) = ax*x + bx + c + + We assume three points are (x_1,y_1), (0,y0) and (x1,y1). + +\*---------------------------------------------------------------------------*/ + +float sample_log_amp_quad_nl( + float w[], /* frequency points */ + float A[], /* for these amplitude samples */ + int np, /* number of frequency points */ + float w_sample /* frequency of new samples */ +) +{ + int m,i; + float a,b,c,x, log_amp, best_dist; + float x_1, x1; + float y_1, y0, y1; + + //printf("w_sample %f\n", w_sample); + assert(w_sample >= 0.0); assert (w_sample <= 1.1*PI); + + /* find closest point to centre quadratic interpolator */ + + best_dist = 1E32; + for (i=0; i<np; i++) + if (fabs(w[i] - w_sample) < best_dist) { + best_dist = fabs(w[i] - w_sample); + m = i; + } + + /* stay one point away from edge of array */ + + if (m < 1) m = 1; + if (m > (np-2)) m = np - 2; + + /* find polynomial coeffs */ + + x_1 = w[m-1]- w[m]; x1 = w[m+1] - w[m]; + y_1 = log10(A[m-1]+1E-6); + y0 = log10(A[m]+1E-6); + y1 = log10(A[m+1]+1E-6); + + c = y0; + a = (y_1*x1 - y1*x_1 + c*x_1 - c*x1)/(x_1*x_1*x1 - x1*x1*x_1); + b = (y1 -a*x1*x1 - c)/x1; + x = w_sample - w[m]; + + //printf("%f %f %f\n", w[0], w[1], w[2]); + //printf("%f %f %f %f %f %f\n", x_1, y_1, 0.0, y0, x1, y1); + log_amp = a*x*x + b*x + c; + //printf("a %f b %f c %f\n", a, b, c); + //printf("m=%d A[m-1] %f A[m] %f A[m+1] %f w_sample %f w[m] %f x %f log_amp %f\n", m, + // A[m-1], + // A[m], A[m+1], w_sample, w[m], x, log_amp); + //exit(0); + return log_amp; +} + +#define M_MAX 40 + +float fres[] = {100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, + 1200, 1400, 1600, 1850, 2100, 2350, 2600, 2900, 3400, 3800}; + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: resample_amp_nl() + AUTHOR......: David Rowe + DATE CREATED: 7 March 2011 + + Converts the current model with L {Am} samples spaced Wo apart to + RES_POINTS samples spaced Wo/RES_POINTS apart. Then subtracts + from the previous frames samples to get the delta. + +\*---------------------------------------------------------------------------*/ + +void resample_amp_fixed(MODEL *model, + float w[], float A[], + float wres[], float Ares[], + float AresdB_prev[], + float AresdB[], + float deltat[]) +{ + int i; + + for(i=1; i<=model->L; i++) { + w[i-1] = i*model->Wo; + A[i-1] = model->A[i]; + } + + for(i=0; i<RES_POINTS; i++) { + wres[i] = fres[i]*PI/4000.0; + } + + for(i=0; i<RES_POINTS; i++) { + Ares[i] = pow(10.0,sample_log_amp_quad_nl(w, A, model->L, wres[i])); + } + + /* work out delta T vector for this frame */ + + for(i=0; i<RES_POINTS; i++) { + AresdB[i] = 20.0*log10(Ares[i]); + deltat[i] = AresdB[i] - AresdB_prev[i]; + } + +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: resample_amp_nl() + AUTHOR......: David Rowe + DATE CREATED: 7 March 2011 + + Converts the current model with L {Am} samples spaced Wo apart to M + samples spaced Wo/M apart. Then converts back to L {Am} samples. + used to prototype constant rate Amplitude encoding ideas. + + Returns the SNR in dB. + +\*---------------------------------------------------------------------------*/ + +float resample_amp_nl(MODEL *model, int m, float AresdB_prev[]) +{ + int i; + float w[MAX_AMP], A[MAX_AMP]; + float wres[MAX_AMP], Ares[MAX_AMP], AresdB[MAX_AMP]; + float signal, noise, snr; + float new_A; + float deltat[MAX_AMP], deltat_q[MAX_AMP], AresdB_q[MAX_AMP]; + + resample_amp_fixed(model, w, A, wres, Ares, AresdB_prev, AresdB, deltat); + + /* quantise delta T vector */ + + for(i=0; i<RES_POINTS; i++) { + noise = 3.0*(1.0 - 2.0*rand()/RAND_MAX); + //noise = 0.0; + deltat_q[i] = deltat[i] + noise; + } + + /* recover Ares vector */ + + for(i=0; i<RES_POINTS; i++) { + AresdB_q[i] = AresdB_prev[i] + deltat_q[i]; + Ares[i] = pow(10.0, AresdB_q[i]/20.0); + //printf("%d %f %f\n", i, AresdB[i], AresdB_q[i]); + } + + /* update memory based on version at decoder */ + + for(i=0; i<RES_POINTS; i++) { + AresdB_prev[i] = AresdB_q[i]; + } + +#ifdef DUMP + dump_resample(wres,Ares,M_MAX); +#endif + + signal = noise = 0.0; + + for(i=1; i<model->L; i++) { + new_A = pow(10.0,sample_log_amp_quad_nl(wres, Ares, RES_POINTS, model->Wo*i)); + signal += pow(model->A[i], 2.0); + noise += pow(model->A[i] - new_A, 2.0); + //printf("%f %f\n", model->A[i], new_A); + model->A[i] = new_A; + } + + snr = 10.0*log10(signal/noise); + printf("snr = %3.2f\n", snr); + //exit(0); + return snr; +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: resample_amp() + AUTHOR......: David Rowe + DATE CREATED: 10 March 2011 + + Converts the current model with L {Am} samples spaced Wo apart to M + samples with a non-linear spacing. Then converts back to L {Am} + samples. used to prototype constant rate Amplitude encoding ideas. + + Returns the SNR in dB. + +\*---------------------------------------------------------------------------*/ + +float resample_amp(MODEL *model, int m) +{ + int i; + MODEL model_m; + float new_A, signal, noise, snr, log_amp_dB; + float n_db = 0.0; + + model_m.Wo = PI/(float)m; + model_m.L = PI/model_m.Wo; + + for(i=1; i<=model_m.L; i++) { + log_amp_dB = 20.0*sample_log_amp_quad(model, i*model_m.Wo); + log_amp_dB += n_db*(1.0 - 2.0*rand()/RAND_MAX); + model_m.A[i] = pow(10,log_amp_dB/20.0); + } + + //dump_resample(&model_m); + + signal = noise = 0.0; + + for(i=1; i<model->L/4; i++) { + new_A = pow(10,sample_log_amp_quad(&model_m, i*model->Wo)); + signal += pow(model->A[i], 2.0); + noise += pow(model->A[i] - new_A, 2.0); + //printf("%f %f\n", model->A[i], new_A); + model->A[i] = new_A; + } + + snr = 10.0*log10(signal/noise); + //printf("snr = %3.2f\n", snr); + //exit(0); + return snr; +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: interp_lsp() + AUTHOR......: David Rowe + DATE CREATED: 10 Nov 2010 + + Given two frames decribed by model parameters 20ms apart, determines + the model parameters of the 10ms frame between them. Assumes + voicing is available for middle (interpolated) frame. Outputs are + amplitudes and Wo for the interpolated frame. + + This version uses interpolation of LSPs, seems to do a better job + with bg noise. + +\*---------------------------------------------------------------------------*/ + +void interpolate_lsp( + MODEL *interp, /* interpolated model params */ + MODEL *prev, /* previous frames model params */ + MODEL *next, /* next frames model params */ + float *prev_lsps, /* previous frames LSPs */ + float prev_e, /* previous frames LPC energy */ + float *next_lsps, /* next frames LSPs */ + float next_e, /* next frames LPC energy */ + float *ak_interp /* interpolated aks for this frame */ + ) +{ + int l,i; + float lsps[LPC_ORD],e; + float snr; + + /* Wo depends on voicing of this and adjacent frames */ + + if (interp->voiced) { + if (prev->voiced && next->voiced) + interp->Wo = (prev->Wo + next->Wo)/2.0; + if (!prev->voiced && next->voiced) + interp->Wo = next->Wo; + if (prev->voiced && !next->voiced) + interp->Wo = prev->Wo; + } + else { + interp->Wo = TWO_PI/P_MAX; + } + interp->L = PI/interp->Wo; + + /* interpolate LSPs */ + + for(i=0; i<LPC_ORD; i++) { + lsps[i] = (prev_lsps[i] + next_lsps[i])/2.0; + } + + /* Interpolate LPC energy in log domain */ + + e = pow(10.0, (log10(prev_e) + log10(next_e))/2.0); + + /* convert back to amplitudes */ + + lsp_to_lpc(lsps, ak_interp, LPC_ORD); + aks_to_M2(ak_interp, LPC_ORD, interp, e, &snr, 0); +} diff --git a/gr-codec2-vocoder/src/lib/codec2/interp.h b/gr-codec2-vocoder/src/lib/codec2/interp.h new file mode 100644 index 000000000..d41eac3f8 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/interp.h @@ -0,0 +1,41 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: interp.h + AUTHOR......: David Rowe + DATE CREATED: 9/10/09 + + Interpolation of 20ms frames to 10ms frames. + +\*---------------------------------------------------------------------------*/ + +/* + 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 Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +#ifndef __INTERP__ +#define __INTERP__ + +#define RES_POINTS 20 + +void interpolate(MODEL *interp, MODEL *prev, MODEL *next); +void interpolate_lsp(MODEL *interp, MODEL *prev, MODEL *next, + float *prev_lsps, float prev_e, + float *next_lsps, float next_e, + float *ak_interp); +float resample_amp(MODEL *model, int m); +float resample_amp_nl(MODEL *model, int m, float Ares_prev[]); + +#endif diff --git a/gr-codec2-vocoder/src/lib/codec2/kiss_fft.c b/gr-codec2-vocoder/src/lib/codec2/kiss_fft.c new file mode 100644 index 000000000..465d6c97a --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/kiss_fft.c @@ -0,0 +1,408 @@ +/* +Copyright (c) 2003-2010, Mark Borgerding + +All rights reserved. + +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 author nor the names of any 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 COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +*/ + + +#include "_kiss_fft_guts.h" +/* The guts header contains all the multiplication and addition macros that are defined for + fixed or floating point complex numbers. It also delares the kf_ internal functions. + */ + +static void kf_bfly2( + kiss_fft_cpx * Fout, + const size_t fstride, + const kiss_fft_cfg st, + int m + ) +{ + kiss_fft_cpx * Fout2; + kiss_fft_cpx * tw1 = st->twiddles; + kiss_fft_cpx t; + Fout2 = Fout + m; + do{ + C_FIXDIV(*Fout,2); C_FIXDIV(*Fout2,2); + + C_MUL (t, *Fout2 , *tw1); + tw1 += fstride; + C_SUB( *Fout2 , *Fout , t ); + C_ADDTO( *Fout , t ); + ++Fout2; + ++Fout; + }while (--m); +} + +static void kf_bfly4( + kiss_fft_cpx * Fout, + const size_t fstride, + const kiss_fft_cfg st, + const size_t m + ) +{ + kiss_fft_cpx *tw1,*tw2,*tw3; + kiss_fft_cpx scratch[6]; + size_t k=m; + const size_t m2=2*m; + const size_t m3=3*m; + + + tw3 = tw2 = tw1 = st->twiddles; + + do { + C_FIXDIV(*Fout,4); C_FIXDIV(Fout[m],4); C_FIXDIV(Fout[m2],4); C_FIXDIV(Fout[m3],4); + + C_MUL(scratch[0],Fout[m] , *tw1 ); + C_MUL(scratch[1],Fout[m2] , *tw2 ); + C_MUL(scratch[2],Fout[m3] , *tw3 ); + + C_SUB( scratch[5] , *Fout, scratch[1] ); + C_ADDTO(*Fout, scratch[1]); + C_ADD( scratch[3] , scratch[0] , scratch[2] ); + C_SUB( scratch[4] , scratch[0] , scratch[2] ); + C_SUB( Fout[m2], *Fout, scratch[3] ); + tw1 += fstride; + tw2 += fstride*2; + tw3 += fstride*3; + C_ADDTO( *Fout , scratch[3] ); + + if(st->inverse) { + Fout[m].r = scratch[5].r - scratch[4].i; + Fout[m].i = scratch[5].i + scratch[4].r; + Fout[m3].r = scratch[5].r + scratch[4].i; + Fout[m3].i = scratch[5].i - scratch[4].r; + }else{ + Fout[m].r = scratch[5].r + scratch[4].i; + Fout[m].i = scratch[5].i - scratch[4].r; + Fout[m3].r = scratch[5].r - scratch[4].i; + Fout[m3].i = scratch[5].i + scratch[4].r; + } + ++Fout; + }while(--k); +} + +static void kf_bfly3( + kiss_fft_cpx * Fout, + const size_t fstride, + const kiss_fft_cfg st, + size_t m + ) +{ + size_t k=m; + const size_t m2 = 2*m; + kiss_fft_cpx *tw1,*tw2; + kiss_fft_cpx scratch[5]; + kiss_fft_cpx epi3; + epi3 = st->twiddles[fstride*m]; + + tw1=tw2=st->twiddles; + + do{ + C_FIXDIV(*Fout,3); C_FIXDIV(Fout[m],3); C_FIXDIV(Fout[m2],3); + + C_MUL(scratch[1],Fout[m] , *tw1); + C_MUL(scratch[2],Fout[m2] , *tw2); + + C_ADD(scratch[3],scratch[1],scratch[2]); + C_SUB(scratch[0],scratch[1],scratch[2]); + tw1 += fstride; + tw2 += fstride*2; + + Fout[m].r = Fout->r - HALF_OF(scratch[3].r); + Fout[m].i = Fout->i - HALF_OF(scratch[3].i); + + C_MULBYSCALAR( scratch[0] , epi3.i ); + + C_ADDTO(*Fout,scratch[3]); + + Fout[m2].r = Fout[m].r + scratch[0].i; + Fout[m2].i = Fout[m].i - scratch[0].r; + + Fout[m].r -= scratch[0].i; + Fout[m].i += scratch[0].r; + + ++Fout; + }while(--k); +} + +static void kf_bfly5( + kiss_fft_cpx * Fout, + const size_t fstride, + const kiss_fft_cfg st, + int m + ) +{ + kiss_fft_cpx *Fout0,*Fout1,*Fout2,*Fout3,*Fout4; + int u; + kiss_fft_cpx scratch[13]; + kiss_fft_cpx * twiddles = st->twiddles; + kiss_fft_cpx *tw; + kiss_fft_cpx ya,yb; + ya = twiddles[fstride*m]; + yb = twiddles[fstride*2*m]; + + Fout0=Fout; + Fout1=Fout0+m; + Fout2=Fout0+2*m; + Fout3=Fout0+3*m; + Fout4=Fout0+4*m; + + tw=st->twiddles; + for ( u=0; u<m; ++u ) { + C_FIXDIV( *Fout0,5); C_FIXDIV( *Fout1,5); C_FIXDIV( *Fout2,5); C_FIXDIV( *Fout3,5); C_FIXDIV( *Fout4,5); + scratch[0] = *Fout0; + + C_MUL(scratch[1] ,*Fout1, tw[u*fstride]); + C_MUL(scratch[2] ,*Fout2, tw[2*u*fstride]); + C_MUL(scratch[3] ,*Fout3, tw[3*u*fstride]); + C_MUL(scratch[4] ,*Fout4, tw[4*u*fstride]); + + C_ADD( scratch[7],scratch[1],scratch[4]); + C_SUB( scratch[10],scratch[1],scratch[4]); + C_ADD( scratch[8],scratch[2],scratch[3]); + C_SUB( scratch[9],scratch[2],scratch[3]); + + Fout0->r += scratch[7].r + scratch[8].r; + Fout0->i += scratch[7].i + scratch[8].i; + + scratch[5].r = scratch[0].r + S_MUL(scratch[7].r,ya.r) + S_MUL(scratch[8].r,yb.r); + scratch[5].i = scratch[0].i + S_MUL(scratch[7].i,ya.r) + S_MUL(scratch[8].i,yb.r); + + scratch[6].r = S_MUL(scratch[10].i,ya.i) + S_MUL(scratch[9].i,yb.i); + scratch[6].i = -S_MUL(scratch[10].r,ya.i) - S_MUL(scratch[9].r,yb.i); + + C_SUB(*Fout1,scratch[5],scratch[6]); + C_ADD(*Fout4,scratch[5],scratch[6]); + + scratch[11].r = scratch[0].r + S_MUL(scratch[7].r,yb.r) + S_MUL(scratch[8].r,ya.r); + scratch[11].i = scratch[0].i + S_MUL(scratch[7].i,yb.r) + S_MUL(scratch[8].i,ya.r); + scratch[12].r = - S_MUL(scratch[10].i,yb.i) + S_MUL(scratch[9].i,ya.i); + scratch[12].i = S_MUL(scratch[10].r,yb.i) - S_MUL(scratch[9].r,ya.i); + + C_ADD(*Fout2,scratch[11],scratch[12]); + C_SUB(*Fout3,scratch[11],scratch[12]); + + ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4; + } +} + +/* perform the butterfly for one stage of a mixed radix FFT */ +static void kf_bfly_generic( + kiss_fft_cpx * Fout, + const size_t fstride, + const kiss_fft_cfg st, + int m, + int p + ) +{ + int u,k,q1,q; + kiss_fft_cpx * twiddles = st->twiddles; + kiss_fft_cpx t; + int Norig = st->nfft; + + kiss_fft_cpx * scratch = (kiss_fft_cpx*)KISS_FFT_TMP_ALLOC(sizeof(kiss_fft_cpx)*p); + + for ( u=0; u<m; ++u ) { + k=u; + for ( q1=0 ; q1<p ; ++q1 ) { + scratch[q1] = Fout[ k ]; + C_FIXDIV(scratch[q1],p); + k += m; + } + + k=u; + for ( q1=0 ; q1<p ; ++q1 ) { + int twidx=0; + Fout[ k ] = scratch[0]; + for (q=1;q<p;++q ) { + twidx += fstride * k; + if (twidx>=Norig) twidx-=Norig; + C_MUL(t,scratch[q] , twiddles[twidx] ); + C_ADDTO( Fout[ k ] ,t); + } + k += m; + } + } + KISS_FFT_TMP_FREE(scratch); +} + +static +void kf_work( + kiss_fft_cpx * Fout, + const kiss_fft_cpx * f, + const size_t fstride, + int in_stride, + int * factors, + const kiss_fft_cfg st + ) +{ + kiss_fft_cpx * Fout_beg=Fout; + const int p=*factors++; /* the radix */ + const int m=*factors++; /* stage's fft length/p */ + const kiss_fft_cpx * Fout_end = Fout + p*m; + +#ifdef _OPENMP + // use openmp extensions at the + // top-level (not recursive) + if (fstride==1 && p<=5) + { + int k; + + // execute the p different work units in different threads +# pragma omp parallel for + for (k=0;k<p;++k) + kf_work( Fout +k*m, f+ fstride*in_stride*k,fstride*p,in_stride,factors,st); + // all threads have joined by this point + + switch (p) { + case 2: kf_bfly2(Fout,fstride,st,m); break; + case 3: kf_bfly3(Fout,fstride,st,m); break; + case 4: kf_bfly4(Fout,fstride,st,m); break; + case 5: kf_bfly5(Fout,fstride,st,m); break; + default: kf_bfly_generic(Fout,fstride,st,m,p); break; + } + return; + } +#endif + + if (m==1) { + do{ + *Fout = *f; + f += fstride*in_stride; + }while(++Fout != Fout_end ); + }else{ + do{ + // recursive call: + // DFT of size m*p performed by doing + // p instances of smaller DFTs of size m, + // each one takes a decimated version of the input + kf_work( Fout , f, fstride*p, in_stride, factors,st); + f += fstride*in_stride; + }while( (Fout += m) != Fout_end ); + } + + Fout=Fout_beg; + + // recombine the p smaller DFTs + switch (p) { + case 2: kf_bfly2(Fout,fstride,st,m); break; + case 3: kf_bfly3(Fout,fstride,st,m); break; + case 4: kf_bfly4(Fout,fstride,st,m); break; + case 5: kf_bfly5(Fout,fstride,st,m); break; + default: kf_bfly_generic(Fout,fstride,st,m,p); break; + } +} + +/* facbuf is populated by p1,m1,p2,m2, ... + where + p[i] * m[i] = m[i-1] + m0 = n */ +static +void kf_factor(int n,int * facbuf) +{ + int p=4; + double floor_sqrt; + floor_sqrt = floor( sqrt((double)n) ); + + /*factor out powers of 4, powers of 2, then any remaining primes */ + do { + while (n % p) { + switch (p) { + case 4: p = 2; break; + case 2: p = 3; break; + default: p += 2; break; + } + if (p > floor_sqrt) + p = n; /* no more factors, skip to end */ + } + n /= p; + *facbuf++ = p; + *facbuf++ = n; + } while (n > 1); +} + +/* + * + * User-callable function to allocate all necessary storage space for the fft. + * + * The return value is a contiguous block of memory, allocated with malloc. As such, + * It can be freed with free(), rather than a kiss_fft-specific function. + * */ +kiss_fft_cfg kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem ) +{ + kiss_fft_cfg st=NULL; + size_t memneeded = sizeof(struct kiss_fft_state) + + sizeof(kiss_fft_cpx)*(nfft-1); /* twiddle factors*/ + + if ( lenmem==NULL ) { + st = ( kiss_fft_cfg)KISS_FFT_MALLOC( memneeded ); + }else{ + if (mem != NULL && *lenmem >= memneeded) + st = (kiss_fft_cfg)mem; + *lenmem = memneeded; + } + if (st) { + int i; + st->nfft=nfft; + st->inverse = inverse_fft; + + for (i=0;i<nfft;++i) { + const double pi=3.141592653589793238462643383279502884197169399375105820974944; + double phase = -2*pi*i / nfft; + if (st->inverse) + phase *= -1; + kf_cexp(st->twiddles+i, phase ); + } + + kf_factor(nfft,st->factors); + } + return st; +} + + +void kiss_fft_stride(kiss_fft_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int in_stride) +{ + if (fin == fout) { + //NOTE: this is not really an in-place FFT algorithm. + //It just performs an out-of-place FFT into a temp buffer + kiss_fft_cpx * tmpbuf = (kiss_fft_cpx*)KISS_FFT_TMP_ALLOC( sizeof(kiss_fft_cpx)*st->nfft); + kf_work(tmpbuf,fin,1,in_stride, st->factors,st); + memcpy(fout,tmpbuf,sizeof(kiss_fft_cpx)*st->nfft); + KISS_FFT_TMP_FREE(tmpbuf); + }else{ + kf_work( fout, fin, 1,in_stride, st->factors,st ); + } +} + +void kiss_fft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout) +{ + kiss_fft_stride(cfg,fin,fout,1); +} + + +void kiss_fft_cleanup(void) +{ + // nothing needed any more +} + +int kiss_fft_next_fast_size(int n) +{ + while(1) { + int m=n; + while ( (m%2) == 0 ) m/=2; + while ( (m%3) == 0 ) m/=3; + while ( (m%5) == 0 ) m/=5; + if (m<=1) + break; /* n is completely factorable by twos, threes, and fives */ + n++; + } + return n; +} diff --git a/gr-codec2-vocoder/src/lib/codec2/kiss_fft.h b/gr-codec2-vocoder/src/lib/codec2/kiss_fft.h new file mode 100644 index 000000000..20621d8b3 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/kiss_fft.h @@ -0,0 +1,125 @@ +#ifndef KISS_FFT_H +#define KISS_FFT_H + +#include <stdlib.h> +#include <stdio.h> +#include <math.h> +#include <string.h> +#include <malloc.h> + +#ifdef __cplusplus +extern "C" { +#endif + +/* + ATTENTION! + If you would like a : + -- a utility that will handle the caching of fft objects + -- real-only (no imaginary time component ) FFT + -- a multi-dimensional FFT + -- a command-line utility to perform ffts + -- a command-line utility to perform fast-convolution filtering + + Then see kfc.h kiss_fftr.h kiss_fftnd.h fftutil.c kiss_fastfir.c + in the tools/ directory. +*/ + +#ifdef USE_SIMD +# include <xmmintrin.h> +# define kiss_fft_scalar __m128 +#define KISS_FFT_MALLOC(nbytes) _mm_malloc(nbytes,16) +#define KISS_FFT_FREE _mm_free +#else +#define KISS_FFT_MALLOC malloc +#define KISS_FFT_FREE free +#endif + + +#ifdef FIXED_POINT +#include <sys/types.h> +# if (FIXED_POINT == 32) +# define kiss_fft_scalar int32_t +# else +# define kiss_fft_scalar int16_t +# endif +#else +# ifndef kiss_fft_scalar +/* default is float */ +# define kiss_fft_scalar float +# endif +#endif + +typedef struct { + kiss_fft_scalar r; + kiss_fft_scalar i; +}kiss_fft_cpx; + +typedef struct kiss_fft_state* kiss_fft_cfg; + +/* + * kiss_fft_alloc + * + * Initialize a FFT (or IFFT) algorithm's cfg/state buffer. + * + * typical usage: kiss_fft_cfg mycfg=kiss_fft_alloc(1024,0,NULL,NULL); + * + * The return value from fft_alloc is a cfg buffer used internally + * by the fft routine or NULL. + * + * If lenmem is NULL, then kiss_fft_alloc will allocate a cfg buffer using malloc. + * The returned value should be free()d when done to avoid memory leaks. + * + * The state can be placed in a user supplied buffer 'mem': + * If lenmem is not NULL and mem is not NULL and *lenmem is large enough, + * then the function places the cfg in mem and the size used in *lenmem + * and returns mem. + * + * If lenmem is not NULL and ( mem is NULL or *lenmem is not large enough), + * then the function returns NULL and places the minimum cfg + * buffer size in *lenmem. + * */ + +kiss_fft_cfg kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem); + +/* + * kiss_fft(cfg,in_out_buf) + * + * Perform an FFT on a complex input buffer. + * for a forward FFT, + * fin should be f[0] , f[1] , ... ,f[nfft-1] + * fout will be F[0] , F[1] , ... ,F[nfft-1] + * Note that each element is complex and can be accessed like + f[k].r and f[k].i + * */ +void kiss_fft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout); + +/* + A more generic version of the above function. It reads its input from every Nth sample. + * */ +void kiss_fft_stride(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int fin_stride); + +/* If kiss_fft_alloc allocated a buffer, it is one contiguous + buffer and can be simply free()d when no longer needed*/ +#define kiss_fft_free free + +/* + Cleans up some memory that gets managed internally. Not necessary to call, but it might clean up + your compiler output to call this before you exit. +*/ +void kiss_fft_cleanup(void); + + +/* + * Returns the smallest integer k, such that k>=n and k has only "fast" factors (2,3,5) + */ +int kiss_fft_next_fast_size(int n); + +/* for real ffts, we need an even size */ +#define kiss_fftr_next_fast_size_real(n) \ + (kiss_fft_next_fast_size( ((n)+1)>>1)<<1) + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/gr-codec2-vocoder/src/lib/codec2/listensim.sh b/gr-codec2-vocoder/src/lib/codec2/listensim.sh new file mode 100755 index 000000000..0b27a1b0c --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/listensim.sh @@ -0,0 +1,9 @@ +#!/bin/sh +# listensim.sh +# David Rowe 10 Sep 2009 +# +# Listen to files processed with sim.sh + +../script/menu.sh ../raw/$1.raw $1_uq.raw $1_phase0.raw $1_lpc10.raw $1_phase0_lpc10.raw $1_phase0_lpc10_dec.raw $1_phase0_lsp_dec.raw $2 $3 + + diff --git a/gr-codec2-vocoder/src/lib/codec2/lpc.c b/gr-codec2-vocoder/src/lib/codec2/lpc.c new file mode 100644 index 000000000..ba8011377 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/lpc.c @@ -0,0 +1,279 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: lpc.c + AUTHOR......: David Rowe + DATE CREATED: 30/9/90 + + Linear Prediction functions written in C. + +\*---------------------------------------------------------------------------*/ + +/* + 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 Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +#define LPC_MAX_N 512 /* maximum no. of samples in frame */ +#define PI 3.141592654 /* mathematical constant */ + +#include <assert.h> +#include <math.h> +#include "defines.h" +#include "lpc.h" + +/*---------------------------------------------------------------------------*\ + + hanning_window() + + Hanning windows a frame of speech samples. + +\*---------------------------------------------------------------------------*/ + +void hanning_window( + float Sn[], /* input frame of speech samples */ + float Wn[], /* output frame of windowed samples */ + int Nsam /* number of samples */ +) +{ + int i; /* loop variable */ + + for(i=0; i<Nsam; i++) + Wn[i] = Sn[i]*(0.5 - 0.5*cos(2*PI*(float)i/(Nsam-1))); +} + +/*---------------------------------------------------------------------------*\ + + autocorrelate() + + Finds the first P autocorrelation values of an array of windowed speech + samples Sn[]. + +\*---------------------------------------------------------------------------*/ + +void autocorrelate( + float Sn[], /* frame of Nsam windowed speech samples */ + float Rn[], /* array of P+1 autocorrelation coefficients */ + int Nsam, /* number of windowed samples to use */ + int order /* order of LPC analysis */ +) +{ + int i,j; /* loop variables */ + + for(j=0; j<order+1; j++) { + Rn[j] = 0.0; + for(i=0; i<Nsam-j; i++) + Rn[j] += Sn[i]*Sn[i+j]; + } +} + +/*---------------------------------------------------------------------------*\ + + autocorrelate_freq() + + Finds the first P autocorrelation values from an array of frequency domain + power samples. + +\*---------------------------------------------------------------------------*/ + +void autocorrelate_freq( + float Pw[], /* Nsam frequency domain power spectrum samples */ + float w[], /* frequency of each sample in Pw[] */ + float R[], /* array of order+1 autocorrelation coefficients */ + int Nsam, /* number of windowed samples to use */ + int order /* order of LPC analysis */ +) +{ + int i,j; /* loop variables */ + + for(j=0; j<order+1; j++) { + R[j] = 0.0; + for(i=0; i<Nsam; i++) + R[j] += Pw[i]*cos(j*w[i]); + } + R[j] /= Nsam; +} + +/*---------------------------------------------------------------------------*\ + + levinson_durbin() + + Given P+1 autocorrelation coefficients, finds P Linear Prediction Coeff. + (LPCs) where P is the order of the LPC all-pole model. The Levinson-Durbin + algorithm is used, and is described in: + + J. Makhoul + "Linear prediction, a tutorial review" + Proceedings of the IEEE + Vol-63, No. 4, April 1975 + +\*---------------------------------------------------------------------------*/ + +void levinson_durbin( + float R[], /* order+1 autocorrelation coeff */ + float lpcs[], /* order+1 LPC's */ + int order /* order of the LPC analysis */ +) +{ + float E[LPC_MAX+1]; + float k[LPC_MAX+1]; + float a[LPC_MAX+1][LPC_MAX+1]; + float sum; + int i,j; /* loop variables */ + + E[0] = R[0]; /* Equation 38a, Makhoul */ + + for(i=1; i<=order; i++) { + sum = 0.0; + for(j=1; j<=i-1; j++) + sum += a[i-1][j]*R[i-j]; + k[i] = -1.0*(R[i] + sum)/E[i-1]; /* Equation 38b, Makhoul */ + if (fabs(k[i]) > 1.0) + k[i] = 0.0; + + a[i][i] = k[i]; + + for(j=1; j<=i-1; j++) + a[i][j] = a[i-1][j] + k[i]*a[i-1][i-j]; /* Equation 38c, Makhoul */ + + E[i] = (1-k[i]*k[i])*E[i-1]; /* Equation 38d, Makhoul */ + } + + for(i=1; i<=order; i++) + lpcs[i] = a[order][i]; + lpcs[0] = 1.0; +} + +/*---------------------------------------------------------------------------*\ + + inverse_filter() + + Inverse Filter, A(z). Produces an array of residual samples from an array + of input samples and linear prediction coefficients. + + The filter memory is stored in the first order samples of the input array. + +\*---------------------------------------------------------------------------*/ + +void inverse_filter( + float Sn[], /* Nsam input samples */ + float a[], /* LPCs for this frame of samples */ + int Nsam, /* number of samples */ + float res[], /* Nsam residual samples */ + int order /* order of LPC */ +) +{ + int i,j; /* loop variables */ + + for(i=0; i<Nsam; i++) { + res[i] = 0.0; + for(j=0; j<=order; j++) + res[i] += Sn[i-j]*a[j]; + } +} + +/*---------------------------------------------------------------------------*\ + + synthesis_filter() + + C version of the Speech Synthesis Filter, 1/A(z). Given an array of + residual or excitation samples, and the the LP filter coefficients, this + function will produce an array of speech samples. This filter structure is + IIR. + + The synthesis filter has memory as well, this is treated in the same way + as the memory for the inverse filter (see inverse_filter() notes above). + The difference is that the memory for the synthesis filter is stored in + the output array, wheras the memory of the inverse filter is stored in the + input array. + + Note: the calling function must update the filter memory. + +\*---------------------------------------------------------------------------*/ + +void synthesis_filter( + float res[], /* Nsam input residual (excitation) samples */ + float a[], /* LPCs for this frame of speech samples */ + int Nsam, /* number of speech samples */ + int order, /* LPC order */ + float Sn_[] /* Nsam output synthesised speech samples */ +) +{ + int i,j; /* loop variables */ + + /* Filter Nsam samples */ + + for(i=0; i<Nsam; i++) { + Sn_[i] = res[i]*a[0]; + for(j=1; j<=order; j++) + Sn_[i] -= Sn_[i-j]*a[j]; + } +} + +/*---------------------------------------------------------------------------*\ + + find_aks() + + This function takes a frame of samples, and determines the linear + prediction coefficients for that frame of samples. + +\*---------------------------------------------------------------------------*/ + +void find_aks( + 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[LPC_MAX_N]; /* windowed frame of Nsam speech samples */ + float R[LPC_MAX+1]; /* order+1 autocorrelation values of Sn[] */ + int i; + + assert(order < LPC_MAX); + assert(Nsam < LPC_MAX_N); + + hanning_window(Sn,Wn,Nsam); + autocorrelate(Wn,R,Nsam,order); + 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; +} + +/*---------------------------------------------------------------------------*\ + + weight() + + Weights a vector of LPCs. + +\*---------------------------------------------------------------------------*/ + +void weight( + float ak[], /* vector of order+1 LPCs */ + float gamma, /* weighting factor */ + int order, /* num LPCs (excluding leading 1.0) */ + float akw[] /* weighted vector of order+1 LPCs */ +) +{ + int i; + + for(i=1; i<=order; i++) + akw[i] = ak[i]*pow(gamma,(float)i); +} + diff --git a/gr-codec2-vocoder/src/lib/codec2/lpc.h b/gr-codec2-vocoder/src/lib/codec2/lpc.h new file mode 100644 index 000000000..ead05e1ba --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/lpc.h @@ -0,0 +1,42 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: lpc.h + AUTHOR......: David Rowe + DATE CREATED: 24/8/09 + + Linear Prediction functions written in C. + +\*---------------------------------------------------------------------------*/ + +/* + 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 Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +#ifndef __LPC__ +#define __LPC__ + +#define LPC_MAX_ORDER 20 + +void hanning_window(float Sn[], float Wn[], int Nsam); +void autocorrelate(float Sn[], float Rn[], int Nsam, int order); +void autocorrelate_freq(float Pw[], float w[], float R[], int Nsam, int order); +void levinson_durbin(float R[], float lpcs[], int order); +void inverse_filter(float Sn[], float a[], int Nsam, float res[], int order); +void synthesis_filter(float res[], float a[], int Nsam, int order, float Sn_[]); +void find_aks(float Sn[], float a[], int Nsam, int order, float *E); +void weight(float ak[], float gamma, int order, float akw[]); + +#endif diff --git a/gr-codec2-vocoder/src/lib/codec2/lsp.c b/gr-codec2-vocoder/src/lib/codec2/lsp.c new file mode 100644 index 000000000..47001c1ef --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/lsp.c @@ -0,0 +1,325 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: lsp.c + AUTHOR......: David Rowe + DATE CREATED: 24/2/93 + + + This file contains functions for LPC to LSP conversion and LSP to + LPC conversion. Note that the LSP coefficients are not in radians + format but in the x domain of the unit circle. + +\*---------------------------------------------------------------------------*/ + +/* + 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 Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +#include "defines.h" +#include "lsp.h" +#include <math.h> +#include <stdio.h> +#include <stdlib.h> + +/* Only 10 gets used, so far. */ +#define LSP_MAX_ORDER 20 + +/*---------------------------------------------------------------------------*\ + + 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}. + +\*---------------------------------------------------------------------------*/ + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: cheb_poly_eva() + AUTHOR......: David Rowe + DATE CREATED: 24/2/93 + + This function evalutes a series of chebyshev polynomials + + FIXME: performing memory allocation at run time is very inefficient, + replace with stack variables of MAX_P size. + +\*---------------------------------------------------------------------------*/ + + +static float +cheb_poly_eva(float *coef,float x,int m) +/* float coef[] coefficients of the polynomial to be evaluated */ +/* float x the point where polynomial is to be evaluated */ +/* int m order of the polynomial */ +{ + int i; + float *t,*u,*v,sum; + float T[(LSP_MAX_ORDER / 2) + 1]; + + /* Initialise pointers */ + + t = T; /* T[i-2] */ + *t++ = 1.0; + u = t--; /* T[i-1] */ + *u++ = x; + v = u--; /* T[i] */ + + /* Evaluate chebyshev series formulation using iterative approach */ + + for(i=2;i<=m/2;i++) + *v++ = (2*x)*(*u++) - *t++; /* T[i] = 2*x*T[i-1] - T[i-2] */ + + sum=0.0; /* initialise sum to zero */ + t = T; /* reset pointer */ + + /* Evaluate polynomial and return value also free memory space */ + + for(i=0;i<=m/2;i++) + sum+=coef[(m/2)-i]**t++; + + return sum; +} + + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: lpc_to_lsp() + AUTHOR......: David Rowe + DATE CREATED: 24/2/93 + + This function converts LPC coefficients to LSP coefficients. + +\*---------------------------------------------------------------------------*/ + +int lpc_to_lsp (float *a, int lpcrdr, float *freq, int nb, float delta) +/* float *a lpc coefficients */ +/* int lpcrdr order of LPC coefficients (10) */ +/* float *freq LSP frequencies in radians */ +/* int nb number of sub-intervals (4) */ +/* float delta grid spacing interval (0.02) */ +{ + float psuml,psumr,psumm,temp_xr,xl,xr,xm = 0; + float temp_psumr; + int i,j,m,flag,k; + float *px; /* ptrs of respective P'(z) & Q'(z) */ + float *qx; + float *p; + float *q; + float *pt; /* ptr used for cheb_poly_eval() + whether P' or Q' */ + int roots=0; /* number of roots found */ + float Q[LSP_MAX_ORDER + 1]; + float P[LSP_MAX_ORDER + 1]; + + flag = 1; + m = lpcrdr/2; /* order of P'(z) & Q'(z) polynimials */ + + /* Allocate memory space for polynomials */ + + /* 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; /* initilaise ptrs */ + qx = Q; + p = px; + q = qx; + *px++ = 1.0; + *qx++ = 1.0; + for(i=1;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++; + } + px = P; /* re-initialise ptrs */ + qx = Q; + + /* 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 = 1.0; /* start at point xl = 1 */ + + + for(j=0;j<lpcrdr;j++){ + if(j%2) /* determines whether P' or Q' is eval. */ + pt = qx; + else + pt = px; + + psuml = cheb_poly_eva(pt,xl,lpcrdr); /* evals poly. at xl */ + flag = 1; + while(flag && (xr >= -1.0)){ + xr = xl - delta ; /* interval spacing */ + psumr = cheb_poly_eva(pt,xr,lpcrdr);/* 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((psumr*psuml)<0.0){ + roots++; + + psumm=psuml; + for(k=0;k<=nb;k++){ + xm = (xl+xr)/2; /* bisect the interval */ + psumm=cheb_poly_eva(pt,xm,lpcrdr); + if(psumm*psuml>0.){ + psuml=psumm; + xl=xm; + } + else{ + psumr=psumm; + xr=xm; + } + } + + /* once zero is found, reset initial interval to xr */ + freq[j] = (xm); + xl = xm; + flag = 0; /* reset flag for next search */ + } + else{ + psuml=temp_psumr; + xl=temp_xr; + } + } + } + + /* convert from x domain to radians */ + + for(i=0; i<lpcrdr; i++) { + freq[i] = acos(freq[i]); + } + + return(roots); +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: lsp_to_lpc() + AUTHOR......: David Rowe + DATE CREATED: 24/2/93 + + This function converts LSP coefficients to LPC coefficients. In the + Speex code we worked out a way to simplify this significantly. + +\*---------------------------------------------------------------------------*/ + +void lsp_to_lpc(float *lsp, float *ak, int lpcrdr) +/* float *freq array of LSP frequencies in radians */ +/* float *ak array of LPC coefficients */ +/* int lpcrdr order of LPC coefficients */ + + +{ + int i,j; + float xout1,xout2,xin1,xin2; + float *pw,*n1,*n2,*n3,*n4 = 0; + int m = lpcrdr/2; + float freq[LSP_MAX_ORDER]; + float Wp[(LSP_MAX_ORDER * 4) + 2]; + + /* convert from radians to the x=cos(w) domain */ + + for(i=0; i<lpcrdr; i++) + freq[i] = cos(lsp[i]); + + 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; + + /* 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++){ + for(i=0;i<m;i++){ + n1 = pw+(i*4); + n2 = n1 + 1; + n3 = n2 + 1; + n4 = n3 + 1; + xout1 = xin1 - 2*(freq[2*i]) * *n1 + *n2; + xout2 = xin2 - 2*(freq[2*i+1]) * *n3 + *n4; + *n2 = *n1; + *n4 = *n3; + *n1 = xin1; + *n3 = xin2; + xin1 = xout1; + xin2 = xout2; + } + xout1 = xin1 + *(n4+1); + xout2 = xin2 - *(n4+2); + ak[j] = (xout1 + xout2)*0.5; + *(n4+1) = xin1; + *(n4+2) = xin2; + + xin1 = 0.0; + xin2 = 0.0; + } +} + diff --git a/gr-codec2-vocoder/src/lib/codec2/lsp.h b/gr-codec2-vocoder/src/lib/codec2/lsp.h new file mode 100644 index 000000000..5acef0184 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/lsp.h @@ -0,0 +1,37 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: lsp.c + AUTHOR......: David Rowe + DATE CREATED: 24/2/93 + + + This file contains functions for LPC to LSP conversion and LSP to + LPC conversion. Note that the LSP coefficients are not in radians + format but in the x domain of the unit circle. + +\*---------------------------------------------------------------------------*/ + +/* + 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 Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +#ifndef __LSP__ +#define __LSP__ + +int lpc_to_lsp (float *a, int lpcrdr, float *freq, int nb, float delta); +void lsp_to_lpc(float *freq, float *ak, int lpcrdr); + +#endif diff --git a/gr-codec2-vocoder/src/lib/codec2/nlp.c b/gr-codec2-vocoder/src/lib/codec2/nlp.c new file mode 100644 index 000000000..42ae90919 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/nlp.c @@ -0,0 +1,364 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: nlp.c + AUTHOR......: David Rowe + DATE CREATED: 23/3/93 + + Non Linear Pitch (NLP) estimation functions. + +\*---------------------------------------------------------------------------*/ + +/* + 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 Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +#include "defines.h" +#include "nlp.h" +#include "dump.h" +#include "fft.h" + +#include <assert.h> +#include <math.h> +#include <stdlib.h> + +/*---------------------------------------------------------------------------*\ + + DEFINES + +\*---------------------------------------------------------------------------*/ + +#define PMAX_M 600 /* maximum NLP analysis window size */ +#define COEFF 0.95 /* notch filter parameter */ +#define PE_FFT_SIZE 512 /* DFT size for pitch estimation */ +#define DEC 5 /* decimation factor */ +#define SAMPLE_RATE 8000 +#define PI 3.141592654 /* mathematical constant */ +#define T 0.1 /* threshold for local minima candidate */ +#define F0_MAX 500 +#define CNLP 0.3 /* post processor constant */ +#define NLP_NTAP 48 /* Decimation LPF order */ + +/*---------------------------------------------------------------------------*\ + + GLOBALS + +\*---------------------------------------------------------------------------*/ + +/* 48 tap 600Hz low pass FIR filter coefficients */ + +const float nlp_fir[] = { + -1.0818124e-03, + -1.1008344e-03, + -9.2768838e-04, + -4.2289438e-04, + 5.5034190e-04, + 2.0029849e-03, + 3.7058509e-03, + 5.1449415e-03, + 5.5924666e-03, + 4.3036754e-03, + 8.0284511e-04, + -4.8204610e-03, + -1.1705810e-02, + -1.8199275e-02, + -2.2065282e-02, + -2.0920610e-02, + -1.2808831e-02, + 3.2204775e-03, + 2.6683811e-02, + 5.5520624e-02, + 8.6305944e-02, + 1.1480192e-01, + 1.3674206e-01, + 1.4867556e-01, + 1.4867556e-01, + 1.3674206e-01, + 1.1480192e-01, + 8.6305944e-02, + 5.5520624e-02, + 2.6683811e-02, + 3.2204775e-03, + -1.2808831e-02, + -2.0920610e-02, + -2.2065282e-02, + -1.8199275e-02, + -1.1705810e-02, + -4.8204610e-03, + 8.0284511e-04, + 4.3036754e-03, + 5.5924666e-03, + 5.1449415e-03, + 3.7058509e-03, + 2.0029849e-03, + 5.5034190e-04, + -4.2289438e-04, + -9.2768838e-04, + -1.1008344e-03, + -1.0818124e-03 +}; + +typedef struct { + float sq[PMAX_M]; /* squared speech samples */ + float mem_x,mem_y; /* memory for notch filter */ + float mem_fir[NLP_NTAP]; /* decimation FIR filter memory */ +} NLP; + +float post_process_mbe(COMP Fw[], int pmin, int pmax, float gmax); +float post_process_sub_multiples(COMP Fw[], + int pmin, int pmax, float gmax, int gmax_bin, + float *prev_Wo); + +/*---------------------------------------------------------------------------*\ + + nlp_create() + + Initialisation function for NLP pitch estimator. + +\*---------------------------------------------------------------------------*/ + +void *nlp_create() +{ + NLP *nlp; + int i; + + nlp = (NLP*)malloc(sizeof(NLP)); + if (nlp == NULL) + return NULL; + + for(i=0; i<PMAX_M; i++) + nlp->sq[i] = 0.0; + nlp->mem_x = 0.0; + nlp->mem_y = 0.0; + for(i=0; i<NLP_NTAP; i++) + nlp->mem_fir[i] = 0.0; + + return (void*)nlp; +} + +/*---------------------------------------------------------------------------*\ + + nlp_destory() + + Initialisation function for NLP pitch estimator. + +\*---------------------------------------------------------------------------*/ + +void nlp_destroy(void *nlp_state) +{ + assert(nlp_state != NULL); + free(nlp_state); +} + +/*---------------------------------------------------------------------------*\ + + nlp() + + Determines the pitch in samples using the Non Linear Pitch (NLP) + algorithm [1]. Returns the fundamental in Hz. Note that the actual + pitch estimate is for the centre of the M sample Sn[] vector, not + the current N sample input vector. This is (I think) a delay of 2.5 + frames with N=80 samples. You should align further analysis using + this pitch estimate to be centred on the middle of Sn[]. + + Two post processors have been tried, the MBE version (as discussed + in [1]), and a post processor that checks sub-multiples. Both + suffer occasional gross pitch errors (i.e. neither are perfect). In + the presence of background noise the sub-multiple algorithm tends + towards low F0 which leads to better sounding background noise than + the MBE post processor. + + A good way to test and develop the NLP pitch estimator is using the + tnlp (codec2/unittest) and the codec2/octave/plnlp.m Octave script. + + A pitch tracker searching a few frames forward and backward in time + would be a useful addition. + + References: + + [1] http://www.itr.unisa.edu.au/~steven/thesis/dgr.pdf Chapter 4 + +\*---------------------------------------------------------------------------*/ + +float nlp( + void *nlp_state, + float Sn[], /* input speech vector */ + int n, /* frames shift (no. new samples in Sn[]) */ + int m, /* analysis window size */ + int pmin, /* minimum pitch value */ + int pmax, /* maximum pitch value */ + float *pitch, /* estimated pitch period in samples */ + COMP Sw[], /* Freq domain version of Sn[] */ + float *prev_Wo +) +{ + NLP *nlp; + float notch; /* current notch filter output */ + COMP Fw[PE_FFT_SIZE]; /* DFT of squared signal */ + float gmax; + int gmax_bin; + int i,j; + float best_f0; + + assert(nlp_state != NULL); + nlp = (NLP*)nlp_state; + + /* Square, notch filter at DC, and LP filter vector */ + + for(i=m-n; i<M; i++) /* square latest speech samples */ + nlp->sq[i] = Sn[i]*Sn[i]; + + for(i=m-n; i<m; i++) { /* notch filter at DC */ + notch = nlp->sq[i] - nlp->mem_x; + notch += COEFF*nlp->mem_y; + nlp->mem_x = nlp->sq[i]; + nlp->mem_y = notch; + nlp->sq[i] = notch; + } + + for(i=m-n; i<m; i++) { /* FIR filter vector */ + + for(j=0; j<NLP_NTAP-1; j++) + nlp->mem_fir[j] = nlp->mem_fir[j+1]; + nlp->mem_fir[NLP_NTAP-1] = nlp->sq[i]; + + nlp->sq[i] = 0.0; + for(j=0; j<NLP_NTAP; j++) + nlp->sq[i] += nlp->mem_fir[j]*nlp_fir[j]; + } + + /* Decimate and DFT */ + + for(i=0; i<PE_FFT_SIZE; i++) { + Fw[i].real = 0.0; + Fw[i].imag = 0.0; + } + for(i=0; i<m/DEC; i++) { + Fw[i].real = nlp->sq[i*DEC]*(0.5 - 0.5*cos(2*PI*i/(m/DEC-1))); + } +#ifdef DUMP + dump_dec(Fw); +#endif + fft(&Fw[0].real,PE_FFT_SIZE,1); + for(i=0; i<PE_FFT_SIZE; i++) + Fw[i].real = Fw[i].real*Fw[i].real + Fw[i].imag*Fw[i].imag; + +#ifdef DUMP + dump_sq(nlp->sq); + dump_Fw(Fw); +#endif + + /* find global peak */ + + gmax = 0.0; + gmax_bin = PE_FFT_SIZE*DEC/pmax; + for(i=PE_FFT_SIZE*DEC/pmax; i<=PE_FFT_SIZE*DEC/pmin; i++) { + if (Fw[i].real > gmax) { + gmax = Fw[i].real; + gmax_bin = i; + } + } + + best_f0 = post_process_sub_multiples(Fw, pmin, pmax, gmax, gmax_bin, + prev_Wo); + + /* Shift samples in buffer to make room for new samples */ + + for(i=0; i<m-n; i++) + nlp->sq[i] = nlp->sq[i+n]; + + /* return pitch and F0 estimate */ + + *pitch = (float)SAMPLE_RATE/best_f0; + return(best_f0); +} + +/*---------------------------------------------------------------------------*\ + + post_process_sub_multiples() + + Given the global maximma of Fw[] we search interger submultiples for + local maxima. If local maxima exist and they are above an + experimentally derived threshold (OK a magic number I pulled out of + the air) we choose the submultiple as the F0 estimate. + + The rational for this is that the lowest frequency peak of Fw[] + should be F0, as Fw[] can be considered the autocorrelation function + of Sw[] (the speech spectrum). However sometimes due to phase + effects the lowest frequency maxima may not be the global maxima. + + This works OK in practice and favours low F0 values in the presence + of background noise which means the sinusoidal codec does an OK job + of synthesising the background noise. High F0 in background noise + tends to sound more periodic introducing annoying artifacts. + +\*---------------------------------------------------------------------------*/ + +float post_process_sub_multiples(COMP Fw[], + int pmin, int pmax, float gmax, int gmax_bin, + float *prev_Wo) +{ + int min_bin, cmax_bin; + int mult; + float thresh, best_f0; + int b, bmin, bmax, lmax_bin; + float lmax, cmax; + int prev_f0_bin; + + /* post process estimate by searching submultiples */ + + mult = 2; + min_bin = PE_FFT_SIZE*DEC/pmax; + cmax_bin = gmax_bin; + prev_f0_bin = *prev_Wo*(4000.0/PI)*(PE_FFT_SIZE*DEC)/SAMPLE_RATE; + + while(gmax_bin/mult >= min_bin) { + + b = gmax_bin/mult; /* determine search interval */ + bmin = 0.8*b; + bmax = 1.2*b; + if (bmin < min_bin) + bmin = min_bin; + + /* lower threshold to favour previous frames pitch estimate, + this is a form of pitch tracking */ + + if ((prev_f0_bin > bmin) && (prev_f0_bin < bmax)) + thresh = CNLP*0.5*gmax; + else + thresh = CNLP*gmax; + + lmax = 0; + lmax_bin = bmin; + for (b=bmin; b<=bmax; b++) /* look for maximum in interval */ + if (Fw[b].real > lmax) { + lmax = Fw[b].real; + lmax_bin = b; + } + + if (lmax > thresh) + if ((lmax > Fw[lmax_bin-1].real) && (lmax > Fw[lmax_bin+1].real)) { + cmax = lmax; + cmax_bin = lmax_bin; + } + + mult++; + } + + best_f0 = (float)cmax_bin*SAMPLE_RATE/(PE_FFT_SIZE*DEC); + + return best_f0; +} + diff --git a/gr-codec2-vocoder/src/lib/codec2/nlp.h b/gr-codec2-vocoder/src/lib/codec2/nlp.h new file mode 100644 index 000000000..88a3733dc --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/nlp.h @@ -0,0 +1,39 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: nlp.c + AUTHOR......: David Rowe + DATE CREATED: 23/3/93 + + Non Linear Pitch (NLP) estimation functions. + +\*---------------------------------------------------------------------------*/ + +/* + 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 Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +#ifndef __NLP__ +#define __NLP__ + +#include "comp.h" + +void *nlp_create(); +void nlp_destroy(void *nlp_state); +float nlp(void *nlp_state, float Sn[], int n, int m, int pmin, int pmax, + float *pitch, COMP Sw[], float *prev_Wo); +float test_candidate_mbe(COMP Sw[], float f0, COMP Sw_[]); + +#endif diff --git a/gr-codec2-vocoder/src/lib/codec2/pack.c b/gr-codec2-vocoder/src/lib/codec2/pack.c new file mode 100644 index 000000000..31551dfc4 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/pack.c @@ -0,0 +1,105 @@ +/* + Copyright (C) 2010 Perens LLC <bruce@perens.com> + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + 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, see <http://www.gnu.org/licenses/>. + + */ +#include "defines.h" +#include "quantise.h" +#include <stdio.h> + +/* Compile-time constants */ +/* Size of unsigned char in bits. Assumes 8 bits-per-char. */ +static const unsigned int WordSize = 8; + +/* Mask to pick the bit component out of bitIndex. */ +static const unsigned int IndexMask = 0x7; + +/* Used to pick the word component out of bitIndex. */ +static const unsigned int ShiftRight = 3; + +/** Pack a bit field into a bit string, encoding the field in Gray code. + * + * The output is an array of unsigned char data. The fields are efficiently + * packed into the bit string. The Gray coding is a naive attempt to reduce + * the effect of single-bit errors, we expect to do a better job as the + * codec develops. + * + * This code would be simpler if it just set one bit at a time in the string, + * but would hit the same cache line more often. I'm not sure the complexity + * gains us anything here. + * + * Although field is currently of int type rather than unsigned for + * compatibility with the rest of the code, indices are always expected to + * be >= 0. + */ +void +pack( + unsigned char * bitArray, /* The output bit string. */ + unsigned int * bitIndex, /* Index into the string in BITS, not bytes.*/ + int field, /* The bit field to be packed. */ + unsigned int fieldWidth/* Width of the field in BITS, not bytes. */ + ) +{ + /* Convert the field to Gray code */ + field = (field >> 1) ^ field; + + do { + unsigned int bI = *bitIndex; + unsigned int bitsLeft = WordSize - (bI & IndexMask); + unsigned int sliceWidth = + bitsLeft < fieldWidth ? bitsLeft : fieldWidth; + unsigned int wordIndex = bI >> ShiftRight; + + bitArray[wordIndex] |= + ((unsigned char)((field >> (fieldWidth - sliceWidth)) + << (bitsLeft - sliceWidth))); + + *bitIndex = bI + sliceWidth; + fieldWidth -= sliceWidth; + } while ( fieldWidth != 0 ); +} + +/** Unpack a field from a bit string, converting from Gray code to binary. + * + */ +int +unpack( + const unsigned char * bitArray, /* The input bit string. */ + unsigned int * bitIndex, /* Index into the string in BITS, not bytes.*/ + unsigned int fieldWidth/* Width of the field in BITS, not bytes. */ + ) +{ + unsigned int field = 0; + unsigned int t; + + do { + unsigned int bI = *bitIndex; + unsigned int bitsLeft = WordSize - (bI & IndexMask); + unsigned int sliceWidth = + bitsLeft < fieldWidth ? bitsLeft : fieldWidth; + + field |= (((bitArray[bI >> ShiftRight] >> (bitsLeft - sliceWidth)) & ((1 << sliceWidth) - 1)) << (fieldWidth - sliceWidth)); + + *bitIndex = bI + sliceWidth; + fieldWidth -= sliceWidth; + } while ( fieldWidth != 0 ); + + /* Convert from Gray code to binary. Works for maximum 8-bit fields. */ + t = field ^ (field >> 8); + t ^= (t >> 4); + t ^= (t >> 2); + t ^= (t >> 1); + return t; +} diff --git a/gr-codec2-vocoder/src/lib/codec2/phase.c b/gr-codec2-vocoder/src/lib/codec2/phase.c new file mode 100644 index 000000000..0e1a14a60 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/phase.c @@ -0,0 +1,262 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: phase.c + AUTHOR......: David Rowe + DATE CREATED: 1/2/09 + + Functions for modelling and synthesising phase. + +\*---------------------------------------------------------------------------*/ + +/* + 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 Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not,see <http://www.gnu.org/licenses/>. +*/ + +#include "defines.h" +#include "phase.h" +#include "fft.h" +#include "comp.h" +#include "glottal.c" + +#include <assert.h> +#include <math.h> +#include <string.h> +#include <stdlib.h> + +#define GLOTTAL_FFT_SIZE 512 + +/*---------------------------------------------------------------------------*\ + + aks_to_H() + + Samples the complex LPC synthesis filter spectrum at the harmonic + frequencies. + +\*---------------------------------------------------------------------------*/ + +void aks_to_H( + MODEL *model, /* model parameters */ + float aks[], /* LPC's */ + float G, /* energy term */ + COMP H[], /* complex LPC spectral samples */ + int order +) +{ + COMP Pw[FFT_DEC]; /* power spectrum */ + int i,m; /* loop variables */ + int am,bm; /* limits of current band */ + float r; /* no. rads/bin */ + float Em; /* energy in band */ + float Am; /* spectral amplitude sample */ + int b; /* centre bin of harmonic */ + float phi_; /* phase of LPC spectra */ + + r = TWO_PI/(FFT_DEC); + + /* Determine DFT of A(exp(jw)) ------------------------------------------*/ + + for(i=0; i<FFT_DEC; i++) { + Pw[i].real = 0.0; + Pw[i].imag = 0.0; + } + + for(i=0; i<=order; i++) + Pw[i].real = aks[i]; + + fft(&Pw[0].real,FFT_DEC,-1); + + /* Sample magnitude and phase at harmonics */ + + for(m=1; m<=model->L; m++) { + am = floor((m - 0.5)*model->Wo/r + 0.5); + bm = floor((m + 0.5)*model->Wo/r + 0.5); + b = floor(m*model->Wo/r + 0.5); + + Em = 0.0; + for(i=am; i<bm; i++) + Em += G/(Pw[i].real*Pw[i].real + Pw[i].imag*Pw[i].imag); + Am = sqrt(fabs(Em/(bm-am))); + + phi_ = -atan2(Pw[b].imag,Pw[b].real); + H[m].real = Am*cos(phi_); + H[m].imag = Am*sin(phi_); + } +} + + +/*---------------------------------------------------------------------------*\ + + phase_synth_zero_order() + + Synthesises phases based on SNR and a rule based approach. No phase + parameters are required apart from the SNR (which can be reduced to a + 1 bit V/UV decision per frame). + + The phase of each harmonic is modelled as the phase of a LPC + synthesis filter excited by an impulse. Unlike the first order + model the position of the impulse is not transmitted, so we create + an excitation pulse train using a rule based approach. + + Consider a pulse train with a pulse starting time n=0, with pulses + repeated at a rate of Wo, the fundamental frequency. A pulse train + in the time domain is equivalent to harmonics in the frequency + domain. We can make an excitation pulse train using a sum of + sinsusoids: + + for(m=1; m<=L; m++) + ex[n] = cos(m*Wo*n) + + Note: the Octave script ../octave/phase.m is an example of this if + you would like to try making a pulse train. + + The phase of each excitation harmonic is: + + arg(E[m]) = mWo + + where E[m] are the complex excitation (freq domain) samples, + arg(x), just returns the phase of a complex sample x. + + As we don't transmit the pulse position for this model, we need to + synthesise it. Now the excitation pulses occur at a rate of Wo. + This means the phase of the first harmonic advances by N samples + over a synthesis frame of N samples. For example if Wo is pi/20 + (200 Hz), then over a 10ms frame (N=80 samples), the phase of the + first harmonic would advance (pi/20)*80 = 4*pi or two complete + cycles. + + We generate the excitation phase of the fundamental (first + harmonic): + + arg[E[1]] = Wo*N; + + We then relate the phase of the m-th excitation harmonic to the + phase of the fundamental as: + + arg(E[m]) = m*arg(E[1]) + + This E[m] then gets passed through the LPC synthesis filter to + determine the final harmonic phase. + + Comparing to speech synthesised using original phases: + + - Through headphones speech synthesised with this model is not as + good. Through a loudspeaker it is very close to original phases. + + - If there are voicing errors, the speech can sound clicky or + staticy. If V speech is mistakenly declared UV, this model tends to + synthesise impulses or clicks, as there is usually very little shift or + dispersion through the LPC filter. + + - When combined with LPC amplitude modelling there is an additional + drop in quality. I am not sure why, theory is interformant energy + is raised making any phase errors more obvious. + + NOTES: + + 1/ This synthesis model is effectively the same as a simple LPC-10 + vocoders, and yet sounds much better. Why? Conventional wisdom + (AMBE, MELP) says mixed voicing is required for high quality + speech. + + 2/ I am pretty sure the Lincoln Lab sinusoidal coding guys (like xMBE + also from MIT) first described this zero phase model, I need to look + up the paper. + + 3/ Note that this approach could cause some discontinuities in + the phase at the edge of synthesis frames, as no attempt is made + to make sure that the phase tracks are continuous (the excitation + phases are continuous, but not the final phases after filtering + by the LPC spectra). Technically this is a bad thing. However + this may actually be a good thing, disturbing the phase tracks a + bit. More research needed, e.g. test a synthesis model that adds + a small delta-W to make phase tracks line up for voiced + harmonics. + +\*---------------------------------------------------------------------------*/ + +void phase_synth_zero_order( + MODEL *model, + float aks[], + float *ex_phase, /* excitation phase of fundamental */ + int order +) +{ + int m; + float new_phi; + COMP Ex[MAX_AMP]; /* excitation samples */ + COMP A_[MAX_AMP]; /* synthesised harmonic samples */ + COMP H[MAX_AMP]; /* LPC freq domain samples */ + float G; + float jitter = 0.0; + float r; + int b; + + G = 1.0; + aks_to_H(model, aks, G, H, order); + + /* + Update excitation fundamental phase track, this sets the position + of each pitch pulse during voiced speech. After much experiment + I found that using just this frame's Wo improved quality for UV + sounds compared to interpolating two frames Wo like this: + + ex_phase[0] += (*prev_Wo+mode->Wo)*N/2; + */ + + ex_phase[0] += (model->Wo)*N; + ex_phase[0] -= TWO_PI*floor(ex_phase[0]/TWO_PI + 0.5); + r = TWO_PI/GLOTTAL_FFT_SIZE; + + for(m=1; m<=model->L; m++) { + + /* generate excitation */ + + if (model->voiced) { + /* I think adding a little jitter helps improve low pitch + males like hts1a. This moves the onset of each harmonic + over at +/- 0.25 of a sample. + */ + jitter = 0.25*(1.0 - 2.0*rand()/RAND_MAX); + b = floor(m*model->Wo/r + 0.5); + if (b > ((GLOTTAL_FFT_SIZE/2)-1)) { + b = (GLOTTAL_FFT_SIZE/2)-1; + } + Ex[m].real = cos(ex_phase[0]*m - jitter*model->Wo*m + glottal[b]); + Ex[m].imag = sin(ex_phase[0]*m - jitter*model->Wo*m + glottal[b]); + } + else { + + /* When a few samples were tested I found that LPC filter + phase is not needed in the unvoiced case, but no harm in + keeping it. + */ + float phi = TWO_PI*(float)rand()/RAND_MAX; + Ex[m].real = cos(phi); + Ex[m].imag = sin(phi); + } + + /* filter using LPC filter */ + + A_[m].real = H[m].real*Ex[m].real - H[m].imag*Ex[m].imag; + A_[m].imag = H[m].imag*Ex[m].real + H[m].real*Ex[m].imag; + + /* modify sinusoidal phase */ + + new_phi = atan2(A_[m].imag, A_[m].real+1E-12); + model->phi[m] = new_phi; + } + +} diff --git a/gr-codec2-vocoder/src/lib/codec2/phase.h b/gr-codec2-vocoder/src/lib/codec2/phase.h new file mode 100644 index 000000000..833bc7cdc --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/phase.h @@ -0,0 +1,34 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: phase.h + AUTHOR......: David Rowe + DATE CREATED: 1/2/09 + + Functions for modelling phase. + +\*---------------------------------------------------------------------------*/ + +/* + 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 Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +#ifndef __PHASE__ +#define __PHASE__ + +void phase_synth_zero_order(MODEL *model, float aks[], float *ex_phase, + int order); + +#endif diff --git a/gr-codec2-vocoder/src/lib/codec2/postfilter.c b/gr-codec2-vocoder/src/lib/codec2/postfilter.c new file mode 100644 index 000000000..6e17eeb87 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/postfilter.c @@ -0,0 +1,133 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: postfilter.c + AUTHOR......: David Rowe + DATE CREATED: 13/09/09 + + Postfilter to improve sound quality for speech with high levels of + background noise. Unlike mixed-excitation models requires no bits + to be transmitted to handle background noise. + +\*---------------------------------------------------------------------------*/ + +/* + 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 Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +#include <stdlib.h> +#include <stdio.h> +#include <math.h> + +#include "defines.h" +#include "comp.h" +#include "dump.h" +#include "postfilter.h" + +/*---------------------------------------------------------------------------*\ + + DEFINES + +\*---------------------------------------------------------------------------*/ + +#define BG_THRESH 40.0 /* only consider low levels signals for bg_est */ +#define BG_BETA 0.1 /* averaging filter constant */ + +/*---------------------------------------------------------------------------*\ + + postfilter() + + The post filter is designed to help with speech corrupted by + background noise. The zero phase model tends to make speech with + background noise sound "clicky". With high levels of background + noise the low level inter-formant parts of the spectrum will contain + noise rather than speech harmonics, so modelling them as voiced + (i.e. a continuous, non-random phase track) is inaccurate. + + Some codecs (like MBE) have a mixed voicing model that breaks the + spectrum into voiced and unvoiced regions. Several bits/frame + (5-12) are required to transmit the frequency selective voicing + information. Mixed excitation also requires accurate voicing + estimation (parameter estimators always break occasionally under + exceptional condition). + + In our case we use a post filter approach which requires no + additional bits to be transmitted. The decoder measures the average + level of the background noise during unvoiced frames. If a harmonic + is less than this level it is made unvoiced by randomising it's + phases. + + This idea is rather experimental. Some potential problems that may + happen: + + 1/ If someone says "aaaaaaaahhhhhhhhh" will background estimator track + up to speech level? This would be a bad thing. + + 2/ If background noise suddenly dissapears from the source speech does + estimate drop quickly? What is noise suddenly re-appears? + + 3/ Background noise with a non-flat sepctrum. Current algorithm just + comsiders scpetrum as a whole, but this could be broken up into + bands, each with their own estimator. + + 4/ Males and females with the same level of background noise. Check + performance the same. Changing Wo affects width of each band, may + affect bg energy estimates. + + 5/ Not sure what happens during long periods of voiced speech + e.g. "sshhhhhhh" + +\*---------------------------------------------------------------------------*/ + +void postfilter( + MODEL *model, + float *bg_est +) +{ + int m, uv; + float e; + + /* determine average energy across spectrum */ + + e = 0.0; + for(m=1; m<=model->L; m++) + e += model->A[m]*model->A[m]; + + e = 10.0*log10(e/model->L); + + /* If beneath threhold, update bg estimate. The idea + of the threshold is to prevent updating during high level + speech. */ + + if ((e < BG_THRESH) && !model->voiced) + *bg_est = *bg_est*(1.0 - BG_BETA) + e*BG_BETA; + + /* now mess with phases during voiced frames to make any harmonics + less then our background estimate unvoiced. + */ + + uv = 0; + if (model->voiced) + for(m=1; m<=model->L; m++) + if (20.0*log10(model->A[m]) < *bg_est) { + model->phi[m] = TWO_PI*(float)rand()/RAND_MAX; + uv++; + } + +#ifdef DUMP + dump_bg(e, *bg_est, 100.0*uv/model->L); +#endif + +} diff --git a/gr-codec2-vocoder/src/lib/codec2/postfilter.h b/gr-codec2-vocoder/src/lib/codec2/postfilter.h new file mode 100644 index 000000000..bf080b1b6 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/postfilter.h @@ -0,0 +1,33 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: postfilter.h + AUTHOR......: David Rowe + DATE CREATED: 13/09/09 + + Postfilter header file. + +\*---------------------------------------------------------------------------*/ + +/* + 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 Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +#ifndef __POSTFILTER__ +#define __POSTFILTER__ + +void postfilter(MODEL *model, float *bg_est); + +#endif diff --git a/gr-codec2-vocoder/src/lib/codec2/quantise.c b/gr-codec2-vocoder/src/lib/codec2/quantise.c new file mode 100644 index 000000000..ff8d156b5 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/quantise.c @@ -0,0 +1,851 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: quantise.c + AUTHOR......: David Rowe + DATE CREATED: 31/5/92 + + Quantisation functions for the sinusoidal coder. + +\*---------------------------------------------------------------------------*/ + +/* + All rights reserved. + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. + +*/ + +#include <assert.h> +#include <ctype.h> +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <math.h> + +#include "defines.h" +#include "dump.h" +#include "quantise.h" +#include "lpc.h" +#include "lsp.h" +#include "fft.h" + +#define LSP_DELTA1 0.01 /* grid spacing for LSP root searches */ + +/*---------------------------------------------------------------------------*\ + + FUNCTION HEADERS + +\*---------------------------------------------------------------------------*/ + +float speech_to_uq_lsps(float lsp[], float ak[], float Sn[], float w[], + int order); + +/*---------------------------------------------------------------------------*\ + + FUNCTIONS + +\*---------------------------------------------------------------------------*/ + +int lsp_bits(int i) { + return lsp_cb[i].log2m; +} + +#if VECTOR_QUANTISATION +/*---------------------------------------------------------------------------*\ + + quantise_uniform + + Simulates uniform quantising of a float. + +\*---------------------------------------------------------------------------*/ + +void quantise_uniform(float *val, float min, float max, int bits) +{ + int levels = 1 << (bits-1); + float norm; + int index; + + /* hard limit to quantiser range */ + + printf("min: %f max: %f val: %f ", min, max, val[0]); + if (val[0] < min) val[0] = min; + if (val[0] > max) val[0] = max; + + norm = (*val - min)/(max-min); + printf("%f norm: %f ", val[0], norm); + index = fabs(levels*norm + 0.5); + + *val = min + index*(max-min)/levels; + + printf("index %d val_: %f\n", index, val[0]); +} + +#endif + +/*---------------------------------------------------------------------------*\ + + quantise_init + + Loads the entire LSP quantiser comprised of several vector quantisers + (codebooks). + +\*---------------------------------------------------------------------------*/ + +void quantise_init() +{ +} + +/*---------------------------------------------------------------------------*\ + + quantise + + Quantises vec by choosing the nearest vector in codebook cb, and + returns the vector index. The squared error of the quantised vector + is added to se. + +\*---------------------------------------------------------------------------*/ + +long quantise(const float * cb, float vec[], float w[], int k, int m, float *se) +/* float cb[][K]; current VQ codebook */ +/* float vec[]; vector to quantise */ +/* float w[]; weighting vector */ +/* int k; dimension of vectors */ +/* int m; size of codebook */ +/* float *se; accumulated squared error */ +{ + float e; /* current error */ + long besti; /* best index so far */ + float beste; /* best error so far */ + long j; + int i; + + besti = 0; + beste = 1E32; + for(j=0; j<m; j++) { + e = 0.0; + for(i=0; i<k; i++) + e += pow((cb[j*k+i]-vec[i])*w[i],2.0); + if (e < beste) { + beste = e; + besti = j; + } + } + + *se += beste; + + return(besti); +} + +/*---------------------------------------------------------------------------*\ + + lspd_quantise + + Scalar lsp difference quantiser. + +\*---------------------------------------------------------------------------*/ + +void lspd_quantise( + float lsp[], + float lsp_[], + int order +) +{ + int i,k,m; + float lsp_hz[LPC_MAX]; + float lsp__hz[LPC_MAX]; + float dlsp[LPC_MAX]; + float dlsp_[LPC_MAX]; + float wt[1]; + const float *cb; + float se; + int indexes[LPC_MAX]; + + /* convert from radians to Hz so we can use human readable + frequencies */ + + for(i=0; i<order; i++) + lsp_hz[i] = (4000.0/PI)*lsp[i]; + + dlsp[0] = lsp_hz[0]; + for(i=1; i<order; i++) + dlsp[i] = lsp_hz[i] - lsp_hz[i-1]; + + /* simple uniform scalar quantisers */ + + wt[0] = 1.0; + for(i=0; i<order; i++) { + if (i) + dlsp[i] = lsp_hz[i] - lsp__hz[i-1]; + else + dlsp[0] = lsp_hz[0]; + + k = lsp_cbd[i].k; + m = lsp_cbd[i].m; + cb = lsp_cbd[i].cb; + indexes[i] = quantise(cb, &dlsp[i], wt, k, m, &se); + dlsp_[i] = cb[indexes[i]*k]; + + if (i) + lsp__hz[i] = lsp__hz[i-1] + dlsp_[i]; + else + lsp__hz[0] = dlsp_[0]; + } + for(; i<order; i++) + lsp__hz[i] = lsp__hz[i-1] + dlsp[i]; + + /* convert back to radians */ + + for(i=0; i<order; i++) + lsp_[i] = (PI/4000.0)*lsp__hz[i]; +} + +/*---------------------------------------------------------------------------*\ + + lspd_vq_quantise + + Vector lsp difference quantiser. + +\*---------------------------------------------------------------------------*/ + +void lspdvq_quantise( + float lsp[], + float lsp_[], + int order +) +{ + int i,k,m,ncb, nlsp; + float dlsp[LPC_MAX]; + float dlsp_[LPC_MAX]; + float wt[LPC_ORD]; + const float *cb; + float se; + int index; + + dlsp[0] = lsp[0]; + for(i=1; i<order; i++) + dlsp[i] = lsp[i] - lsp[i-1]; + + for(i=0; i<order; i++) + dlsp_[i] = dlsp[i]; + + for(i=0; i<order; i++) + wt[i] = 1.0; + + /* scalar quantise dLSPs 1,2,3,4,5 */ + + for(i=0; i<5; i++) { + if (i) + dlsp[i] = (lsp[i] - lsp_[i-1])*4000.0/PI; + else + dlsp[0] = lsp[0]*4000.0/PI; + + k = lsp_cbdvq[i].k; + m = lsp_cbdvq[i].m; + cb = lsp_cbdvq[i].cb; + index = quantise(cb, &dlsp[i], wt, k, m, &se); + dlsp_[i] = cb[index*k]*PI/4000.0; + + if (i) + lsp_[i] = lsp_[i-1] + dlsp_[i]; + else + lsp_[0] = dlsp_[0]; + } + dlsp[i] = lsp[i] - lsp_[i-1]; + dlsp_[i] = dlsp[i]; + + //printf("lsp[0] %f lsp_[0] %f\n", lsp[0], lsp_[0]); + //printf("lsp[1] %f lsp_[1] %f\n", lsp[1], lsp_[1]); + +#ifdef TT + /* VQ dLSPs 3,4,5 */ + + ncb = 2; + nlsp = 2; + k = lsp_cbdvq[ncb].k; + m = lsp_cbdvq[ncb].m; + cb = lsp_cbdvq[ncb].cb; + index = quantise(cb, &dlsp[nlsp], wt, k, m, &se); + dlsp_[nlsp] = cb[index*k]; + dlsp_[nlsp+1] = cb[index*k+1]; + dlsp_[nlsp+2] = cb[index*k+2]; + + lsp_[0] = dlsp_[0]; + for(i=1; i<5; i++) + lsp_[i] = lsp_[i-1] + dlsp_[i]; + dlsp[i] = lsp[i] - lsp_[i-1]; + dlsp_[i] = dlsp[i]; +#endif + /* VQ dLSPs 6,7,8,9,10 */ + + ncb = 5; + nlsp = 5; + k = lsp_cbdvq[ncb].k; + m = lsp_cbdvq[ncb].m; + cb = lsp_cbdvq[ncb].cb; + index = quantise(cb, &dlsp[nlsp], wt, k, m, &se); + dlsp_[nlsp] = cb[index*k]; + dlsp_[nlsp+1] = cb[index*k+1]; + dlsp_[nlsp+2] = cb[index*k+2]; + dlsp_[nlsp+3] = cb[index*k+3]; + dlsp_[nlsp+4] = cb[index*k+4]; + + /* rebuild LSPs for dLSPs */ + + lsp_[0] = dlsp_[0]; + for(i=1; i<order; i++) + lsp_[i] = lsp_[i-1] + dlsp_[i]; +} + +void check_lsp_order(float lsp[], int lpc_order) +{ + int i; + float tmp; + + for(i=1; i<lpc_order; i++) + if (lsp[i] < lsp[i-1]) { + printf("swap %d\n",i); + tmp = lsp[i-1]; + lsp[i-1] = lsp[i]-0.05; + lsp[i] = tmp+0.05; + } +} + +void force_min_lsp_dist(float lsp[], int lpc_order) +{ + int i; + + for(i=1; i<lpc_order; i++) + if ((lsp[i]-lsp[i-1]) < 0.01) { + lsp[i] += 0.01; + } +} + +/*---------------------------------------------------------------------------*\ + + lpc_model_amplitudes + + Derive a LPC model for amplitude samples then estimate amplitude samples + from this model with optional LSP quantisation. + + Returns the spectral distortion for this frame. + +\*---------------------------------------------------------------------------*/ + +float lpc_model_amplitudes( + float Sn[], /* Input frame of speech samples */ + float w[], + MODEL *model, /* sinusoidal model parameters */ + int order, /* LPC model order */ + int lsp_quant, /* optional LSP quantisation if non-zero */ + float ak[] /* output aks */ +) +{ + float Wn[M]; + float R[LPC_MAX+1]; + float E; + int i,j; + float snr; + float lsp[LPC_MAX]; + float lsp_hz[LPC_MAX]; + float lsp_[LPC_MAX]; + int roots; /* number of LSP roots found */ + int index; + float se; + int k,m; + const float * cb; + float wt[LPC_MAX]; + + for(i=0; i<M; i++) + Wn[i] = Sn[i]*w[i]; + autocorrelate(Wn,R,M,order); + levinson_durbin(R,ak,order); + + E = 0.0; + for(i=0; i<=order; i++) + E += ak[i]*R[i]; + + for(i=0; i<order; i++) + wt[i] = 1.0; + + if (lsp_quant) { + roots = lpc_to_lsp(ak, order, lsp, 5, LSP_DELTA1); + if (roots != order) + printf("LSP roots not found\n"); + + /* convert from radians to Hz to make quantisers more + human readable */ + + for(i=0; i<order; i++) + lsp_hz[i] = (4000.0/PI)*lsp[i]; + + /* simple uniform scalar quantisers */ + + for(i=0; i<10; i++) { + k = lsp_cb[i].k; + m = lsp_cb[i].m; + cb = lsp_cb[i].cb; + index = quantise(cb, &lsp_hz[i], wt, k, m, &se); + lsp_hz[i] = cb[index*k]; + } + + /* experiment: simulating uniform quantisation error + for(i=0; i<order; i++) + lsp[i] += PI*(12.5/4000.0)*(1.0 - 2.0*(float)rand()/RAND_MAX); + */ + + for(i=0; i<order; i++) + lsp[i] = (PI/4000.0)*lsp_hz[i]; + + /* Bandwidth Expansion (BW). Prevents any two LSPs getting too + close together after quantisation. We know from experiment + that LSP quantisation errors < 12.5Hz (25Hz setp size) are + inaudible so we use that as the minimum LSP separation. + */ + + for(i=1; i<5; i++) { + if (lsp[i] - lsp[i-1] < PI*(12.5/4000.0)) + lsp[i] = lsp[i-1] + PI*(12.5/4000.0); + } + + /* as quantiser gaps increased, larger BW expansion was required + to prevent twinkly noises */ + + for(i=5; i<8; i++) { + if (lsp[i] - lsp[i-1] < PI*(25.0/4000.0)) + lsp[i] = lsp[i-1] + PI*(25.0/4000.0); + } + for(i=8; i<order; i++) { + if (lsp[i] - lsp[i-1] < PI*(75.0/4000.0)) + lsp[i] = lsp[i-1] + PI*(75.0/4000.0); + } + + for(j=0; j<order; j++) + lsp_[j] = lsp[j]; + + lsp_to_lpc(lsp_, ak, order); +#ifdef DUMP + dump_lsp(lsp); +#endif + } + +#ifdef DUMP + dump_E(E); +#endif + #ifdef SIM_QUANT + /* simulated LPC energy quantisation */ + { + float e = 10.0*log10(E); + e += 2.0*(1.0 - 2.0*(float)rand()/RAND_MAX); + E = pow(10.0,e/10.0); + } + #endif + + aks_to_M2(ak,order,model,E,&snr, 1); /* {ak} -> {Am} LPC decode */ + + return snr; +} + +/*---------------------------------------------------------------------------*\ + + aks_to_M2() + + Transforms the linear prediction coefficients to spectral amplitude + samples. This function determines A(m) from the average energy per + band using an FFT. + +\*---------------------------------------------------------------------------*/ + +void aks_to_M2( + float ak[], /* LPC's */ + int order, + MODEL *model, /* sinusoidal model parameters for this frame */ + float E, /* energy term */ + float *snr, /* signal to noise ratio for this frame in dB */ + int dump /* true to dump sample to dump file */ +) +{ + COMP Pw[FFT_DEC]; /* power spectrum */ + int i,m; /* loop variables */ + int am,bm; /* limits of current band */ + float r; /* no. rads/bin */ + float Em; /* energy in band */ + float Am; /* spectral amplitude sample */ + float signal, noise; + + r = TWO_PI/(FFT_DEC); + + /* Determine DFT of A(exp(jw)) --------------------------------------------*/ + + for(i=0; i<FFT_DEC; i++) { + Pw[i].real = 0.0; + Pw[i].imag = 0.0; + } + + for(i=0; i<=order; i++) + Pw[i].real = ak[i]; + fft(&Pw[0].real,FFT_DEC,1); + + /* Determine power spectrum P(w) = E/(A(exp(jw))^2 ------------------------*/ + + for(i=0; i<FFT_DEC/2; i++) + Pw[i].real = E/(Pw[i].real*Pw[i].real + Pw[i].imag*Pw[i].imag); +#ifdef DUMP + if (dump) + dump_Pw(Pw); +#endif + + /* Determine magnitudes by linear interpolation of P(w) -------------------*/ + + signal = noise = 0.0; + for(m=1; m<=model->L; m++) { + am = floor((m - 0.5)*model->Wo/r + 0.5); + bm = floor((m + 0.5)*model->Wo/r + 0.5); + Em = 0.0; + + for(i=am; i<bm; i++) + Em += Pw[i].real; + Am = sqrt(Em); + + signal += pow(model->A[m],2.0); + noise += pow(model->A[m] - Am,2.0); + model->A[m] = Am; + } + *snr = 10.0*log10(signal/noise); +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: encode_Wo() + AUTHOR......: David Rowe + DATE CREATED: 22/8/2010 + + Encodes Wo using a WO_LEVELS quantiser. + +\*---------------------------------------------------------------------------*/ + +int encode_Wo(float Wo) +{ + int index; + float Wo_min = TWO_PI/P_MAX; + float Wo_max = TWO_PI/P_MIN; + float norm; + + norm = (Wo - Wo_min)/(Wo_max - Wo_min); + index = floor(WO_LEVELS * norm + 0.5); + if (index < 0 ) index = 0; + if (index > (WO_LEVELS-1)) index = WO_LEVELS-1; + + return index; +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: decode_Wo() + AUTHOR......: David Rowe + DATE CREATED: 22/8/2010 + + Decodes Wo using a WO_LEVELS quantiser. + +\*---------------------------------------------------------------------------*/ + +float decode_Wo(int index) +{ + float Wo_min = TWO_PI/P_MAX; + float Wo_max = TWO_PI/P_MIN; + float step; + float Wo; + + step = (Wo_max - Wo_min)/WO_LEVELS; + Wo = Wo_min + step*(index); + + return Wo; +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: speech_to_uq_lsps() + AUTHOR......: David Rowe + DATE CREATED: 22/8/2010 + + Analyse a windowed frame of time domain speech to determine LPCs + which are the converted to LSPs for quantisation and transmission + over the channel. + +\*---------------------------------------------------------------------------*/ + +float speech_to_uq_lsps(float lsp[], + float ak[], + float Sn[], + float w[], + int order +) +{ + int i, roots; + float Wn[M]; + float R[LPC_MAX+1]; + float E; + + for(i=0; i<M; i++) + Wn[i] = Sn[i]*w[i]; + autocorrelate(Wn, R, M, order); + levinson_durbin(R, ak, order); + + E = 0.0; + for(i=0; i<=order; i++) + E += ak[i]*R[i]; + + roots = lpc_to_lsp(ak, order, lsp, 5, LSP_DELTA1); + if (roots != order) { + /* for some reason LSP roots could not be found */ + /* some alpha testers are reporting this condition */ + fprintf(stderr, "LSP roots not found!\nroots = %d\n", roots); + for(i=0; i<=order; i++) + fprintf(stderr, "a[%d] = %f\n", i, ak[i]); + + /* some benign LSP values we can use instead */ + for(i=0; i<order; i++) + lsp[i] = (PI/order)*(float)i; + } + + return E; +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: encode_lsps() + AUTHOR......: David Rowe + DATE CREATED: 22/8/2010 + + From a vector of unquantised (floating point) LSPs finds the quantised + LSP indexes. + +\*---------------------------------------------------------------------------*/ + +void encode_lsps(int indexes[], float lsp[], int order) +{ + int i,k,m; + float wt[1]; + float lsp_hz[LPC_MAX]; + const float * cb; + float se; + + /* convert from radians to Hz so we can use human readable + frequencies */ + + for(i=0; i<order; i++) + lsp_hz[i] = (4000.0/PI)*lsp[i]; + + /* simple uniform scalar quantisers */ + + wt[0] = 1.0; + for(i=0; i<order; i++) { + k = lsp_cb[i].k; + m = lsp_cb[i].m; + cb = lsp_cb[i].cb; + indexes[i] = quantise(cb, &lsp_hz[i], wt, k, m, &se); + } +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: decode_lsps() + AUTHOR......: David Rowe + DATE CREATED: 22/8/2010 + + From a vector of quantised LSP indexes, returns the quantised + (floating point) LSPs. + +\*---------------------------------------------------------------------------*/ + +void decode_lsps(float lsp[], int indexes[], int order) +{ + int i,k; + float lsp_hz[LPC_MAX]; + const float * cb; + + for(i=0; i<order; i++) { + k = lsp_cb[i].k; + cb = lsp_cb[i].cb; + lsp_hz[i] = cb[indexes[i]*k]; + } + + /* convert back to radians */ + + for(i=0; i<order; i++) + lsp[i] = (PI/4000.0)*lsp_hz[i]; +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: bw_expand_lsps() + AUTHOR......: David Rowe + DATE CREATED: 22/8/2010 + + Applies Bandwidth Expansion (BW) to a vector of LSPs. Prevents any + two LSPs getting too close together after quantisation. We know + from experiment that LSP quantisation errors < 12.5Hz (25Hz setp + size) are inaudible so we use that as the minimum LSP separation. + +\*---------------------------------------------------------------------------*/ + +void bw_expand_lsps(float lsp[], + int order +) +{ + int i; + + for(i=1; i<5; i++) { + if (lsp[i] - lsp[i-1] < PI*(12.5/4000.0)) + lsp[i] = lsp[i-1] + PI*(12.5/4000.0); + } + + /* As quantiser gaps increased, larger BW expansion was required + to prevent twinkly noises. This may need more experiment for + different quanstisers. + */ + + for(i=5; i<8; i++) { + if (lsp[i] - lsp[i-1] < PI*(25.0/4000.0)) + lsp[i] = lsp[i-1] + PI*(25.0/4000.0); + } + for(i=8; i<order; i++) { + if (lsp[i] - lsp[i-1] < PI*(75.0/4000.0)) + lsp[i] = lsp[i-1] + PI*(75.0/4000.0); + } +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: apply_lpc_correction() + AUTHOR......: David Rowe + DATE CREATED: 22/8/2010 + + Apply first harmonic LPC correction at decoder. This helps improve + low pitch males after LPC modelling, like hts1a and morig. + +\*---------------------------------------------------------------------------*/ + +void apply_lpc_correction(MODEL *model) +{ + if (model->Wo < (PI*150.0/4000)) { + model->A[1] *= 0.032; + } +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: encode_energy() + AUTHOR......: David Rowe + DATE CREATED: 22/8/2010 + + Encodes LPC energy using an E_LEVELS quantiser. + +\*---------------------------------------------------------------------------*/ + +int encode_energy(float e) +{ + int index; + float e_min = E_MIN_DB; + float e_max = E_MAX_DB; + float norm; + + e = 10.0*log10(e); + norm = (e - e_min)/(e_max - e_min); + index = floor(E_LEVELS * norm + 0.5); + if (index < 0 ) index = 0; + if (index > (E_LEVELS-1)) index = E_LEVELS-1; + + return index; +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: decode_energy() + AUTHOR......: David Rowe + DATE CREATED: 22/8/2010 + + Decodes energy using a WO_BITS quantiser. + +\*---------------------------------------------------------------------------*/ + +float decode_energy(int index) +{ + float e_min = E_MIN_DB; + float e_max = E_MAX_DB; + float step; + float e; + + step = (e_max - e_min)/E_LEVELS; + e = e_min + step*(index); + e = pow(10.0,e/10.0); + + return e; +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: encode_amplitudes() + AUTHOR......: David Rowe + DATE CREATED: 22/8/2010 + + Time domain LPC is used model the amplitudes which are then + converted to LSPs and quantised. So we don't actually encode the + amplitudes directly, rather we derive an equivalent representation + from the time domain speech. + +\*---------------------------------------------------------------------------*/ + +void encode_amplitudes(int lsp_indexes[], + int *energy_index, + MODEL *model, + float Sn[], + float w[]) +{ + float lsps[LPC_ORD]; + float ak[LPC_ORD+1]; + float e; + + e = speech_to_uq_lsps(lsps, ak, Sn, w, LPC_ORD); + encode_lsps(lsp_indexes, lsps, LPC_ORD); + *energy_index = encode_energy(e); +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: decode_amplitudes() + AUTHOR......: David Rowe + DATE CREATED: 22/8/2010 + + Given the amplitude quantiser indexes recovers the harmonic + amplitudes. + +\*---------------------------------------------------------------------------*/ + +float decode_amplitudes(MODEL *model, + float ak[], + int lsp_indexes[], + int energy_index, + float lsps[], + float *e +) +{ + float snr; + + decode_lsps(lsps, lsp_indexes, LPC_ORD); + bw_expand_lsps(lsps, LPC_ORD); + lsp_to_lpc(lsps, ak, LPC_ORD); + *e = decode_energy(energy_index); + aks_to_M2(ak, LPC_ORD, model, *e, &snr, 1); + apply_lpc_correction(model); + + return snr; +} diff --git a/gr-codec2-vocoder/src/lib/codec2/quantise.h b/gr-codec2-vocoder/src/lib/codec2/quantise.h new file mode 100644 index 000000000..90a3661ff --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/quantise.h @@ -0,0 +1,83 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: quantise.h + AUTHOR......: David Rowe + DATE CREATED: 31/5/92 + + Quantisation functions for the sinusoidal coder. + +\*---------------------------------------------------------------------------*/ + +/* + All rights reserved. + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +#ifndef __QUANTISE__ +#define __QUANTISE__ + +#define WO_BITS 7 +#define WO_LEVELS (1<<WO_BITS) +#define E_BITS 5 +#define E_LEVELS (1<<E_BITS) +#define E_MIN_DB -10.0 +#define E_MAX_DB 40.0 + +void quantise_init(); +float lpc_model_amplitudes(float Sn[], float w[], MODEL *model, int order, + int lsp,float ak[]); +void aks_to_M2(float ak[], int order, MODEL *model, float E, float *snr, + int dump); + +int encode_Wo(float Wo); +float decode_Wo(int index); + +void encode_lsps(int indexes[], float lsp[], int order); +void decode_lsps(float lsp[], int indexes[], int order); +void lspd_quantise(float lsp[], float lsp_[], int order); +void lspdvq_quantise(float lsp[], float lsp_[], int order); + +int encode_energy(float e); +float decode_energy(int index); + +void encode_amplitudes(int lsp_indexes[], + int *energy_index, + MODEL *model, + float Sn[], + float w[]); + +float decode_amplitudes(MODEL *model, + float ak[], + int lsp_indexes[], + int energy_index, + float lsps[], + float *e); + +void pack(unsigned char * bits, unsigned int *nbit, int index, unsigned int index_bits); +int unpack(const unsigned char * bits, unsigned int *nbit, unsigned int index_bits); + +int lsp_bits(int i); + +void apply_lpc_correction(MODEL *model); +float speech_to_uq_lsps(float lsp[], + float ak[], + float Sn[], + float w[], + int order + ); +void bw_expand_lsps(float lsp[], + int order + ); +void decode_lsps(float lsp[], int indexes[], int order); + +#endif diff --git a/gr-codec2-vocoder/src/lib/codec2/sim.sh b/gr-codec2-vocoder/src/lib/codec2/sim.sh new file mode 100755 index 000000000..10152d979 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/sim.sh @@ -0,0 +1,22 @@ +#!/bin/sh +# sim.sh +# David Rowe 10 Sep 2009 + +# Process a source file using the codec 2 simulation. An output +# speech file is generated for each major processing step, from the +# unquantised siusoidal model to fully quantised. This way we can +# listen to the effect of each processing step. Use listensim.sh to +# test the output files. + +../src/c2sim ../raw/$1.raw -o $1_uq.raw +../src/c2sim ../raw/$1.raw --phase0 -o $1_phase0.raw --postfilter +../src/c2sim ../raw/$1.raw --lpc 10 -o $1_lpc10.raw --postfilter +../src/c2sim ../raw/$1.raw --phase0 --lpc 10 -o $1_phase0_lpc10.raw --postfilter +../src/c2sim ../raw/$1.raw --phase0 --lpc 10 --dec -o $1_phase0_lpc10_dec.raw --postfilter +../src/c2sim ../raw/$1.raw --phase0 --lpc 10 --lsp --dec -o $1_phase0_lsp_dec.raw --postfilter + +#../src/c2sim ../raw/$1.raw --lpc 10 --lsp -o $1_lsp.raw +#../src/c2sim ../raw/$1.raw --phase0 --lpc 10 -o $1_phase0_lpc10.raw --postfilter +#../src/c2sim ../raw/$1.raw --phase0 --lpc 10 --lsp -o $1_phase0_lsp.raw --postfilter +#../src/c2sim ../raw/$1.raw --phase0 --lpc 10 --lsp -o $1_phase0_lsp_dec.raw --postfilter --dec + diff --git a/gr-codec2-vocoder/src/lib/codec2/sine.c b/gr-codec2-vocoder/src/lib/codec2/sine.c new file mode 100644 index 000000000..45cc9de71 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/sine.c @@ -0,0 +1,638 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: sine.c + AUTHOR......: David Rowe + DATE CREATED: 19/8/2010 + + Sinusoidal analysis and synthesis functions. + +\*---------------------------------------------------------------------------*/ + +/* + Copyright (C) 1990-2010 David Rowe + + All rights reserved. + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +/*---------------------------------------------------------------------------*\ + + INCLUDES + +\*---------------------------------------------------------------------------*/ + +#include <stdlib.h> +#include <stdio.h> +#include <math.h> + +#include "defines.h" +#include "sine.h" +#include "fft.h" + +#define HPF_BETA 0.125 + +/*---------------------------------------------------------------------------*\ + + HEADERS + +\*---------------------------------------------------------------------------*/ + +void hs_pitch_refinement(MODEL *model, COMP Sw[], float pmin, float pmax, + float pstep); + +/*---------------------------------------------------------------------------*\ + + FUNCTIONS + +\*---------------------------------------------------------------------------*/ + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: make_analysis_window + AUTHOR......: David Rowe + DATE CREATED: 11/5/94 + + Init function that generates the time domain analysis window and it's DFT. + +\*---------------------------------------------------------------------------*/ + +void make_analysis_window(float w[],COMP W[]) +{ + float m; + COMP temp; + int i,j; + + /* + Generate Hamming window centered on M-sample pitch analysis window + + 0 M/2 M-1 + |-------------|-------------| + |-------|-------| + NW samples + + All our analysis/synthsis is centred on the M/2 sample. + */ + + m = 0.0; + for(i=0; i<M/2-NW/2; i++) + w[i] = 0.0; + for(i=M/2-NW/2,j=0; i<M/2+NW/2; i++,j++) { + w[i] = 0.5 - 0.5*cos(TWO_PI*j/(NW-1)); + m += w[i]*w[i]; + } + for(i=M/2+NW/2; i<M; i++) + w[i] = 0.0; + + /* Normalise - makes freq domain amplitude estimation straight + forward */ + + m = 1.0/sqrt(m*FFT_ENC); + for(i=0; i<M; i++) { + w[i] *= m; + } + + /* + Generate DFT of analysis window, used for later processing. Note + we modulo FFT_ENC shift the time domain window w[], this makes the + imaginary part of the DFT W[] equal to zero as the shifted w[] is + even about the n=0 time axis if NW is odd. Having the imag part + of the DFT W[] makes computation easier. + + 0 FFT_ENC-1 + |-------------------------| + + ----\ /---- + \ / + \ / <- shifted version of window w[n] + \ / + \ / + ------- + + |---------| |---------| + NW/2 NW/2 + */ + + for(i=0; i<FFT_ENC; i++) { + W[i].real = 0.0; + W[i].imag = 0.0; + } + for(i=0; i<NW/2; i++) + W[i].real = w[i+M/2]; + for(i=FFT_ENC-NW/2,j=M/2-NW/2; i<FFT_ENC; i++,j++) + W[i].real = w[j]; + + fft(&W[0].real,FFT_ENC,-1); /* "Numerical Recipes in C" FFT */ + + /* + Re-arrange W[] to be symmetrical about FFT_ENC/2. Makes later + analysis convenient. + + Before: + + + 0 FFT_ENC-1 + |----------|---------| + __ _ + \ / + \_______________/ + + After: + + 0 FFT_ENC-1 + |----------|---------| + ___ + / \ + ________/ \_______ + + */ + + + for(i=0; i<FFT_ENC/2; i++) { + temp.real = W[i].real; + temp.imag = W[i].imag; + W[i].real = W[i+FFT_ENC/2].real; + W[i].imag = W[i+FFT_ENC/2].imag; + W[i+FFT_ENC/2].real = temp.real; + W[i+FFT_ENC/2].imag = temp.imag; + } + +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: hpf + AUTHOR......: David Rowe + DATE CREATED: 16 Nov 2010 + + High pass filter with a -3dB point of about 160Hz. + + y(n) = -HPF_BETA*y(n-1) + x(n) - x(n-1) + +\*---------------------------------------------------------------------------*/ + +float hpf(float x, float states[]) +{ + states[0] += -HPF_BETA*states[0] + x - states[1]; + states[1] = x; + + return states[0]; +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: dft_speech + AUTHOR......: David Rowe + DATE CREATED: 27/5/94 + + Finds the DFT of the current speech input speech frame. + +\*---------------------------------------------------------------------------*/ + +void dft_speech(COMP Sw[], float Sn[], float w[]) +{ + int i; + + for(i=0; i<FFT_ENC; i++) { + Sw[i].real = 0.0; + Sw[i].imag = 0.0; + } + + /* Centre analysis window on time axis, we need to arrange input + to FFT this way to make FFT phases correct */ + + /* move 2nd half to start of FFT input vector */ + + for(i=0; i<NW/2; i++) + Sw[i].real = Sn[i+M/2]*w[i+M/2]; + + /* move 1st half to end of FFT input vector */ + + for(i=0; i<NW/2; i++) + Sw[FFT_ENC-NW/2+i].real = Sn[i+M/2-NW/2]*w[i+M/2-NW/2]; + + fft(&Sw[0].real,FFT_ENC,-1); +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: two_stage_pitch_refinement + AUTHOR......: David Rowe + DATE CREATED: 27/5/94 + + Refines the current pitch estimate using the harmonic sum pitch + estimation technique. + +\*---------------------------------------------------------------------------*/ + +void two_stage_pitch_refinement(MODEL *model, COMP Sw[]) +{ + float pmin,pmax,pstep; /* pitch refinment minimum, maximum and step */ + + /* Coarse refinement */ + + pmax = TWO_PI/model->Wo + 5; + pmin = TWO_PI/model->Wo - 5; + pstep = 1.0; + hs_pitch_refinement(model,Sw,pmin,pmax,pstep); + + /* Fine refinement */ + + pmax = TWO_PI/model->Wo + 1; + pmin = TWO_PI/model->Wo - 1; + pstep = 0.25; + hs_pitch_refinement(model,Sw,pmin,pmax,pstep); + + /* Limit range */ + + if (model->Wo < TWO_PI/P_MAX) + model->Wo = TWO_PI/P_MAX; + if (model->Wo > TWO_PI/P_MIN) + model->Wo = TWO_PI/P_MIN; + + model->L = floor(PI/model->Wo); +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: hs_pitch_refinement + AUTHOR......: David Rowe + DATE CREATED: 27/5/94 + + Harmonic sum pitch refinement function. + + pmin pitch search range minimum + pmax pitch search range maximum + step pitch search step size + model current pitch estimate in model.Wo + + model refined pitch estimate in model.Wo + +\*---------------------------------------------------------------------------*/ + +void hs_pitch_refinement(MODEL *model, COMP Sw[], float pmin, float pmax, float pstep) +{ + int m; /* loop variable */ + int b; /* bin for current harmonic centre */ + float E; /* energy for current pitch*/ + float Wo; /* current "test" fundamental freq. */ + float Wom; /* Wo that maximises E */ + float Em; /* mamimum energy */ + float r; /* number of rads/bin */ + float p; /* current pitch */ + + /* Initialisation */ + + model->L = PI/model->Wo; /* use initial pitch est. for L */ + Wom = model->Wo; + Em = 0.0; + r = TWO_PI/FFT_ENC; + + /* Determine harmonic sum for a range of Wo values */ + + for(p=pmin; p<=pmax; p+=pstep) { + E = 0.0; + Wo = TWO_PI/p; + + /* Sum harmonic magnitudes */ + + for(m=1; m<=model->L; m++) { + b = floor(m*Wo/r + 0.5); + E += Sw[b].real*Sw[b].real + Sw[b].imag*Sw[b].imag; + } + + /* Compare to see if this is a maximum */ + + if (E > Em) { + Em = E; + Wom = Wo; + } + } + + model->Wo = Wom; +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: estimate_amplitudes + AUTHOR......: David Rowe + DATE CREATED: 27/5/94 + + Estimates the complex amplitudes of the harmonics. + +\*---------------------------------------------------------------------------*/ + +void estimate_amplitudes(MODEL *model, COMP Sw[], COMP W[]) +{ + int i,m; /* loop variables */ + int am,bm; /* bounds of current harmonic */ + int b; /* DFT bin of centre of current harmonic */ + float den; /* denominator of amplitude expression */ + float r; /* number of rads/bin */ + int offset; + COMP Am; + + r = TWO_PI/FFT_ENC; + + for(m=1; m<=model->L; m++) { + den = 0.0; + am = floor((m - 0.5)*model->Wo/r + 0.5); + bm = floor((m + 0.5)*model->Wo/r + 0.5); + b = floor(m*model->Wo/r + 0.5); + + /* Estimate ampltude of harmonic */ + + den = 0.0; + Am.real = Am.imag = 0.0; + for(i=am; i<bm; i++) { + den += Sw[i].real*Sw[i].real + Sw[i].imag*Sw[i].imag; + offset = i + FFT_ENC/2 - floor(m*model->Wo/r + 0.5); + Am.real += Sw[i].real*W[offset].real; + Am.imag += Sw[i].imag*W[offset].real; + } + + model->A[m] = sqrt(den); + + /* Estimate phase of harmonic */ + + model->phi[m] = atan2(Sw[b].imag,Sw[b].real); + } +} + +/*---------------------------------------------------------------------------*\ + + est_voicing_mbe() + + Returns the error of the MBE cost function for a fiven F0. + + Note: I think a lot of the operations below can be simplified as + W[].imag = 0 and has been normalised such that den always equals 1. + +\*---------------------------------------------------------------------------*/ + +float est_voicing_mbe( + MODEL *model, + COMP Sw[], + COMP W[], + COMP Sw_[], /* DFT of all voiced synthesised signal */ + /* useful for debugging/dump file */ + COMP Ew[], /* DFT of error */ + float prev_Wo) +{ + int i,l,al,bl,m; /* loop variables */ + COMP Am; /* amplitude sample for this band */ + int offset; /* centers Hw[] about current harmonic */ + float den; /* denominator of Am expression */ + float error; /* accumulated error between original and synthesised */ + float Wo; + float sig, snr; + float elow, ehigh, eratio; + float dF0, sixty; + + sig = 0.0; + for(l=1; l<=model->L/4; l++) { + sig += model->A[l]*model->A[l]; + } + for(i=0; i<FFT_ENC; i++) { + Sw_[i].real = 0.0; + Sw_[i].imag = 0.0; + Ew[i].real = 0.0; + Ew[i].imag = 0.0; + } + + Wo = model->Wo; + error = 0.0; + + /* Just test across the harmonics in the first 1000 Hz (L/4) */ + + for(l=1; l<=model->L/4; l++) { + Am.real = 0.0; + Am.imag = 0.0; + den = 0.0; + al = ceil((l - 0.5)*Wo*FFT_ENC/TWO_PI); + bl = ceil((l + 0.5)*Wo*FFT_ENC/TWO_PI); + + /* Estimate amplitude of harmonic assuming harmonic is totally voiced */ + + for(m=al; m<bl; m++) { + offset = FFT_ENC/2 + m - l*Wo*FFT_ENC/TWO_PI + 0.5; + Am.real += Sw[m].real*W[offset].real + Sw[m].imag*W[offset].imag; + Am.imag += Sw[m].imag*W[offset].real - Sw[m].real*W[offset].imag; + den += W[offset].real*W[offset].real + W[offset].imag*W[offset].imag; + } + + Am.real = Am.real/den; + Am.imag = Am.imag/den; + + /* Determine error between estimated harmonic and original */ + + for(m=al; m<bl; m++) { + offset = FFT_ENC/2 + m - l*Wo*FFT_ENC/TWO_PI + 0.5; + Sw_[m].real = Am.real*W[offset].real - Am.imag*W[offset].imag; + Sw_[m].imag = Am.real*W[offset].imag + Am.imag*W[offset].real; + Ew[m].real = Sw[m].real - Sw_[m].real; + Ew[m].imag = Sw[m].imag - Sw_[m].imag; + error += Ew[m].real*Ew[m].real; + error += Ew[m].imag*Ew[m].imag; + } + } + + snr = 10.0*log10(sig/error); + if (snr > V_THRESH) + model->voiced = 1; + else + model->voiced = 0; + + /* post processing, helps clean up some voicing errors ------------------*/ + + /* + Determine the ratio of low freancy to high frequency energy, + voiced speech tends to be dominated by low frequency energy, + unvoiced by high frequency. This measure can be used to + determine if we have made any gross errors. + */ + + elow = ehigh = 0.0; + for(l=1; l<=model->L/2; l++) { + elow += model->A[l]*model->A[l]; + } + for(l=model->L/2; l<=model->L; l++) { + ehigh += model->A[l]*model->A[l]; + } + eratio = 10.0*log10(elow/ehigh); + dF0 = 0.0; + + /* Look for Type 1 errors, strongly V speech that has been + accidentally declared UV */ + + if (model->voiced == 0) + if (eratio > 10.0) + model->voiced = 1; + + /* Look for Type 2 errors, strongly UV speech that has been + accidentally declared V */ + + if (model->voiced == 1) { + if (eratio < -10.0) + model->voiced = 0; + + /* If pitch is jumping about it's likely this is UV */ + + dF0 = (model->Wo - prev_Wo)*FS/TWO_PI; + if (fabs(dF0) > 15.0) + model->voiced = 0; + + /* A common source of Type 2 errors is the pitch estimator + gives a low (50Hz) estimate for UV speech, which gives a + good match with noise due to the close harmoonic spacing. + These errors are much more common than people with 50Hz + pitch, so we have just a small eratio threshold. */ + + sixty = 60.0*TWO_PI/FS; + if ((eratio < -4.0) && (model->Wo <= sixty)) + model->voiced = 0; + } + //printf(" v: %d snr: %f eratio: %3.2f %f\n",model->voiced,snr,eratio,dF0); + + return snr; +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: make_synthesis_window + AUTHOR......: David Rowe + DATE CREATED: 11/5/94 + + Init function that generates the trapezoidal (Parzen) sythesis window. + +\*---------------------------------------------------------------------------*/ + +void make_synthesis_window(float Pn[]) +{ + int i; + float win; + + /* Generate Parzen window in time domain */ + + win = 0.0; + for(i=0; i<N/2-TW; i++) + Pn[i] = 0.0; + win = 0.0; + for(i=N/2-TW; i<N/2+TW; win+=1.0/(2*TW), i++ ) + Pn[i] = win; + for(i=N/2+TW; i<3*N/2-TW; i++) + Pn[i] = 1.0; + win = 1.0; + for(i=3*N/2-TW; i<3*N/2+TW; win-=1.0/(2*TW), i++) + Pn[i] = win; + for(i=3*N/2+TW; i<2*N; i++) + Pn[i] = 0.0; +} + +/*---------------------------------------------------------------------------*\ + + FUNCTION....: synthesise + AUTHOR......: David Rowe + DATE CREATED: 20/2/95 + + Synthesise a speech signal in the frequency domain from the + sinusodal model parameters. Uses overlap-add with a trapezoidal + window to smoothly interpolate betwen frames. + +\*---------------------------------------------------------------------------*/ + +void synthesise( + float Sn_[], /* time domain synthesised signal */ + MODEL *model, /* ptr to model parameters for this frame */ + float Pn[], /* time domain Parzen window */ + int shift /* flag used to handle transition frames */ +) +{ + int i,l,j,b; /* loop variables */ + COMP Sw_[FFT_DEC]; /* DFT of synthesised signal */ + + if (shift) { + /* Update memories */ + + for(i=0; i<N-1; i++) { + Sn_[i] = Sn_[i+N]; + } + Sn_[N-1] = 0.0; + } + + for(i=0; i<FFT_DEC; i++) { + Sw_[i].real = 0.0; + Sw_[i].imag = 0.0; + } + + /* + Nov 2010 - found that synthesis using time domain cos() functions + gives better results for synthesis frames greater than 10ms. Inverse + FFT synthesis using a 512 pt FFT works well for 10ms window. I think + (but am not sure) that the problem is realted to the quantisation of + the harmonic frequencies to the FFT bin size, e.g. there is a + 8000/512 Hz step between FFT bins. For some reason this makes + the speech from longer frame > 10ms sound poor. The effect can also + be seen when synthesising test signals like single sine waves, some + sort of amplitude modulation at the frame rate. + + Another possibility is using a larger FFT size (1024 or 2048). + */ + +#define FFT_SYNTHESIS +#ifdef FFT_SYNTHESIS + /* Now set up frequency domain synthesised speech */ + for(l=1; l<=model->L; l++) { + b = floor(l*model->Wo*FFT_DEC/TWO_PI + 0.5); + if (b > ((FFT_DEC/2)-1)) { + b = (FFT_DEC/2)-1; + } + Sw_[b].real = model->A[l]*cos(model->phi[l]); + Sw_[b].imag = model->A[l]*sin(model->phi[l]); + Sw_[FFT_DEC-b].real = Sw_[b].real; + Sw_[FFT_DEC-b].imag = -Sw_[b].imag; + } + + /* Perform inverse DFT */ + + fft(&Sw_[0].real,FFT_DEC,1); +#else + /* + Direct time domain synthesis using the cos() function. Works + well at 10ms and 20ms frames rates. Note synthesis window is + still used to handle overlap-add between adjacent frames. This + could be simplified as we don't need to synthesise where Pn[] + is zero. + */ + for(l=1; l<=model->L; l++) { + for(i=0,j=-N+1; i<N-1; i++,j++) { + Sw_[FFT_DEC-N+1+i].real += 2.0*model->A[l]*cos(j*model->Wo*l + model->phi[l]); + } + for(i=N-1,j=0; i<2*N; i++,j++) + Sw_[j].real += 2.0*model->A[l]*cos(j*model->Wo*l + model->phi[l]); + } +#endif + + /* Overlap add to previous samples */ + + for(i=0; i<N-1; i++) { + Sn_[i] += Sw_[FFT_DEC-N+1+i].real*Pn[i]; + } + + if (shift) + for(i=N-1,j=0; i<2*N; i++,j++) + Sn_[i] = Sw_[j].real*Pn[i]; + else + for(i=N-1,j=0; i<2*N; i++,j++) + Sn_[i] += Sw_[j].real*Pn[i]; +} + diff --git a/gr-codec2-vocoder/src/lib/codec2/sine.h b/gr-codec2-vocoder/src/lib/codec2/sine.h new file mode 100644 index 000000000..ae578bf70 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2/sine.h @@ -0,0 +1,44 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: sine.h + AUTHOR......: David Rowe + DATE CREATED: 1/11/94 + + Header file for sinusoidal analysis and synthesis functions. + +\*---------------------------------------------------------------------------*/ + +/* + 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 Lesser General Public License version 2.1, 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 Lesser General Public License + along with this program; if not, see <http://www.gnu.org/licenses/>. +*/ + +#ifndef __SINE__ +#define __SINE__ + +#include "defines.h" +#include "comp.h" + +void make_analysis_window(float w[], COMP W[]); +float hpf(float x, float states[]); +void dft_speech(COMP Sw[], float Sn[], float w[]); +void two_stage_pitch_refinement(MODEL *model, COMP Sw[]); +void estimate_amplitudes(MODEL *model, COMP Sw[], COMP W[]); +float est_voicing_mbe(MODEL *model, COMP Sw[], COMP W[], COMP Sw_[],COMP Ew[], + float prev_Wo); +void make_synthesis_window(float Pn[]); +void synthesise(float Sn_[], MODEL *model, float Pn[], int shift); + +#endif diff --git a/gr-codec2-vocoder/src/lib/codec2_decode_ps.cc b/gr-codec2-vocoder/src/lib/codec2_decode_ps.cc new file mode 100644 index 000000000..797a9ba24 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2_decode_ps.cc @@ -0,0 +1,72 @@ +/* -*- c++ -*- */ +/* + * Copyright 2005,2010 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * GNU Radio is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3, or (at your option) + * any later version. + * + * GNU Radio 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 GNU Radio; see the file COPYING. If not, write to + * the Free Software Foundation, Inc., 51 Franklin Street, + * Boston, MA 02110-1301, USA. + */ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif +#include "codec2_decode_ps.h" +extern "C"{ +#include "codec2/codec2.h" +} +#include <gr_io_signature.h> +#include <stdexcept> +#include <assert.h> + +codec2_decode_ps_sptr +codec2_make_decode_ps () +{ + return gnuradio::get_initial_sptr(new codec2_decode_ps ()); +} + +codec2_decode_ps::codec2_decode_ps () + : gr_sync_interpolator ("codec2_decode_ps", + gr_make_io_signature (1, 1, CODEC2_BITS_PER_FRAME * sizeof (char)), + gr_make_io_signature (1, 1, sizeof (short)), + CODEC2_SAMPLES_PER_FRAME) +{ + if ((d_codec2 = codec2_create ()) == 0) + throw std::runtime_error ("codec2_decode_ps: codec2_create failed"); +} + +codec2_decode_ps::~codec2_decode_ps () +{ + codec2_destroy(d_codec2); +} + +int +codec2_decode_ps::work (int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items) +{ + const unsigned char *in = (const unsigned char *) input_items[0]; + short *out = (short *) output_items[0]; + + assert ((noutput_items % CODEC2_SAMPLES_PER_FRAME) == 0); + + for (int i = 0; i < noutput_items; i += CODEC2_SAMPLES_PER_FRAME){ + codec2_decode (d_codec2, out, const_cast<unsigned char*>(in)); + in += CODEC2_BITS_PER_FRAME * sizeof (char); + out += CODEC2_SAMPLES_PER_FRAME; + } + + return noutput_items; +} diff --git a/gr-codec2-vocoder/src/lib/codec2_decode_ps.h b/gr-codec2-vocoder/src/lib/codec2_decode_ps.h new file mode 100644 index 000000000..06c850a78 --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2_decode_ps.h @@ -0,0 +1,50 @@ +/* -*- c++ -*- */ +/* + * Copyright 2005 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * GNU Radio is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3, or (at your option) + * any later version. + * + * GNU Radio 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 GNU Radio; see the file COPYING. If not, write to + * the Free Software Foundation, Inc., 51 Franklin Street, + * Boston, MA 02110-1301, USA. + */ +#ifndef INCLUDED_CODEC2_DECODE_PS_H +#define INCLUDED_CODEC2_DECODE_PS_H + +#include <gr_sync_interpolator.h> + +class codec2_decode_ps; +typedef boost::shared_ptr<codec2_decode_ps> codec2_decode_ps_sptr; + +codec2_decode_ps_sptr codec2_make_decode_ps (); + +/*! + * \brief CODEC2 Vocoder Decoder + * \ingroup vocoder_blk + */ +class codec2_decode_ps : public gr_sync_interpolator { + void *d_codec2; + + friend codec2_decode_ps_sptr codec2_make_decode_ps (); + codec2_decode_ps (); + +public: + ~codec2_decode_ps (); + + int work (int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items); +}; + +#endif /* INCLUDED_CODEC2_DECODE_PS_H */ diff --git a/gr-codec2-vocoder/src/lib/codec2_encode_sp.cc b/gr-codec2-vocoder/src/lib/codec2_encode_sp.cc new file mode 100644 index 000000000..2f348693c --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2_encode_sp.cc @@ -0,0 +1,69 @@ +/* -*- c++ -*- */ +/* + * Copyright 2005,2010 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * GNU Radio is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3, or (at your option) + * any later version. + * + * GNU Radio 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 GNU Radio; see the file COPYING. If not, write to + * the Free Software Foundation, Inc., 51 Franklin Street, + * Boston, MA 02110-1301, USA. + */ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif +#include "codec2_encode_sp.h" +extern "C"{ +#include "codec2/codec2.h" +} +#include <gr_io_signature.h> +#include <stdexcept> + +codec2_encode_sp_sptr +codec2_make_encode_sp () +{ + return gnuradio::get_initial_sptr(new codec2_encode_sp ()); +} + +codec2_encode_sp::codec2_encode_sp () + : gr_sync_decimator ("codec2_encode_sp", + gr_make_io_signature (1, 1, sizeof (short)), + gr_make_io_signature (1, 1, CODEC2_BITS_PER_FRAME * sizeof (char)), + CODEC2_SAMPLES_PER_FRAME) +{ + if ((d_codec2 = codec2_create ()) == 0) + throw std::runtime_error ("codec2_encode_sp: codec2_create failed"); +} + +codec2_encode_sp::~codec2_encode_sp () +{ + codec2_destroy(d_codec2); +} + +int +codec2_encode_sp::work (int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items) +{ + const short *in = (const short *) input_items[0]; + unsigned char *out = (unsigned char *) output_items[0]; + + for (int i = 0; i < noutput_items; i++){ + codec2_encode (d_codec2, out, const_cast<short*>(in)); + in += CODEC2_SAMPLES_PER_FRAME; + out += CODEC2_BITS_PER_FRAME * sizeof (char); + } + + return noutput_items; +} diff --git a/gr-codec2-vocoder/src/lib/codec2_encode_sp.h b/gr-codec2-vocoder/src/lib/codec2_encode_sp.h new file mode 100644 index 000000000..bec90181d --- /dev/null +++ b/gr-codec2-vocoder/src/lib/codec2_encode_sp.h @@ -0,0 +1,50 @@ +/* -*- c++ -*- */ +/* + * Copyright 2005 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * GNU Radio is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 3, or (at your option) + * any later version. + * + * GNU Radio 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 GNU Radio; see the file COPYING. If not, write to + * the Free Software Foundation, Inc., 51 Franklin Street, + * Boston, MA 02110-1301, USA. + */ +#ifndef INCLUDED_CODEC2_ENCODE_SP_H +#define INCLUDED_CODEC2_ENCODE_SP_H + +#include <gr_sync_decimator.h> + +class codec2_encode_sp; +typedef boost::shared_ptr<codec2_encode_sp> codec2_encode_sp_sptr; + +codec2_encode_sp_sptr codec2_make_encode_sp (); + +/*! + * \brief CODEC2 Vocoder Encoder + * \ingroup vocoder_blk + */ +class codec2_encode_sp : public gr_sync_decimator { + void *d_codec2; + + friend codec2_encode_sp_sptr codec2_make_encode_sp (); + codec2_encode_sp (); + +public: + ~codec2_encode_sp (); + + int work (int noutput_items, + gr_vector_const_void_star &input_items, + gr_vector_void_star &output_items); +}; + +#endif /* INCLUDED_CODEC2_ENCODE_SP_H */ |