From: drowe67 Date: Tue, 18 Sep 2012 04:46:29 +0000 (+0000) Subject: modified to correct for MAX_AMP+1 overflow error, added mbest VQ search X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=5b090e7b8caec286053a4c7fe2b39c05d4b49b0f;p=freetel-svn-tracking.git modified to correct for MAX_AMP+1 overflow error, added mbest VQ search git-svn-id: https://svn.code.sf.net/p/freetel/code@723 01035d8c-6547-0410-b346-abe4f91aad63 --- diff --git a/codec2-dev/README b/codec2-dev/README index 4a87406a..c8398883 100644 --- a/codec2-dev/README +++ b/codec2-dev/README @@ -64,19 +64,21 @@ Debugging Directories ----------- - fltk - FLTK GUI programs(s) - octave - Octave scripts used for visualising internal signals - during development - script - shell scripts for playing and converting raw files - src - C source code - raw - speech files in raw format (16 bits signed linear 8 kHz) - unittest - unit test source code - voicing - hand-estimated voicing files, used for development - wav - speech files in wave file format - win32 - Support for building Windows DLL version of Codec 2 and FDMDV libraries + fltk - FLTK GUI programs(s) + octave - Octave scripts used for visualising internal signals + during development + portaudio - Portaudio test programs + script - shell scripts for playing and converting raw files + src - C source code + raw - speech files in raw format (16 bits signed linear 8 kHz) + unittest - unit test source code + voicing - hand-estimated voicing files, used for development + wav - speech files in wave file format + win32 - Support for building Windows DLL version of Codec 2 and FDMDV libraries TODO ---- -[ ] Get win32/Makefile integarted into Automake system, such that if +[ ] Get win32/Makefile integrated into Automake system, such that if i586-mingw32msvc exists the Win32 code gets automatically built. +[ ] Same for fltk & portaudio, build these conditionally if libs exist \ No newline at end of file diff --git a/codec2-dev/fltk/fl_fdmdv.cxx b/codec2-dev/fltk/fl_fdmdv.cxx index ae8232c0..30c67aab 100644 --- a/codec2-dev/fltk/fl_fdmdv.cxx +++ b/codec2-dev/fltk/fl_fdmdv.cxx @@ -638,6 +638,15 @@ void new_data(float mag_dB[]) { The ouput of the demod is codec voice data so it's OK if we miss or repeat a frame every now and again. + TODOs: + + + this might work with arbitrary input and output buffer lengths, + 0,1, or 2 only apply if we are inputting the nominal number of + samples on every call. + + + so the I/O buffer sizes might not matter, as long as they of + reasonable size (say twice the nominal size). + \*------------------------------------------------------------------*/ void per_frame_rx_processing(short output_buf[], /* output buf of decoded speech samples */ diff --git a/codec2-dev/octave/lspwarp.m b/codec2-dev/octave/lspwarp.m new file mode 100644 index 00000000..2d2f2c99 --- /dev/null +++ b/codec2-dev/octave/lspwarp.m @@ -0,0 +1,40 @@ +% lspwarp.m +% David Rowe Sep 2012 +% +% Experimenting with non-linear LSP frequency axis for LSP quantisation +% Plots a scaled mel axis. + +1; + +function mel = freq2mel(f) + mel = 70*log10(1 + f/700); +endfunction + +function freq = mel2freq(m) + freq = 700*(10 ^ (m/70) - 1); +endfunction + +x = []; y = []; + +for freq = 100:25:4000 + mel = freq2mel(freq); + x = [x freq]; + y = [y mel]; +end + +plot(x,y) +grid + +mel_start = floor(freq2mel(100)); +mel_end = floor(freq2mel(4000)); + +x = []; y = []; +for mel=mel_start:mel_end + freq = mel2freq(mel); + x = [x freq]; + y = [y mel]; +end + +hold on; +plot(x,y, '+') +hold off; diff --git a/codec2-dev/octave/pllpcpf.m b/codec2-dev/octave/pllpcpf.m index b3912f8a..924e045a 100644 --- a/codec2-dev/octave/pllpcpf.m +++ b/codec2-dev/octave/pllpcpf.m @@ -9,10 +9,13 @@ function pllpcpf(samname, f) % switch some stuff off to unclutter display + plot_Am = 0; + plot_Amq = 0; + plot_err = 0; plot_lsp = 0; plot_snr = 0; plot_vsnr = 0; - plot_sw = 1; + plot_sw = 0; plot_pw = 1; plot_pwb = 1; plot_rw = 1; @@ -69,10 +72,13 @@ function pllpcpf(samname, f) axis([1 length(s) -20000 20000]); figure(2); + clf; Wo = model(f,1); L = model(f,2); Am = model(f,3:(L+2)); - plot((1:L)*Wo*4000/pi, 20*log10(Am),";Am;r"); + if plot_Am + plot((1:L)*Wo*4000/pi, 20*log10(Am),";Am;r"); + end axis([1 4000 -10 80]); hold on; if plot_sw @@ -82,25 +88,29 @@ function pllpcpf(samname, f) if (file_in_path(".",modelq_name)) Amq = modelq(f,3:(L+2)); - plot((1:L)*Wo*4000/pi, 20*log10(Amq),";Amq;g" ); + if plot_Amq + plot((1:L)*Wo*4000/pi, 20*log10(Amq),";Amq;g" ); + end if (file_in_path(".",pwb_name) && plot_pwb) - plot((0:255)*4000/256, 10*log10(Pwb(f,:)),";Pwb;c"); + plot((0:255)*4000/256, 10*log10(Pwb(f,:)),";Pwb;r"); endif if (file_in_path(".",rw_name) && plot_rw) - plot((0:255)*4000/256, 10*log10(Rw(f,:)),";Rw;c"); + plot((0:255)*4000/256, 10*log10(Rw(f,:)),";Rw;b"); endif if (file_in_path(".",pw_name) && plot_pw) - plot((0:255)*4000/256, 10*log10(Pw(f,:)),";Pw;c."); + plot((0:255)*4000/256, 10*log10(Pw(f,:)),";Pw;g."); endif signal = Am * Am'; noise = (Am-Amq) * (Am-Amq)'; snr1 = 10*log10(signal/noise); Am_err_label = sprintf(";Am error SNR %4.2f dB;m",snr1); - plot((1:L)*Wo*4000/pi, 20*log10(Amq) - 20*log10(Am), Am_err_label); + if plot_err + plot((1:L)*Wo*4000/pi, 20*log10(Amq) - 20*log10(Am), Am_err_label); + end endif diff --git a/codec2-dev/octave/sd.m b/codec2-dev/octave/sd.m index 112af614..4c3d13c9 100644 --- a/codec2-dev/octave/sd.m +++ b/codec2-dev/octave/sd.m @@ -26,6 +26,7 @@ function sd(raw_filename, dump_file_prefix, f) A2(i,:) = -20*log10(abs(fft(ak2(i,:),Ndft))); sd(i) = sum((A1(i,:) - A2(i,:)).^2)/Ndft; end + printf("sd av %3.2f dB*dB\n", sum(sd)/frames); figure(1); clf; @@ -46,6 +47,11 @@ function sd(raw_filename, dump_file_prefix, f) lsp1 = load(lsp1_filename); lsp2 = load(lsp2_filename); + weights_filename = sprintf("%s_weights.txt", dump_file_prefix); + if file_in_path(".",weights_filename) + weights = load(weights_filename); + end + k = ' '; do figure(2); @@ -61,7 +67,10 @@ function sd(raw_filename, dump_file_prefix, f) axis([1 4000 -20 40]); hold on; plot((1:Ndft/2)*4000/(Ndft/2), A2(f,1:(Ndft/2)),";A2;"); - + if file_in_path(".",weights_filename) + plot(lsp1(f,:)*4000/pi, weights(f,:),";weights;g+"); + end + for l=1:10 plot([lsp1(f,l)*4000/pi lsp1(f,l)*4000/pi], [0 -10], 'r'); plot([lsp2(f,l)*4000/pi lsp2(f,l)*4000/pi], [-10 -20], 'b'); diff --git a/codec2-dev/portaudio/Makefile b/codec2-dev/portaudio/Makefile index 6c1e2c12..a1876de9 100644 --- a/codec2-dev/portaudio/Makefile +++ b/codec2-dev/portaudio/Makefile @@ -5,7 +5,7 @@ CFLAGS = -g -Wall -I../src LIBS = -lm -lportaudio -pthread SRC = ../src/fdmdv.c ../src/kiss_fft.c -all: pa_rec pa_play pa_recplay +all: pa_rec pa_play pa_recplay pa_impresp pa_rec: Makefile pa_rec.c $(SRC) gcc $(CFLAGS) pa_rec.c $(SRC) -o pa_rec $(LIBS) @@ -16,5 +16,8 @@ pa_play: Makefile pa_play.c $(SRC) pa_recplay: Makefile pa_recplay.c $(SRC) gcc $(CFLAGS) pa_recplay.c $(SRC) -o pa_recplay $(LIBS) +pa_impresp: Makefile pa_impresp.c $(SRC) + gcc $(CFLAGS) pa_impresp.c $(SRC) -o pa_impresp $(LIBS) + clean: rm -f pa_rec pa_play pa_recplay diff --git a/codec2-dev/portaudio/pa_impresp.c b/codec2-dev/portaudio/pa_impresp.c new file mode 100644 index 00000000..7317e07e --- /dev/null +++ b/codec2-dev/portaudio/pa_impresp.c @@ -0,0 +1,254 @@ +/* + pa_impresp.c + David Rowe + August 29 2012 + + Measures the impulse reponse of the path between the speaker and + microphone. Used to explore why Codec audio quality is + different through a speaker and headphones. + + Modified from pa_playrec.c +*/ + +/* + * $Id: paex_record.c 1752 2011-09-08 03:21:55Z philburk $ + * + * This program uses the PortAudio Portable Audio Library. + * For more information see: http://www.portaudio.com + * Copyright (c) 1999-2000 Ross Bencina and Phil Burk + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files + * (the "Software"), to deal in the Software without restriction, + * including without limitation the rights to use, copy, modify, merge, + * publish, distribute, sublicense, and/or sell copies of the Software, + * and to permit persons to whom the Software is furnished to do so, + * subject to the following conditions: + * + * The above copyright notice and this permission notice shall be + * included in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. + * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR + * ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF + * CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#include +#include +#include +#include "portaudio.h" +#include "fdmdv.h" + +#define SAMPLE_RATE 48000 /* 48 kHz sampling rate rec. as we + can trust accuracy of sound + card */ +#define N8 160 /* processing buffer size at 8 kHz */ +#define N48 (N8*FDMDV_OS) /* processing buffer size at 48 kHz */ +#define MEM8 (FDMDV_OS_TAPS/FDMDV_OS) +#define NUM_CHANNELS 2 /* I think most sound cards prefer + stereo, we will convert to mono + as we sample */ + +#define IMPULSE_AMP 16384 /* amplitide of impulse */ +#define IMPULSE_PERIOD 0.1 /* period (dly between impulses) in secs */ + +/* state information passed to call back */ + +typedef struct { + float in48k[FDMDV_OS_TAPS + N48]; + float in8k[MEM8 + N8]; + FILE *fimp; + float *impulse_buf; + int impulse_buf_length; + int impulse_sample_count; + int framesLeft; +} paTestData; + + +/* + This routine will be called by the PortAudio engine when audio is + required. It may be called at interrupt level on some machines so + don't do anything that could mess up the system like calling + malloc() or free(). +*/ + +static int callback( const void *inputBuffer, void *outputBuffer, + unsigned long framesPerBuffer, + const PaStreamCallbackTimeInfo* timeInfo, + PaStreamCallbackFlags statusFlags, + void *userData ) +{ + paTestData *data = (paTestData*)userData; + int i; + short *rptr = (short*)inputBuffer; + short *wptr = (short*)outputBuffer; + float *in8k = data->in8k; + float *in48k = data->in48k; + float out8k[N8]; + float out48k[N48]; + short out48k_short[N48]; + short out8k_short[N8]; + + (void) timeInfo; + (void) statusFlags; + + assert(inputBuffer != NULL); + + /* just use left channel */ + + for(i=0; iin48k[i+FDMDV_OS_TAPS] = *rptr; + + /* downsample and update filter memory */ + + fdmdv_48_to_8(out8k, &in48k[FDMDV_OS_TAPS], N8); + for(i=0; ifimp); + + /* play side, read from impulse buffer */ + + for(i=0; iimpulse_buf[data->impulse_sample_count]; + data->impulse_sample_count++; + if (data->impulse_sample_count == data->impulse_buf_length) + data->impulse_sample_count = 0; + } + + /* upsample and update filter memory */ + + fdmdv_8_to_48(out48k, &in8k[MEM8], N8); + for(i=0; iframesLeft -= framesPerBuffer; + if (data->framesLeft > 0) + return paContinue; + else + return paComplete; +} + +int main(int argc, char *argv[]) +{ + PaStreamParameters inputParameters, outputParameters; + PaStream* stream; + PaError err = paNoError; + paTestData data; + int i, numSecs; + + if (argc != 3) { + printf("usage: %s impulseRawFile time(s)\n", argv[0]); + exit(0); + } + + data.fimp = fopen(argv[1], "wb"); + if (data.fimp == NULL) { + printf("Error opening impulse output file %s\n", argv[1]); + exit(1); + } + + numSecs = atoi(argv[2]); + data.framesLeft = numSecs * SAMPLE_RATE; + + /* init filter states */ + + for(i=0; idefaultLowInputLatency; + inputParameters.hostApiSpecificStreamInfo = NULL; + + outputParameters.device = Pa_GetDefaultOutputDevice(); /* default output device */ + if (outputParameters.device == paNoDevice) { + fprintf(stderr,"Error: No default output device.\n"); + goto done; + } + outputParameters.channelCount = NUM_CHANNELS; /* stereo output */ + outputParameters.sampleFormat = paInt16; + outputParameters.suggestedLatency = Pa_GetDeviceInfo( outputParameters.device )->defaultLowOutputLatency; + outputParameters.hostApiSpecificStreamInfo = NULL; + + /* Play some audio --------------------------------------------- */ + + err = Pa_OpenStream( + &stream, + &inputParameters, + &outputParameters, + SAMPLE_RATE, + N48, + paClipOff, + callback, + &data ); + if( err != paNoError ) goto done; + + err = Pa_StartStream( stream ); + if( err != paNoError ) goto done; + + while( ( err = Pa_IsStreamActive( stream ) ) == 1 ) + { + Pa_Sleep(100); + } + if( err < 0 ) goto done; + + err = Pa_CloseStream( stream ); + if( err != paNoError ) goto done; + + +done: + Pa_Terminate(); + if( err != paNoError ) + { + fprintf( stderr, "An error occured while using the portaudio stream\n" ); + fprintf( stderr, "Error number: %d\n", err ); + fprintf( stderr, "Error message: %s\n", Pa_GetErrorText( err ) ); + err = 1; /* Always return 0 or 1, but no other return codes. */ + } + + fclose(data.fimp); + + return err; +} + diff --git a/codec2-dev/src/Makefile.am b/codec2-dev/src/Makefile.am index 097f4cce..3d2e0c89 100644 --- a/codec2-dev/src/Makefile.am +++ b/codec2-dev/src/Makefile.am @@ -127,7 +127,6 @@ interp.c \ lsp.c \ phase.c \ quantise.c \ -ampexp.c \ pack.c \ codebook.c \ codebookd.c \ @@ -159,7 +158,7 @@ c2dec_SOURCES = c2dec.c c2dec_LDADD = $(lib_LTLIBRARIES) c2dec_LDFLAGS = $(LIBS) -c2sim_SOURCES = c2sim.c +c2sim_SOURCES = c2sim.c ampexp.c phaseexp.c c2sim_LDADD = $(lib_LTLIBRARIES) c2sim_LDFLAGS = $(LIBS) diff --git a/codec2-dev/src/Makefile.in b/codec2-dev/src/Makefile.in index f186e93b..73dcd428 100644 --- a/codec2-dev/src/Makefile.in +++ b/codec2-dev/src/Makefile.in @@ -69,11 +69,11 @@ am_libcodec2_la_OBJECTS = libcodec2_la-dump.lo libcodec2_la-lpc.lo \ libcodec2_la-fdmdv.lo libcodec2_la-kiss_fft.lo \ libcodec2_la-interp.lo libcodec2_la-lsp.lo \ libcodec2_la-phase.lo libcodec2_la-quantise.lo \ - libcodec2_la-ampexp.lo libcodec2_la-pack.lo \ - libcodec2_la-codebook.lo libcodec2_la-codebookd.lo \ - libcodec2_la-codebookvq.lo libcodec2_la-codebookjnd.lo \ - libcodec2_la-codebookjvm.lo libcodec2_la-codebookvqanssi.lo \ - libcodec2_la-codebookdt.lo libcodec2_la-codebookge.lo + libcodec2_la-pack.lo libcodec2_la-codebook.lo \ + libcodec2_la-codebookd.lo libcodec2_la-codebookvq.lo \ + libcodec2_la-codebookjnd.lo libcodec2_la-codebookjvm.lo \ + libcodec2_la-codebookvqanssi.lo libcodec2_la-codebookdt.lo \ + libcodec2_la-codebookge.lo libcodec2_la_OBJECTS = $(am_libcodec2_la_OBJECTS) binPROGRAMS_INSTALL = $(INSTALL_PROGRAM) PROGRAMS = $(bin_PROGRAMS) $(noinst_PROGRAMS) @@ -87,7 +87,7 @@ c2demo_DEPENDENCIES = $(am__DEPENDENCIES_1) am_c2enc_OBJECTS = c2enc.$(OBJEXT) c2enc_OBJECTS = $(am_c2enc_OBJECTS) c2enc_DEPENDENCIES = $(am__DEPENDENCIES_1) -am_c2sim_OBJECTS = c2sim.$(OBJEXT) +am_c2sim_OBJECTS = c2sim.$(OBJEXT) ampexp.$(OBJEXT) phaseexp.$(OBJEXT) c2sim_OBJECTS = $(am_c2sim_OBJECTS) c2sim_DEPENDENCIES = $(am__DEPENDENCIES_1) am_fdmdv_demod_OBJECTS = fdmdv_demod.$(OBJEXT) fdmdv.$(OBJEXT) \ @@ -333,7 +333,6 @@ interp.c \ lsp.c \ phase.c \ quantise.c \ -ampexp.c \ pack.c \ codebook.c \ codebookd.c \ @@ -357,7 +356,7 @@ c2enc_LDFLAGS = $(LIBS) c2dec_SOURCES = c2dec.c c2dec_LDADD = $(lib_LTLIBRARIES) c2dec_LDFLAGS = $(LIBS) -c2sim_SOURCES = c2sim.c +c2sim_SOURCES = c2sim.c ampexp.c phaseexp.c c2sim_LDADD = $(lib_LTLIBRARIES) c2sim_LDFLAGS = $(LIBS) fdmdv_get_test_bits_SOURCES = fdmdv_get_test_bits.c fdmdv.c kiss_fft.c @@ -507,6 +506,7 @@ mostlyclean-compile: distclean-compile: -rm -f *.tab.c +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/ampexp.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/c2dec.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/c2demo.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/c2enc.Po@am__quote@ @@ -520,7 +520,6 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/generate_codebook.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/genlspdtcb.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/kiss_fft.Po@am__quote@ -@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libcodec2_la-ampexp.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libcodec2_la-codebook.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libcodec2_la-codebookd.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libcodec2_la-codebookdt.Plo@am__quote@ @@ -543,6 +542,7 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libcodec2_la-quantise.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libcodec2_la-sine.Plo@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/octave.Po@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/phaseexp.Po@am__quote@ .c.o: @am__fastdepCC_TRUE@ if $(COMPILE) -MT $@ -MD -MP -MF "$(DEPDIR)/$*.Tpo" -c -o $@ $<; \ @@ -649,13 +649,6 @@ libcodec2_la-quantise.lo: quantise.c @AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@ @am__fastdepCC_FALSE@ $(LIBTOOL) --tag=CC --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libcodec2_la_CFLAGS) $(CFLAGS) -c -o libcodec2_la-quantise.lo `test -f 'quantise.c' || echo '$(srcdir)/'`quantise.c -libcodec2_la-ampexp.lo: ampexp.c -@am__fastdepCC_TRUE@ if $(LIBTOOL) --tag=CC --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libcodec2_la_CFLAGS) $(CFLAGS) -MT libcodec2_la-ampexp.lo -MD -MP -MF "$(DEPDIR)/libcodec2_la-ampexp.Tpo" -c -o libcodec2_la-ampexp.lo `test -f 'ampexp.c' || echo '$(srcdir)/'`ampexp.c; \ -@am__fastdepCC_TRUE@ then mv -f "$(DEPDIR)/libcodec2_la-ampexp.Tpo" "$(DEPDIR)/libcodec2_la-ampexp.Plo"; else rm -f "$(DEPDIR)/libcodec2_la-ampexp.Tpo"; exit 1; fi -@AMDEP_TRUE@@am__fastdepCC_FALSE@ source='ampexp.c' object='libcodec2_la-ampexp.lo' libtool=yes @AMDEPBACKSLASH@ -@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@ -@am__fastdepCC_FALSE@ $(LIBTOOL) --tag=CC --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libcodec2_la_CFLAGS) $(CFLAGS) -c -o libcodec2_la-ampexp.lo `test -f 'ampexp.c' || echo '$(srcdir)/'`ampexp.c - libcodec2_la-pack.lo: pack.c @am__fastdepCC_TRUE@ if $(LIBTOOL) --tag=CC --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libcodec2_la_CFLAGS) $(CFLAGS) -MT libcodec2_la-pack.lo -MD -MP -MF "$(DEPDIR)/libcodec2_la-pack.Tpo" -c -o libcodec2_la-pack.lo `test -f 'pack.c' || echo '$(srcdir)/'`pack.c; \ @am__fastdepCC_TRUE@ then mv -f "$(DEPDIR)/libcodec2_la-pack.Tpo" "$(DEPDIR)/libcodec2_la-pack.Plo"; else rm -f "$(DEPDIR)/libcodec2_la-pack.Tpo"; exit 1; fi diff --git a/codec2-dev/src/c2sim.c b/codec2-dev/src/c2sim.c index 721631da..0d9b41c4 100644 --- a/codec2-dev/src/c2sim.c +++ b/codec2-dev/src/c2sim.c @@ -47,6 +47,7 @@ #include "postfilter.h" #include "interp.h" #include "ampexp.h" +#include "phaseexp.h" void synth_one_frame(kiss_fft_cfg fft_inv_cfg, short buf[], MODEL *model, float Sn_[], float Pn[]); void print_help(const struct option *long_options, int num_opts, char* argv[]); @@ -436,7 +437,7 @@ int main(int argc, char *argv[]) /* just to make sure we are not cheating - kill all phases */ - for(i=0; ibg_est = 0.0; c2->ex_phase = 0.0; - for(l=1; lprev_model_dec.A[l] = 0.0; c2->prev_model_dec.Wo = TWO_PI/P_MAX; c2->prev_model_dec.L = PI/c2->prev_model_dec.Wo; diff --git a/codec2-dev/src/defines.h b/codec2-dev/src/defines.h index ed6b435d..4870770c 100644 --- a/codec2-dev/src/defines.h +++ b/codec2-dev/src/defines.h @@ -68,8 +68,8 @@ 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 */ + float A[MAX_AMP+1]; /* amplitiude of each harmonic */ + float phi[MAX_AMP+1]; /* phase of each harmonic */ int voiced; /* non-zero if this frame is voiced */ } MODEL; diff --git a/codec2-dev/src/phase.c b/codec2-dev/src/phase.c index 27229886..d41303b7 100644 --- a/codec2-dev/src/phase.c +++ b/codec2-dev/src/phase.c @@ -200,9 +200,9 @@ void phase_synth_zero_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 */ + COMP Ex[MAX_AMP+1]; /* excitation samples */ + COMP A_[MAX_AMP+1]; /* synthesised harmonic samples */ + COMP H[MAX_AMP+1]; /* LPC freq domain samples */ float G; float jitter = 0.0; float r; @@ -271,1539 +271,3 @@ void phase_synth_zero_order( } -/* Bruce Perens' funcs to load codebook files */ - -struct codebook { - unsigned int k; - unsigned int log2m; - unsigned int m; - COMP *cb; - unsigned int offset; -}; - -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"; - -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(const char * name) -{ - FILE *file; - char line[2048]; - char *cursor = line; - struct codebook *b = malloc(sizeof(struct codebook)); - int i; - int size; - float angle; - - file = fopen(name, "rt"); - assert(file != NULL); - - *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 = (COMP *)malloc(size * sizeof(COMP)); - - for ( i = 0; i < size; i++ ) { - angle = get_float(file, name, &cursor, line, sizeof(line)); - b->cb[i].real = cos(angle); - b->cb[i].imag = sin(angle); - } - - fclose(file); - - return b; -} - - -/* states for phase experiments */ - -struct PEXP { - float phi1; - float phi_prev[MAX_AMP]; - float Wo_prev; - int frames; - float snr; - float var; - int var_n; - struct codebook *vq1,*vq2,*vq3,*vq4,*vq5; - float vq_var; - int vq_var_n; - MODEL prev_model; - int state; -}; - - -/*---------------------------------------------------------------------------* \ - - phase_experiment_create() - - Inits states for phase quantisation experiments. - -\*---------------------------------------------------------------------------*/ - -struct PEXP * phase_experiment_create() { - struct PEXP *pexp; - int i; - - pexp = (struct PEXP *)malloc(sizeof(struct PEXP)); - assert (pexp != NULL); - - pexp->phi1 = 0; - for(i=0; iphi_prev[i] = 0.0; - pexp->Wo_prev = 0.0; - pexp->frames = 0; - pexp->snr = 0.0; - pexp->var = 0.0; - pexp->var_n = 0; - - /* smoothed 10th order for 1st 1 khz */ - //pexp->vq1 = load("../unittest/ph1_10_1024.txt"); - //pexp->vq1->offset = 0; - - /* load experimental phase VQ */ - - //pexp->vq1 = load("../unittest/testn1_20_1024.txt"); - pexp->vq1 = load("../unittest/test.txt"); - //pexp->vq2 = load("../unittest/testn21_40_1024.txt"); - pexp->vq2 = load("../unittest/test11_20_1024.txt"); - pexp->vq3 = load("../unittest/test21_30_1024.txt"); - pexp->vq4 = load("../unittest/test31_40_1024.txt"); - pexp->vq5 = load("../unittest/test41_60_1024.txt"); - pexp->vq1->offset = 0; - pexp->vq2->offset = 10; - pexp->vq3->offset = 20; - pexp->vq4->offset = 30; - pexp->vq5->offset = 40; - - pexp->vq_var = 0.0; - pexp->vq_var_n = 0; - - pexp->state = 0; - - return pexp; -} - - -/*---------------------------------------------------------------------------* \ - - phase_experiment_destroy() - -\*---------------------------------------------------------------------------*/ - -void phase_experiment_destroy(struct PEXP *pexp) { - assert(pexp != NULL); - if (pexp->snr != 0.0) - printf("snr: %4.2f dB\n", pexp->snr/pexp->frames); - if (pexp->var != 0.0) - printf("var...: %4.3f std dev...: %4.3f (%d non zero phases)\n", - pexp->var/pexp->var_n, sqrt(pexp->var/pexp->var_n), pexp->var_n); - if (pexp->vq_var != 0.0) - printf("vq var: %4.3f vq std dev: %4.3f (%d non zero phases)\n", - pexp->vq_var/pexp->vq_var_n, sqrt(pexp->vq_var/pexp->vq_var_n), pexp->vq_var_n); - free(pexp); -} - - -/*---------------------------------------------------------------------------* \ - - Various test and experimental functions ................ - -\*---------------------------------------------------------------------------*/ - -/* Bubblesort to find highest amplitude harmonics */ - -struct AMPINDEX { - float amp; - int index; -}; - -static void bubbleSort(struct AMPINDEX numbers[], int array_size) -{ - int i, j; - struct AMPINDEX temp; - - for (i = (array_size - 1); i > 0; i--) - { - for (j = 1; j <= i; j++) - { - //printf("i %d j %d %f %f \n", i, j, numbers[j-1].amp, numbers[j].amp); - if (numbers[j-1].amp < numbers[j].amp) - { - temp = numbers[j-1]; - numbers[j-1] = numbers[j]; - numbers[j] = temp; - } - } - } -} - - -static void print_pred_error(struct PEXP *pexp, MODEL *model, int start, int end, float mag_thresh) { - int i; - float mag; - - mag = 0.0; - for(i=start; i<=end; i++) - mag += model->A[i]*model->A[i]; - mag = 10*log10(mag/(end-start)); - - if (mag > mag_thresh) { - for(i=start; i<=end; i++) { - float pred = pexp->phi_prev[i] + N*i*(model->Wo + pexp->Wo_prev)/2.0; - float err = pred - model->phi[i]; - err = atan2(sin(err),cos(err)); - printf("%f\n",err); - } - //printf("\n"); - } - -} - - -static void predict_phases(struct PEXP *pexp, MODEL *model, int start, int end) { - int i; - - for(i=start; i<=end; i++) { - model->phi[i] = pexp->phi_prev[i] + N*i*model->Wo; - } - -} -static float refine_Wo(struct PEXP *pexp, - MODEL *model, - int start, - int end); - -/* Fancy state based phase prediction. Actually works OK on most utterances, - but could use some tuning. Breaks down a bit on mmt1. */ - -static void predict_phases_state(struct PEXP *pexp, MODEL *model, int start, int end) { - int i, next_state; - float best_Wo, dWo; - - //best_Wo = refine_Wo(pexp, model, start, end); - //best_Wo = (model->Wo + pexp->Wo_prev)/2.0; - best_Wo = model->Wo; - - dWo = fabs(model->Wo - pexp->Wo_prev)/model->Wo; - next_state = pexp->state; - switch(pexp->state) { - case 0: - if (dWo < 0.1) { - /* UV -> V transition, so start with phases in lock. They will - drift a bit over voiced track which is kinda what we want, so - we don't get clicky speech. - */ - next_state = 1; - for(i=start; i<=end; i++) - pexp->phi_prev[i] = i*pexp->phi1; - } - - break; - case 1: - if (dWo > 0.1) - next_state = 0; - break; - } - pexp->state = next_state; - - if (pexp->state == 0) - for(i=start; i<=end; i++) { - model->phi[i] = PI*(1.0 - 2.0*rand()/RAND_MAX); - } - else - for(i=start; i<=end; i++) { - model->phi[i] = pexp->phi_prev[i] + N*i*best_Wo; - } - printf("state %d\n", pexp->state); -} - -static void struct_phases(struct PEXP *pexp, MODEL *model, int start, int end) { - int i; - - for(i=start; i<=end; i++) - model->phi[i] = pexp->phi1*i; - -} - - -static void predict_phases2(struct PEXP *pexp, MODEL *model, int start, int end) { - int i; - float pred, str, diff; - - for(i=start; i<=end; i++) { - pred = pexp->phi_prev[i] + N*i*model->Wo; - str = pexp->phi1*i; - diff = str - pred; - diff = atan2(sin(diff), cos(diff)); - if (diff > 0) - pred += PI/16; - else - pred -= PI/16; - model->phi[i] = pred; - } - -} - -static void skip_phases(struct PEXP *pexp, MODEL *model, int start, int end) { - int i; - - for(i=start; i<=end; i+=2) - model->phi[i] = model->phi[i-1] - model->phi[i-2]; - -} - -static void rand_phases(MODEL *model, int start, int end) { - int i; - - for(i=start; i<=end; i++) - model->phi[i] = PI*(1.0 - 2.0*(float)rand()/RAND_MAX); - -} - -static void quant_phase(float *phase, float min, float max, int bits) { - int levels = 1 << bits; - int index; - float norm, step; - - norm = (*phase - min)/(max - min); - index = floor(levels*norm); - - //printf("phase %f norm %f index %d ", *phase, norm, index); - if (index < 0 ) index = 0; - if (index > (levels-1)) index = levels-1; - //printf("index %d ", index); - step = (max - min)/levels; - *phase = min + step*index + 0.5*step; - //printf("step %f phase %f\n", step, *phase); -} - -static void quant_phases(MODEL *model, int start, int end, int bits) { - int i; - - for(i=start; i<=end; i++) { - quant_phase(&model->phi[i], -PI, PI, bits); - } -} - -static void fixed_bits_per_frame(struct PEXP *pexp, MODEL *model, int m, int budget) { - int res, finished; - - res = 3; - finished = 0; - - while(!finished) { - if (m > model->L/2) - res = 2; - if (((budget - res) < 0) || (m > model->L)) - finished = 1; - else { - quant_phase(&model->phi[m], -PI, PI, res); - budget -= res; - m++; - } - } - printf("m: %d L: %d budget: %d\n", m, model->L, budget); - predict_phases(pexp, model, m, model->L); - //rand_phases(model, m, model->L); -} - -/* used to plot histogram of quantisation error, for 3 bits, 8 levels, - should be uniform between +/- PI/8 */ - -static void check_phase_quant(MODEL *model, float tol) -{ - int m; - float phi_before[MAX_AMP]; - - for(m=1; m<=model->L; m++) - phi_before[m] = model->phi[m]; - - quant_phases(model, 1, model->L, 3); - - for(m=1; m<=model->L; m++) { - float err = phi_before[m] - model->phi[m]; - printf("%f\n", err); - if (fabs(err) > tol) - exit(0); - } -} - - -static float est_phi1(MODEL *model, int start, int end) -{ - int m; - float delta, s, c, phi1_est; - - if (end > model->L) - end = model->L; - - s = c = 0.0; - for(m=start; mphi[m+1] - model->phi[m]; - s += sin(delta); - c += cos(delta); - } - - phi1_est = atan2(s,c); - - return phi1_est; -} - -static void print_phi1_pred_error(MODEL *model, int start, int end) -{ - int m; - float phi1_est; - - phi1_est = est_phi1(model, start, end); - - for(m=start; mphi[m+1] - model->phi[m] - phi1_est; - err = atan2(sin(err),cos(err)); - printf("%f\n", err); - } -} - - -static void first_order_band(MODEL *model, int start, int end, float phi1_est) -{ - int m; - float pred_err, av_pred_err; - float c,s; - - s = c = 0.0; - for(m=start; mphi[m] - phi1_est*m; - s += sin(pred_err); - c += cos(pred_err); - } - - av_pred_err = atan2(s,c); - for(m=start; mphi[m] = av_pred_err + phi1_est*m; - model->phi[m] = atan2(sin(model->phi[m]), cos(model->phi[m])); - } - -} - - -static void sub_linear(MODEL *model, int start, int end, float phi1_est) -{ - int m; - - for(m=start; mphi[m] = m*phi1_est; - } -} - - -static void top_amp(struct PEXP *pexp, MODEL *model, int start, int end, int n_harm, int pred) -{ - int removed = 0, not_removed = 0; - int top, i, j; - struct AMPINDEX sorted[MAX_AMP]; - - /* sort into ascending order of amplitude */ - - printf("\n"); - for(i=start,j=0; iA[i]; - sorted[j].index = i; - printf("%f ", model->A[i]); - } - bubbleSort(sorted, end-start); - - printf("\n"); - for(j=0; jA[i] == sorted[j].amp) { - top = 1; - assert(i == sorted[j].index); - } - } - - #define ALTTOP - #ifdef ALTTOP - model->phi[i] = 0.0; /* make sure */ - if (top) { - model->phi[i] = i*pexp->phi1; - removed++; - } - else { - model->phi[i] = PI*(1.0 - 2.0*(float)rand()/RAND_MAX); // note: try rand for higher harms - removed++; - } - #else - if (!top) { - model->phi[i] = 0.0; /* make sure */ - if (pred) { - //model->phi[i] = pexp->phi_prev[i] + i*N*(model->Wo + pexp->Wo_prev)/2.0; - model->phi[i] = i*model->phi[1]; - } - else - model->phi[i] = PI*(1.0 - 2.0*(float)rand()/RAND_MAX); // note: try rand for higher harms - removed++; - } - else { - /* need to make this work thru budget of bits */ - quant_phase(&model->phi[i], -PI, PI, 3); - not_removed++; - } - #endif - } - printf("dim: %d rem %d not_rem %d\n", end-start, removed, not_removed); - -} - - -static void limit_prediction_error(struct PEXP *pexp, MODEL *model, int start, int end, float limit) -{ - int i; - float pred, pred_error, error; - - for(i=start; i<=end; i++) { - pred = pexp->phi_prev[i] + N*i*(model->Wo + pexp->Wo_prev)/2.0; - pred_error = pred - model->phi[i]; - pred_error -= TWO_PI*floor((pred_error+PI)/TWO_PI); - quant_phase(&pred_error, -limit, limit, 2); - - error = pred - pred_error - model->phi[i]; - error -= TWO_PI*floor((error+PI)/TWO_PI); - printf("%f\n", pred_error); - model->phi[i] = pred - pred_error; - } -} - - -static void quant_prediction_error(struct PEXP *pexp, MODEL *model, int start, int end, float limit) -{ - int i; - float pred, pred_error; - - for(i=start; i<=end; i++) { - pred = pexp->phi_prev[i] + N*i*(model->Wo + pexp->Wo_prev)/2.0; - pred_error = pred - model->phi[i]; - pred_error -= TWO_PI*floor((pred_error+PI)/TWO_PI); - - printf("%f\n", pred_error); - model->phi[i] = pred - pred_error; - } -} - - -static void print_sparse_pred_error(struct PEXP *pexp, MODEL *model, int start, int end, float mag_thresh) -{ - int i, index; - float mag, pred, error; - float sparse_pe[MAX_AMP]; - - mag = 0.0; - for(i=start; i<=end; i++) - mag += model->A[i]*model->A[i]; - mag = 10*log10(mag/(end-start)); - - if (mag > mag_thresh) { - for(i=0; iphi_prev[i] + N*i*(model->Wo + pexp->Wo_prev)/2.0; - error = pred - model->phi[i]; - error = atan2(sin(error),cos(error)); - - index = MAX_AMP*i*model->Wo/PI; - assert(index < MAX_AMP); - sparse_pe[index] = error; - } - - /* dump spare phase vector in polar format */ - - for(i=0; iL; m++) { - signal += model->A[m]*model->A[m]; - diff = cos(model->phi[m]) - cos(before[m]); - noise += pow(model->A[m]*diff, 2.0); - diff = sin(model->phi[m]) - sin(before[m]); - noise += pow(model->A[m]*diff, 2.0); - //printf("%f %f\n", before[m], model->phi[m]); - } - //printf("%f %f snr = %f\n", signal, noise, 10.0*log10(signal/noise)); - pexp->snr += 10.0*log10(signal/noise); -} - - -static void update_variance_calc(struct PEXP *pexp, MODEL *model, float before[]) -{ - int m; - float diff; - - for(m=1; mL; m++) { - diff = model->phi[m] - before[m]; - diff = atan2(sin(diff), cos(diff)); - pexp->var += diff*diff; - } - pexp->var_n += model->L; -} - -void print_vec(COMP cb[], int d, int e) -{ - int i,j; - - for(j=0; jWo + pexp->Wo_prev)/2.0; - best_var = 1E32; - for(Wo=0.97*Wo_est; Wo<=1.03*Wo_est; Wo+=0.001*Wo_est) { - - /* predict phase and sum differences between harmonics */ - - var = 0.0; - for(i=start; i<=end; i++) { - pred = pexp->phi_prev[i] + N*i*Wo; - error = pred - model->phi[i]; - error = atan2(sin(error),cos(error)); - var += error*error; - } - - if (var < best_var) { - best_var = var; - best_Wo = Wo; - } - } - - return best_Wo; -} - - -static void split_vq(COMP sparse_pe_out[], struct PEXP *pexp, struct codebook *vq, float weights[], COMP sparse_pe_in[]) -{ - int i, j, non_zero, vq_ind; - - //printf("\n offset %d k %d m %d j: ", vq->offset, vq->k, vq->m); - vq_ind = vq_phase(vq->cb, &sparse_pe_in[vq->offset], &weights[vq->offset], vq->k, vq->m, &pexp->vq_var); - - non_zero = 0; - for(i=0, j=vq->offset; ik; i++,j++) { - //printf("%f ", atan2(sparse_pe[i].imag, sparse_pe[i].real)); - if ((sparse_pe_in[j].real != 0.0) && (sparse_pe_in[j].imag != 0.0)) { - //printf("%d ", j); - sparse_pe_out[j] = vq->cb[vq->k * vq_ind + i]; - non_zero++; - } - } - pexp->vq_var_n += non_zero; -} - - -static void sparse_vq_pred_error(struct PEXP *pexp, - MODEL *model -) -{ - int i, index; - float pred, error, error_q_angle, best_Wo; - COMP sparse_pe_in[MAX_AMP], sparse_pe_out[MAX_AMP]; - float weights[MAX_AMP]; - COMP error_q_rect; - - best_Wo = refine_Wo(pexp, model, 1, model->L); - //best_Wo = (model->Wo + pexp->Wo_prev)/2.0; - - /* transform to sparse pred error vector */ - - for(i=0; iL; i++) { - pred = pexp->phi_prev[i] + N*i*best_Wo; - error = pred - model->phi[i]; - - index = MAX_AMP*i*model->Wo/PI; - assert(index < MAX_AMP); - sparse_pe_in[index].real = cos(error); - sparse_pe_in[index].imag = sin(error); - sparse_pe_out[index] = sparse_pe_in[index]; - weights[index] = model->A[i]; - //printf("%d ", index); - } - - /* vector quantise */ - - split_vq(sparse_pe_out, pexp, pexp->vq1, weights, sparse_pe_in); - split_vq(sparse_pe_out, pexp, pexp->vq2, weights, sparse_pe_in); - split_vq(sparse_pe_out, pexp, pexp->vq3, weights, sparse_pe_in); - split_vq(sparse_pe_out, pexp, pexp->vq4, weights, sparse_pe_in); - split_vq(sparse_pe_out, pexp, pexp->vq5, weights, sparse_pe_in); - - /* transform quantised phases back */ - - for(i=1; i<=model->L; i++) { - pred = pexp->phi_prev[i] + N*i*best_Wo; - - index = MAX_AMP*i*model->Wo/PI; - assert(index < MAX_AMP); - error_q_rect = sparse_pe_out[index]; - error_q_angle = atan2(error_q_rect.imag, error_q_rect.real); - model->phi[i] = pred - error_q_angle; - model->phi[i] = atan2(sin(model->phi[i]), cos(model->phi[i])); - } -} - - -/* - est delta (in Wo) - see if this gives us a better (smaller variance) prediction error -*/ - -static void print_pred_error_sparse_wo_correction(struct PEXP *pexp, - MODEL *model, - int start, int end, - float mag_thresh) -{ - int i, index; - float mag, pred, error[MAX_AMP], diff, c, s, delta, err; - float sparse_pe[MAX_AMP]; - - mag = 0.0; - for(i=start; i<=end; i++) - mag += model->A[i]*model->A[i]; - mag = 10*log10(mag/(end-start)); - - if (mag > mag_thresh) { - for(i=0; iphi[i] = pexp->phi_prev[i] + N*i*(model->Wo + pexp->Wo_prev)/2.0 + 0.01*N*i; - pred = pexp->phi_prev[i] + N*i*(model->Wo + pexp->Wo_prev)/2.0; - error[i] = pred - model->phi[i]; - } - - /* estimate delta Wo */ - - c = s = 0.0; - for(i=start+1; i<=end; i++) { - diff = error[i] - error[i-1]; - c += log(model->A[i])*cos(diff); - s += log(model->A[i])*sin(diff); - } - delta = atan2(s,c)/N; - //printf("delta %f\n",delta); - delta = 0; - /* now predict phases using corrected Wo */ - - for(i=start; i<=end; i++) { - pred = pexp->phi_prev[i] + N*i*(model->Wo + pexp->Wo_prev)/2.0 - N*i*delta; - err = pred - model->phi[i]; - err = atan2(sin(err),cos(err)); - - index = MAX_AMP*i*model->Wo/PI; - assert(index < MAX_AMP); - sparse_pe[index] = err; - } - - /* dump spare phase vector in polar format */ - - for(i=0; iA[i]*model->A[i]; - mag = 10*log10(mag/(end-start)); - - if (mag > mag_thresh) { - - best_Wo = refine_Wo(pexp, model, start, end); - - /* now predict phases using corrected Wo */ - - for(i=0; iphi_prev[i] + N*i*best_Wo; - err = pred - model->phi[i]; - err = atan2(sin(err),cos(err)); - - index = MAX_AMP*i*model->Wo/PI; - assert(index < MAX_AMP); - sparse_pe[index] = err; - } - - /* dump spare phase vector in polar format */ - - for(i=0; iL); - - for(i=start; i<=end; i++) { - model->phi[i] = pexp->phi_prev[i] + N*i*best_Wo; - } -} - - -/* - This functions tests theory that some bands can be combined together - due to less frequency resolution at higher frequencies. This will - reduce the amount of information we need to encode. -*/ - -void smooth_phase(struct PEXP *pexp, MODEL *model, int mode) -{ - int m, i, j, index, step, v, en, nav, st; - COMP sparse_pe_in[MAX_AMP], av; - COMP sparse_pe_out[MAX_AMP]; - COMP smoothed[MAX_AMP]; - float best_Wo, pred, err; - float weights[MAX_AMP]; - float avw, smoothed_weights[MAX_AMP]; - COMP smoothed_in[MAX_AMP], smoothed_out[MAX_AMP]; - - best_Wo = refine_Wo(pexp, model, 1, model->L); - - for(m=0; mL; m++) { - pred = pexp->phi_prev[m] + N*m*best_Wo; - err = model->phi[m] - pred; - err = atan2(sin(err),cos(err)); - - index = MAX_AMP*m*model->Wo/PI; - assert(index < MAX_AMP); - sparse_pe_in[index].real = model->A[m]*cos(err); - sparse_pe_in[index].imag = model->A[m]*sin(err); - sparse_pe_out[index] = sparse_pe_in[index]; - weights[index] = model->A[m]; - } - - /* now combine samples at high frequencies to reduce dimension */ - - step = 2; - st = 0; - for(i=st,v=0; i (MAX_AMP-1)) - en = MAX_AMP-1; - for(j=i; jvq1, smoothed_weights, smoothed_in); - for(i=0; i (MAX_AMP-1)) - en = MAX_AMP-1; - for(j=i; jL; m++) { - index = MAX_AMP*m*model->Wo/PI; - assert(index < MAX_AMP); - pred = pexp->phi_prev[m] + N*m*best_Wo; - err = atan2(sparse_pe_out[index].imag, sparse_pe_out[index].real); - model->phi[m] = pred + err; - } - -} - -/* - Another version of a functions that tests the theory that some bands - can be combined together due to less frequency resolution at higher - frequencies. This will reduce the amount of information we need to - encode. -*/ - -void smooth_phase2(struct PEXP *pexp, MODEL *model) { - float m; - float step; - int a,b,h,i; - float best_Wo, pred, err, s,c, phi1_; - - best_Wo = refine_Wo(pexp, model, 1, model->L); - - step = (float)model->L/30; - printf("\nL: %d step: %3.2f am,bm: ", model->L, step); - for(m=(float)model->L/4; m<=model->L; m+=step) { - a = floor(m); - b = floor(m+step); - if (b > model->L) b = model->L; - h = b-a; - - printf("%d,%d,(%d) ", a, b, h); - c = s = 0.0; - if (h>1) { - for(i=a; iphi_prev[i] + N*i*best_Wo; - err = model->phi[i] - pred; - c += cos(err); s += sin(err); - } - phi1_ = atan2(s,c); - for(i=a; iphi_prev[i] + N*i*best_Wo; - printf("%d: %4.3f -> ", i, model->phi[i]); - model->phi[i] = pred + phi1_; - model->phi[i] = atan2(sin(model->phi[i]),cos(model->phi[i])); - printf("%4.3f ", model->phi[i]); - } - } - } -} - - -#define MAX_BINS 40 -//static float bins[] = {2600.0, 2800.0, 3000.0, 3200.0, 3400.0, 3600.0, 3800.0, 4000.0}; -static float bins[] = {/* - - 1000.0, 1100.0, 1200.0, 1300.0, 1400.0, - 1500.0, 1600.0, 1700.0, 1800.0, 1900.0,*/ - - 2000.0, 2400.0, 2800.0, - 3000.0, 3400.0, 3600.0, 4000.0}; - -void smooth_phase3(struct PEXP *pexp, MODEL *model) { - int m, i; - int nbins; - int b; - float f, best_Wo, pred, err; - COMP av[MAX_BINS]; - - nbins = sizeof(bins)/sizeof(float); - best_Wo = refine_Wo(pexp, model, 1, model->L); - - /* clear all bins */ - - for(i=0; iL; m++) { - f = m*model->Wo*FS/TWO_PI; - if (f > bins[0]) { - - /* find bin */ - - for(i=0; i bins[i]) && (f <= bins[i+1])) - b = i; - assert(b < MAX_BINS); - - /* est predicted phase from average */ - - pred = pexp->phi_prev[m] + N*m*best_Wo; - err = model->phi[m] - pred; - av[b].real += cos(err); av[b].imag += sin(err); - } - - } - - /* use averages to est phases */ - - for(m=1; m<=model->L; m++) { - f = m*model->Wo*FS/TWO_PI; - if (f > bins[0]) { - - /* find bin */ - - for(i=0; i bins[i]) && (f <= bins[i+1])) - b = i; - assert(b < MAX_BINS); - - /* add predicted phase error to this bin */ - - printf("L %d m %d f %4.f b %d\n", model->L, m, f, b); - - pred = pexp->phi_prev[m] + N*m*best_Wo; - err = atan2(av[b].imag, av[b].real); - printf(" %d: %4.3f -> ", m, model->phi[m]); - model->phi[m] = pred + err; - model->phi[m] = atan2(sin(model->phi[m]),cos(model->phi[m])); - printf("%4.3f\n", model->phi[m]); - } - } - printf("\n"); -} - - -/* - Try to code the phase of the largest amplitude in each band. Randomise the - phase of the other harmonics. The theory is that only the largest harmonic - will be audible. -*/ - -void cb_phase1(struct PEXP *pexp, MODEL *model) { - int m, i; - int nbins; - int b; - float f, best_Wo; - float max_val[MAX_BINS]; - int max_ind[MAX_BINS]; - - nbins = sizeof(bins)/sizeof(float); - best_Wo = refine_Wo(pexp, model, 1, model->L); - - for(i=0; iL; m++) { - f = m*model->Wo*FS/TWO_PI; - if (f > bins[0]) { - - /* find bin */ - - for(i=0; i bins[i]) && (f <= bins[i+1])) - b = i; - assert(b < MAX_BINS); - - if (model->A[m] > max_val[b]) { - max_val[b] = model->A[m]; - max_ind[b] = m; - } - } - - } - - /* randomise phase of other harmonics */ - - for(m=1; m<=model->L; m++) { - f = m*model->Wo*FS/TWO_PI; - if (f > bins[0]) { - - /* find bin */ - - for(i=0; i bins[i]) && (f <= bins[i+1])) - b = i; - assert(b < MAX_BINS); - - if (m != max_ind[b]) - model->phi[m] = pexp->phi_prev[m] + N*m*best_Wo; - } - } -} - - -/* - Theory is only the phase of the envelope of signal matters within a - Critical Band. So we estimate the position of an impulse that - approximates the envelope of the signal. -*/ - -void cb_phase2(struct PEXP *pexp, MODEL *model) { - int st, m, i, a, b, step; - float diff,w,c,s,phi1_; - float A[MAX_AMP]; - - for(m=1; m<=model->L; m++) { - A[m] = model->A[m]; - model->A[m] = 0; - } - - st = 2*model->L/4; - step = 3; - model->phi[1] = pexp->phi_prev[1] + (pexp->Wo_prev+model->Wo)*N/2.0; - - printf("L=%d ", model->L); - for(m=st; m model->L) - b = model->L; - - c = s = 0; - for(i=a; iphi[i+1] - model->phi[i]; - //w = (model->A[i+1] + model->A[i])/2; - w = 1.0; - c += w*cos(diff); s += w*sin(diff); - } - phi1_ = atan2(s,c); - printf("replacing: "); - for(i=a; iphi[i] = i*phi1_; - //model->phi[i] = i*model->phi[1]; - //model->phi[i] = m*(pexp->Wo_prev+model->Wo)*N/2.0; - model->A[m] = A[m]; - printf("%d ", i); - } - printf(" . "); - } - printf("\n"); -} - - -static void smooth_phase4(MODEL *model) { - int m; - float phi_m, phi_m_1; - - if (model->L > 25) { - printf("\nL %d\n", model->L); - for(m=model->L/2; m<=model->L; m+=2) { - if ((m+1) <= model->L) { - phi_m = (model->phi[m] - model->phi[m+1])/2.0; - phi_m_1 = (model->phi[m+1] - model->phi[m])/2.0; - model->phi[m] = phi_m; - model->phi[m+1] = phi_m_1; - printf("%d %4.3f %4.3f ", m, phi_m, phi_m_1); - } - } - } - -} - -/* try repeating last frame, just advance phases to account for time shift */ - -static void repeat_phases(struct PEXP *pexp, MODEL *model) { - int m; - - *model = pexp->prev_model; - for(m=1; m<=model->L; m++) - model->phi[m] += N*m*model->Wo; - -} - -/*---------------------------------------------------------------------------*\ - - phase_experiment() - - Phase quantisation experiments. - -\*---------------------------------------------------------------------------*/ - -void phase_experiment(struct PEXP *pexp, MODEL *model, char *arg) { - int m; - float before[MAX_AMP], best_Wo; - - assert(pexp != NULL); - memcpy(before, &model->phi[0], sizeof(float)*MAX_AMP); - - if (strcmp(arg,"q3") == 0) { - quant_phases(model, 1, model->L, 3); - update_snr_calc(pexp, model, before); - update_variance_calc(pexp, model, before); - } - - if (strcmp(arg,"dec2") == 0) { - if ((pexp->frames % 2) != 0) { - predict_phases(pexp, model, 1, model->L); - update_snr_calc(pexp, model, before); - update_variance_calc(pexp, model, before); - } - } - - if (strcmp(arg,"repeat") == 0) { - if ((pexp->frames % 2) != 0) { - repeat_phases(pexp, model); - update_snr_calc(pexp, model, before); - update_variance_calc(pexp, model, before); - } - } - - if (strcmp(arg,"vq") == 0) { - sparse_vq_pred_error(pexp, model); - update_snr_calc(pexp, model, before); - update_variance_calc(pexp, model, before); - } - - if (strcmp(arg,"pred") == 0) - predict_phases_state(pexp, model, 1, model->L); - - if (strcmp(arg,"pred1k") == 0) - predict_phases(pexp, model, 1, model->L/4); - - if (strcmp(arg,"smooth") == 0) { - smooth_phase(pexp, model,0); - update_snr_calc(pexp, model, before); - } - if (strcmp(arg,"smoothtrain") == 0) - smooth_phase(pexp, model,1); - if (strcmp(arg,"smoothvq") == 0) { - smooth_phase(pexp, model,2); - update_snr_calc(pexp, model, before); - } - - if (strcmp(arg,"smooth2") == 0) - smooth_phase2(pexp, model); - if (strcmp(arg,"smooth3") == 0) - smooth_phase3(pexp, model); - if (strcmp(arg,"smooth4") == 0) - smooth_phase4(model); - if (strcmp(arg,"vqsmooth3") == 0) { - sparse_vq_pred_error(pexp, model); - smooth_phase3(pexp, model); - } - - if (strcmp(arg,"cb1") == 0) { - cb_phase1(pexp, model); - update_snr_calc(pexp, model, before); - } - - if (strcmp(arg,"top") == 0) { - //top_amp(pexp, model, 1, model->L/4, 4, 1); - //top_amp(pexp, model, model->L/4, model->L/3, 4, 1); - //top_amp(pexp, model, model->L/3+1, model->L/2, 4, 1); - //top_amp(pexp, model, model->L/2, model->L, 6, 1); - //rand_phases(model, model->L/2, 3*model->L/4); - //struct_phases(pexp, model, model->L/2, 3*model->L/4); - //update_snr_calc(pexp, model, before); - } - - if (strcmp(arg,"pred23") == 0) { - predict_phases2(pexp, model, model->L/2, model->L); - update_snr_calc(pexp, model, before); - } - - if (strcmp(arg,"struct23") == 0) { - struct_phases(pexp, model, model->L/2, 3*model->L/4 ); - update_snr_calc(pexp, model, before); - } - - if (strcmp(arg,"addnoise") == 0) { - int m; - float max; - - max = 0; - for(m=1; m<=model->L; m++) - if (model->A[m] > max) - max = model->A[m]; - max = 20.0*log10(max); - for(m=1; m<=model->L; m++) - if (20.0*log10(model->A[m]) < (max-20)) { - model->phi[m] += (PI/4)*(1.0 -2.0*rand()/RAND_MAX); - //printf("m %d\n", m); - } - } - - /* normalise phases */ - - for(m=1; m<=model->L; m++) - model->phi[m] = atan2(sin(model->phi[m]), cos(model->phi[m])); - - /* update states */ - - //best_Wo = refine_Wo(pexp, model, model->L/2, model->L); - pexp->phi1 += N*model->Wo; - - for(m=1; m<=model->L; m++) - pexp->phi_prev[m] = model->phi[m]; - pexp->Wo_prev = model->Wo; - pexp->frames++; - pexp->prev_model = *model; -} - -#ifdef OLD_STUFF - //quant_phases(model, 1, model->L, 3); - //update_variance_calc(pexp, model, before); - //print_sparse_pred_error(pexp, model, 1, model->L, 40.0); - - //sparse_vq_pred_error(pexp, model); - - //quant_phases(model, model->L/4+1, model->L, 3); - - //predict_phases1(pexp, model, 1, model->L/4); - //quant_phases(model, model->L/4+1, model->L, 3); - - //quant_phases(model, 1, model->L/8, 3); - - //update_snr_calc(pexp, model, before); - //update_variance_calc(pexp, model, before); - - //fixed_bits_per_frame(pexp, model, 40); - //struct_phases(pexp, model, 1, model->L/4); - //rand_phases(model, 10, model->L); - //for(m=1; m<=model->L; m++) - // model->A[m] = 0.0; - //model->A[model->L/2] = 1000; - //repeat_phases(model, 20); - //predict_phases(pexp, model, 1, model->L/4); - //quant_phases(model, 1, 10, 3); - //quant_phases(model, 10, 20, 2); - //repeat_phases(model, 20); - //rand_phases(model, 3*model->L/4, model->L); - // print_phi1_pred_error(model, 1, model->L); - //predict_phases(pexp, model, 1, model->L/4); - //first_order_band(model, model->L/4, model->L/2); - //first_order_band(model, model->L/2, 3*model->L/4); - //if (fabs(model->Wo - pexp->Wo_prev)< 0.1*model->Wo) - - //print_pred_error(pexp, model, 1, model->L, 40.0); - //print_sparse_pred_error(pexp, model, 1, model->L, 40.0); - - //phi1_est = est_phi1(model, 1, model->L/4); - //print_phi1_pred_error(model, 1, model->L/4); - - //first_order_band(model, 1, model->L/4, phi1_est); - //sub_linear(model, 1, model->L/4, phi1_est); - - //top_amp(pexp, model, 1, model->L/4, 4); - //top_amp(pexp, model, model->L/4, model->L/2, 4); - - //first_order_band(model, 1, model->L/4, phi1_est); - //first_order_band(model, model->L/4, model->L/2, phi1_est); - - //if (fabs(model->Wo - pexp->Wo_prev) > 0.2*model->Wo) - // rand_phases(model, model->L/2, model->L); - - //top_amp(pexp, model, 1, model->L/4, 4); - //top_amp(pexp, model, model->L/4, model->L/2, 8); - //top_amp(pexp, model, model->L/4+1, model->L/2, 10, 1); - //top_amp(pexp, model, 1, model->L/4, 10, 1); - //top_amp(pexp, model, model->L/4+1, 3*model->L/4, 10, 1); - //top_amp(pexp, model, 1, 3*model->L/4, 20, 1); - - #ifdef REAS_CAND1 - predict_phases(pexp, model, 1, model->L/4); - top_amp(pexp, model, model->L/4+1, 3*model->L/4, 10, 1); - rand_phases(model, 3*model->L/4+1, model->L); - #endif - - #ifdef REAS_CAND2 - if ((pexp->frames % 2) == 0) { - //printf("quant\n"); - predict_phases(pexp, model, 1, model->L/4); - //top_amp(pexp, model, model->L/4+1, 3*model->L/4, 20, 1); - top_amp(pexp, model, model->L/4+1, 7*model->L/8, 20, 1); - rand_phases(model, 7*model->L/8+1, model->L); - } - else { - //printf("predict\n"); - predict_phases(pexp, model, 1, model->L); - } - #endif - - //#define REAS_CAND3 - #ifdef REAS_CAND3 - if ((pexp->frames % 3) != 0) { - printf("pred\n"); - predict_phases(pexp, model, 1, model->L); - } - else { - predict_phases(pexp, model, 1, model->L/4); - fixed_bits_per_frame(pexp, model, model->L/4+1, 60); - } - #endif - //predict_phases(pexp, model, model->L/4, model->L); - - - //print_pred_error(pexp, model, 1, model->L); - //limit_prediction_error(pexp, model, model->L/2, model->L, PI/2); -#endif diff --git a/codec2-dev/src/phase.h b/codec2-dev/src/phase.h index 8619c3fc..367948df 100644 --- a/codec2-dev/src/phase.h +++ b/codec2-dev/src/phase.h @@ -30,13 +30,10 @@ #include "kiss_fft.h" -void phase_synth_zero_order(kiss_fft_cfg fft_dec_cfg, MODEL *model, float aks[], - float *ex_phase, int order); - -struct PEXP; - -struct PEXP * phase_experiment_create(); -void phase_experiment_destroy(struct PEXP *pexp); -void phase_experiment(struct PEXP *pexp, MODEL *model, char *arg); +void phase_synth_zero_order(kiss_fft_cfg fft_dec_cfg, + MODEL *model, + float aks[], + float *ex_phase, + int order); #endif diff --git a/codec2-dev/src/phaseexp.c b/codec2-dev/src/phaseexp.c new file mode 100644 index 00000000..2cd88b0e --- /dev/null +++ b/codec2-dev/src/phaseexp.c @@ -0,0 +1,1575 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: phaseexp.c + AUTHOR......: David Rowe + DATE CREATED: June 2012 + + Experimental functions for quantising, modelling and synthesising phase. + +\*---------------------------------------------------------------------------*/ + +/* + Copyright (C) 2012 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 . +*/ + +#include "defines.h" +#include "phase.h" +#include "kiss_fft.h" +#include "comp.h" +#include "glottal.c" + +#include +#include +#include +#include +#include + +/* Bruce Perens' funcs to load codebook files */ + +struct codebook { + unsigned int k; + unsigned int log2m; + unsigned int m; + COMP *cb; + unsigned int offset; +}; + +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"; + +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(const char * name) +{ + FILE *file; + char line[2048]; + char *cursor = line; + struct codebook *b = malloc(sizeof(struct codebook)); + int i; + int size; + float angle; + + file = fopen(name, "rt"); + assert(file != NULL); + + *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 = (COMP *)malloc(size * sizeof(COMP)); + + for ( i = 0; i < size; i++ ) { + angle = get_float(file, name, &cursor, line, sizeof(line)); + b->cb[i].real = cos(angle); + b->cb[i].imag = sin(angle); + } + + fclose(file); + + return b; +} + + +/* states for phase experiments */ + +struct PEXP { + float phi1; + float phi_prev[MAX_AMP]; + float Wo_prev; + int frames; + float snr; + float var; + int var_n; + struct codebook *vq1,*vq2,*vq3,*vq4,*vq5; + float vq_var; + int vq_var_n; + MODEL prev_model; + int state; +}; + + +/*---------------------------------------------------------------------------* \ + + phase_experiment_create() + + Inits states for phase quantisation experiments. + +\*---------------------------------------------------------------------------*/ + +struct PEXP * phase_experiment_create() { + struct PEXP *pexp; + int i; + + pexp = (struct PEXP *)malloc(sizeof(struct PEXP)); + assert (pexp != NULL); + + pexp->phi1 = 0; + for(i=0; iphi_prev[i] = 0.0; + pexp->Wo_prev = 0.0; + pexp->frames = 0; + pexp->snr = 0.0; + pexp->var = 0.0; + pexp->var_n = 0; + + /* smoothed 10th order for 1st 1 khz */ + //pexp->vq1 = load("../unittest/ph1_10_1024.txt"); + //pexp->vq1->offset = 0; + + /* load experimental phase VQ */ + + //pexp->vq1 = load("../unittest/testn1_20_1024.txt"); + pexp->vq1 = load("../unittest/test.txt"); + //pexp->vq2 = load("../unittest/testn21_40_1024.txt"); + pexp->vq2 = load("../unittest/test11_20_1024.txt"); + pexp->vq3 = load("../unittest/test21_30_1024.txt"); + pexp->vq4 = load("../unittest/test31_40_1024.txt"); + pexp->vq5 = load("../unittest/test41_60_1024.txt"); + pexp->vq1->offset = 0; + pexp->vq2->offset = 10; + pexp->vq3->offset = 20; + pexp->vq4->offset = 30; + pexp->vq5->offset = 40; + + pexp->vq_var = 0.0; + pexp->vq_var_n = 0; + + pexp->state = 0; + + return pexp; +} + + +/*---------------------------------------------------------------------------* \ + + phase_experiment_destroy() + +\*---------------------------------------------------------------------------*/ + +void phase_experiment_destroy(struct PEXP *pexp) { + assert(pexp != NULL); + if (pexp->snr != 0.0) + printf("snr: %4.2f dB\n", pexp->snr/pexp->frames); + if (pexp->var != 0.0) + printf("var...: %4.3f std dev...: %4.3f (%d non zero phases)\n", + pexp->var/pexp->var_n, sqrt(pexp->var/pexp->var_n), pexp->var_n); + if (pexp->vq_var != 0.0) + printf("vq var: %4.3f vq std dev: %4.3f (%d non zero phases)\n", + pexp->vq_var/pexp->vq_var_n, sqrt(pexp->vq_var/pexp->vq_var_n), pexp->vq_var_n); + free(pexp); +} + + +/*---------------------------------------------------------------------------* \ + + Various test and experimental functions ................ + +\*---------------------------------------------------------------------------*/ + +/* Bubblesort to find highest amplitude harmonics */ + +struct AMPINDEX { + float amp; + int index; +}; + +static void bubbleSort(struct AMPINDEX numbers[], int array_size) +{ + int i, j; + struct AMPINDEX temp; + + for (i = (array_size - 1); i > 0; i--) + { + for (j = 1; j <= i; j++) + { + //printf("i %d j %d %f %f \n", i, j, numbers[j-1].amp, numbers[j].amp); + if (numbers[j-1].amp < numbers[j].amp) + { + temp = numbers[j-1]; + numbers[j-1] = numbers[j]; + numbers[j] = temp; + } + } + } +} + + +static void print_pred_error(struct PEXP *pexp, MODEL *model, int start, int end, float mag_thresh) { + int i; + float mag; + + mag = 0.0; + for(i=start; i<=end; i++) + mag += model->A[i]*model->A[i]; + mag = 10*log10(mag/(end-start)); + + if (mag > mag_thresh) { + for(i=start; i<=end; i++) { + float pred = pexp->phi_prev[i] + N*i*(model->Wo + pexp->Wo_prev)/2.0; + float err = pred - model->phi[i]; + err = atan2(sin(err),cos(err)); + printf("%f\n",err); + } + //printf("\n"); + } + +} + + +static void predict_phases(struct PEXP *pexp, MODEL *model, int start, int end) { + int i; + + for(i=start; i<=end; i++) { + model->phi[i] = pexp->phi_prev[i] + N*i*model->Wo; + } + +} +static float refine_Wo(struct PEXP *pexp, + MODEL *model, + int start, + int end); + +/* Fancy state based phase prediction. Actually works OK on most utterances, + but could use some tuning. Breaks down a bit on mmt1. */ + +static void predict_phases_state(struct PEXP *pexp, MODEL *model, int start, int end) { + int i, next_state; + float best_Wo, dWo; + + //best_Wo = refine_Wo(pexp, model, start, end); + //best_Wo = (model->Wo + pexp->Wo_prev)/2.0; + best_Wo = model->Wo; + + dWo = fabs(model->Wo - pexp->Wo_prev)/model->Wo; + next_state = pexp->state; + switch(pexp->state) { + case 0: + if (dWo < 0.1) { + /* UV -> V transition, so start with phases in lock. They will + drift a bit over voiced track which is kinda what we want, so + we don't get clicky speech. + */ + next_state = 1; + for(i=start; i<=end; i++) + pexp->phi_prev[i] = i*pexp->phi1; + } + + break; + case 1: + if (dWo > 0.1) + next_state = 0; + break; + } + pexp->state = next_state; + + if (pexp->state == 0) + for(i=start; i<=end; i++) { + model->phi[i] = PI*(1.0 - 2.0*rand()/RAND_MAX); + } + else + for(i=start; i<=end; i++) { + model->phi[i] = pexp->phi_prev[i] + N*i*best_Wo; + } + printf("state %d\n", pexp->state); +} + +static void struct_phases(struct PEXP *pexp, MODEL *model, int start, int end) { + int i; + + for(i=start; i<=end; i++) + model->phi[i] = pexp->phi1*i; + +} + + +static void predict_phases2(struct PEXP *pexp, MODEL *model, int start, int end) { + int i; + float pred, str, diff; + + for(i=start; i<=end; i++) { + pred = pexp->phi_prev[i] + N*i*model->Wo; + str = pexp->phi1*i; + diff = str - pred; + diff = atan2(sin(diff), cos(diff)); + if (diff > 0) + pred += PI/16; + else + pred -= PI/16; + model->phi[i] = pred; + } + +} + +static void skip_phases(struct PEXP *pexp, MODEL *model, int start, int end) { + int i; + + for(i=start; i<=end; i+=2) + model->phi[i] = model->phi[i-1] - model->phi[i-2]; + +} + +static void rand_phases(MODEL *model, int start, int end) { + int i; + + for(i=start; i<=end; i++) + model->phi[i] = PI*(1.0 - 2.0*(float)rand()/RAND_MAX); + +} + +static void quant_phase(float *phase, float min, float max, int bits) { + int levels = 1 << bits; + int index; + float norm, step; + + norm = (*phase - min)/(max - min); + index = floor(levels*norm); + + //printf("phase %f norm %f index %d ", *phase, norm, index); + if (index < 0 ) index = 0; + if (index > (levels-1)) index = levels-1; + //printf("index %d ", index); + step = (max - min)/levels; + *phase = min + step*index + 0.5*step; + //printf("step %f phase %f\n", step, *phase); +} + +static void quant_phases(MODEL *model, int start, int end, int bits) { + int i; + + for(i=start; i<=end; i++) { + quant_phase(&model->phi[i], -PI, PI, bits); + } +} + +static void fixed_bits_per_frame(struct PEXP *pexp, MODEL *model, int m, int budget) { + int res, finished; + + res = 3; + finished = 0; + + while(!finished) { + if (m > model->L/2) + res = 2; + if (((budget - res) < 0) || (m > model->L)) + finished = 1; + else { + quant_phase(&model->phi[m], -PI, PI, res); + budget -= res; + m++; + } + } + printf("m: %d L: %d budget: %d\n", m, model->L, budget); + predict_phases(pexp, model, m, model->L); + //rand_phases(model, m, model->L); +} + +/* used to plot histogram of quantisation error, for 3 bits, 8 levels, + should be uniform between +/- PI/8 */ + +static void check_phase_quant(MODEL *model, float tol) +{ + int m; + float phi_before[MAX_AMP]; + + for(m=1; m<=model->L; m++) + phi_before[m] = model->phi[m]; + + quant_phases(model, 1, model->L, 3); + + for(m=1; m<=model->L; m++) { + float err = phi_before[m] - model->phi[m]; + printf("%f\n", err); + if (fabs(err) > tol) + exit(0); + } +} + + +static float est_phi1(MODEL *model, int start, int end) +{ + int m; + float delta, s, c, phi1_est; + + if (end > model->L) + end = model->L; + + s = c = 0.0; + for(m=start; mphi[m+1] - model->phi[m]; + s += sin(delta); + c += cos(delta); + } + + phi1_est = atan2(s,c); + + return phi1_est; +} + +static void print_phi1_pred_error(MODEL *model, int start, int end) +{ + int m; + float phi1_est; + + phi1_est = est_phi1(model, start, end); + + for(m=start; mphi[m+1] - model->phi[m] - phi1_est; + err = atan2(sin(err),cos(err)); + printf("%f\n", err); + } +} + + +static void first_order_band(MODEL *model, int start, int end, float phi1_est) +{ + int m; + float pred_err, av_pred_err; + float c,s; + + s = c = 0.0; + for(m=start; mphi[m] - phi1_est*m; + s += sin(pred_err); + c += cos(pred_err); + } + + av_pred_err = atan2(s,c); + for(m=start; mphi[m] = av_pred_err + phi1_est*m; + model->phi[m] = atan2(sin(model->phi[m]), cos(model->phi[m])); + } + +} + + +static void sub_linear(MODEL *model, int start, int end, float phi1_est) +{ + int m; + + for(m=start; mphi[m] = m*phi1_est; + } +} + + +static void top_amp(struct PEXP *pexp, MODEL *model, int start, int end, int n_harm, int pred) +{ + int removed = 0, not_removed = 0; + int top, i, j; + struct AMPINDEX sorted[MAX_AMP]; + + /* sort into ascending order of amplitude */ + + printf("\n"); + for(i=start,j=0; iA[i]; + sorted[j].index = i; + printf("%f ", model->A[i]); + } + bubbleSort(sorted, end-start); + + printf("\n"); + for(j=0; jA[i] == sorted[j].amp) { + top = 1; + assert(i == sorted[j].index); + } + } + + #define ALTTOP + #ifdef ALTTOP + model->phi[i] = 0.0; /* make sure */ + if (top) { + model->phi[i] = i*pexp->phi1; + removed++; + } + else { + model->phi[i] = PI*(1.0 - 2.0*(float)rand()/RAND_MAX); // note: try rand for higher harms + removed++; + } + #else + if (!top) { + model->phi[i] = 0.0; /* make sure */ + if (pred) { + //model->phi[i] = pexp->phi_prev[i] + i*N*(model->Wo + pexp->Wo_prev)/2.0; + model->phi[i] = i*model->phi[1]; + } + else + model->phi[i] = PI*(1.0 - 2.0*(float)rand()/RAND_MAX); // note: try rand for higher harms + removed++; + } + else { + /* need to make this work thru budget of bits */ + quant_phase(&model->phi[i], -PI, PI, 3); + not_removed++; + } + #endif + } + printf("dim: %d rem %d not_rem %d\n", end-start, removed, not_removed); + +} + + +static void limit_prediction_error(struct PEXP *pexp, MODEL *model, int start, int end, float limit) +{ + int i; + float pred, pred_error, error; + + for(i=start; i<=end; i++) { + pred = pexp->phi_prev[i] + N*i*(model->Wo + pexp->Wo_prev)/2.0; + pred_error = pred - model->phi[i]; + pred_error -= TWO_PI*floor((pred_error+PI)/TWO_PI); + quant_phase(&pred_error, -limit, limit, 2); + + error = pred - pred_error - model->phi[i]; + error -= TWO_PI*floor((error+PI)/TWO_PI); + printf("%f\n", pred_error); + model->phi[i] = pred - pred_error; + } +} + + +static void quant_prediction_error(struct PEXP *pexp, MODEL *model, int start, int end, float limit) +{ + int i; + float pred, pred_error; + + for(i=start; i<=end; i++) { + pred = pexp->phi_prev[i] + N*i*(model->Wo + pexp->Wo_prev)/2.0; + pred_error = pred - model->phi[i]; + pred_error -= TWO_PI*floor((pred_error+PI)/TWO_PI); + + printf("%f\n", pred_error); + model->phi[i] = pred - pred_error; + } +} + + +static void print_sparse_pred_error(struct PEXP *pexp, MODEL *model, int start, int end, float mag_thresh) +{ + int i, index; + float mag, pred, error; + float sparse_pe[MAX_AMP]; + + mag = 0.0; + for(i=start; i<=end; i++) + mag += model->A[i]*model->A[i]; + mag = 10*log10(mag/(end-start)); + + if (mag > mag_thresh) { + for(i=0; iphi_prev[i] + N*i*(model->Wo + pexp->Wo_prev)/2.0; + error = pred - model->phi[i]; + error = atan2(sin(error),cos(error)); + + index = MAX_AMP*i*model->Wo/PI; + assert(index < MAX_AMP); + sparse_pe[index] = error; + } + + /* dump spare phase vector in polar format */ + + for(i=0; iL; m++) { + signal += model->A[m]*model->A[m]; + diff = cos(model->phi[m]) - cos(before[m]); + noise += pow(model->A[m]*diff, 2.0); + diff = sin(model->phi[m]) - sin(before[m]); + noise += pow(model->A[m]*diff, 2.0); + //printf("%f %f\n", before[m], model->phi[m]); + } + //printf("%f %f snr = %f\n", signal, noise, 10.0*log10(signal/noise)); + pexp->snr += 10.0*log10(signal/noise); +} + + +static void update_variance_calc(struct PEXP *pexp, MODEL *model, float before[]) +{ + int m; + float diff; + + for(m=1; mL; m++) { + diff = model->phi[m] - before[m]; + diff = atan2(sin(diff), cos(diff)); + pexp->var += diff*diff; + } + pexp->var_n += model->L; +} + +void print_vec(COMP cb[], int d, int e) +{ + int i,j; + + for(j=0; jWo + pexp->Wo_prev)/2.0; + best_var = 1E32; + for(Wo=0.97*Wo_est; Wo<=1.03*Wo_est; Wo+=0.001*Wo_est) { + + /* predict phase and sum differences between harmonics */ + + var = 0.0; + for(i=start; i<=end; i++) { + pred = pexp->phi_prev[i] + N*i*Wo; + error = pred - model->phi[i]; + error = atan2(sin(error),cos(error)); + var += error*error; + } + + if (var < best_var) { + best_var = var; + best_Wo = Wo; + } + } + + return best_Wo; +} + + +static void split_vq(COMP sparse_pe_out[], struct PEXP *pexp, struct codebook *vq, float weights[], COMP sparse_pe_in[]) +{ + int i, j, non_zero, vq_ind; + + //printf("\n offset %d k %d m %d j: ", vq->offset, vq->k, vq->m); + vq_ind = vq_phase(vq->cb, &sparse_pe_in[vq->offset], &weights[vq->offset], vq->k, vq->m, &pexp->vq_var); + + non_zero = 0; + for(i=0, j=vq->offset; ik; i++,j++) { + //printf("%f ", atan2(sparse_pe[i].imag, sparse_pe[i].real)); + if ((sparse_pe_in[j].real != 0.0) && (sparse_pe_in[j].imag != 0.0)) { + //printf("%d ", j); + sparse_pe_out[j] = vq->cb[vq->k * vq_ind + i]; + non_zero++; + } + } + pexp->vq_var_n += non_zero; +} + + +static void sparse_vq_pred_error(struct PEXP *pexp, + MODEL *model +) +{ + int i, index; + float pred, error, error_q_angle, best_Wo; + COMP sparse_pe_in[MAX_AMP], sparse_pe_out[MAX_AMP]; + float weights[MAX_AMP]; + COMP error_q_rect; + + best_Wo = refine_Wo(pexp, model, 1, model->L); + //best_Wo = (model->Wo + pexp->Wo_prev)/2.0; + + /* transform to sparse pred error vector */ + + for(i=0; iL; i++) { + pred = pexp->phi_prev[i] + N*i*best_Wo; + error = pred - model->phi[i]; + + index = MAX_AMP*i*model->Wo/PI; + assert(index < MAX_AMP); + sparse_pe_in[index].real = cos(error); + sparse_pe_in[index].imag = sin(error); + sparse_pe_out[index] = sparse_pe_in[index]; + weights[index] = model->A[i]; + //printf("%d ", index); + } + + /* vector quantise */ + + split_vq(sparse_pe_out, pexp, pexp->vq1, weights, sparse_pe_in); + split_vq(sparse_pe_out, pexp, pexp->vq2, weights, sparse_pe_in); + split_vq(sparse_pe_out, pexp, pexp->vq3, weights, sparse_pe_in); + split_vq(sparse_pe_out, pexp, pexp->vq4, weights, sparse_pe_in); + split_vq(sparse_pe_out, pexp, pexp->vq5, weights, sparse_pe_in); + + /* transform quantised phases back */ + + for(i=1; i<=model->L; i++) { + pred = pexp->phi_prev[i] + N*i*best_Wo; + + index = MAX_AMP*i*model->Wo/PI; + assert(index < MAX_AMP); + error_q_rect = sparse_pe_out[index]; + error_q_angle = atan2(error_q_rect.imag, error_q_rect.real); + model->phi[i] = pred - error_q_angle; + model->phi[i] = atan2(sin(model->phi[i]), cos(model->phi[i])); + } +} + + +/* + est delta (in Wo) + see if this gives us a better (smaller variance) prediction error +*/ + +static void print_pred_error_sparse_wo_correction(struct PEXP *pexp, + MODEL *model, + int start, int end, + float mag_thresh) +{ + int i, index; + float mag, pred, error[MAX_AMP], diff, c, s, delta, err; + float sparse_pe[MAX_AMP]; + + mag = 0.0; + for(i=start; i<=end; i++) + mag += model->A[i]*model->A[i]; + mag = 10*log10(mag/(end-start)); + + if (mag > mag_thresh) { + for(i=0; iphi[i] = pexp->phi_prev[i] + N*i*(model->Wo + pexp->Wo_prev)/2.0 + 0.01*N*i; + pred = pexp->phi_prev[i] + N*i*(model->Wo + pexp->Wo_prev)/2.0; + error[i] = pred - model->phi[i]; + } + + /* estimate delta Wo */ + + c = s = 0.0; + for(i=start+1; i<=end; i++) { + diff = error[i] - error[i-1]; + c += log(model->A[i])*cos(diff); + s += log(model->A[i])*sin(diff); + } + delta = atan2(s,c)/N; + //printf("delta %f\n",delta); + delta = 0; + /* now predict phases using corrected Wo */ + + for(i=start; i<=end; i++) { + pred = pexp->phi_prev[i] + N*i*(model->Wo + pexp->Wo_prev)/2.0 - N*i*delta; + err = pred - model->phi[i]; + err = atan2(sin(err),cos(err)); + + index = MAX_AMP*i*model->Wo/PI; + assert(index < MAX_AMP); + sparse_pe[index] = err; + } + + /* dump spare phase vector in polar format */ + + for(i=0; iA[i]*model->A[i]; + mag = 10*log10(mag/(end-start)); + + if (mag > mag_thresh) { + + best_Wo = refine_Wo(pexp, model, start, end); + + /* now predict phases using corrected Wo */ + + for(i=0; iphi_prev[i] + N*i*best_Wo; + err = pred - model->phi[i]; + err = atan2(sin(err),cos(err)); + + index = MAX_AMP*i*model->Wo/PI; + assert(index < MAX_AMP); + sparse_pe[index] = err; + } + + /* dump spare phase vector in polar format */ + + for(i=0; iL); + + for(i=start; i<=end; i++) { + model->phi[i] = pexp->phi_prev[i] + N*i*best_Wo; + } +} + + +/* + This functions tests theory that some bands can be combined together + due to less frequency resolution at higher frequencies. This will + reduce the amount of information we need to encode. +*/ + +void smooth_phase(struct PEXP *pexp, MODEL *model, int mode) +{ + int m, i, j, index, step, v, en, nav, st; + COMP sparse_pe_in[MAX_AMP], av; + COMP sparse_pe_out[MAX_AMP]; + COMP smoothed[MAX_AMP]; + float best_Wo, pred, err; + float weights[MAX_AMP]; + float avw, smoothed_weights[MAX_AMP]; + COMP smoothed_in[MAX_AMP], smoothed_out[MAX_AMP]; + + best_Wo = refine_Wo(pexp, model, 1, model->L); + + for(m=0; mL; m++) { + pred = pexp->phi_prev[m] + N*m*best_Wo; + err = model->phi[m] - pred; + err = atan2(sin(err),cos(err)); + + index = MAX_AMP*m*model->Wo/PI; + assert(index < MAX_AMP); + sparse_pe_in[index].real = model->A[m]*cos(err); + sparse_pe_in[index].imag = model->A[m]*sin(err); + sparse_pe_out[index] = sparse_pe_in[index]; + weights[index] = model->A[m]; + } + + /* now combine samples at high frequencies to reduce dimension */ + + step = 2; + st = 0; + for(i=st,v=0; i (MAX_AMP-1)) + en = MAX_AMP-1; + for(j=i; jvq1, smoothed_weights, smoothed_in); + for(i=0; i (MAX_AMP-1)) + en = MAX_AMP-1; + for(j=i; jL; m++) { + index = MAX_AMP*m*model->Wo/PI; + assert(index < MAX_AMP); + pred = pexp->phi_prev[m] + N*m*best_Wo; + err = atan2(sparse_pe_out[index].imag, sparse_pe_out[index].real); + model->phi[m] = pred + err; + } + +} + +/* + Another version of a functions that tests the theory that some bands + can be combined together due to less frequency resolution at higher + frequencies. This will reduce the amount of information we need to + encode. +*/ + +void smooth_phase2(struct PEXP *pexp, MODEL *model) { + float m; + float step; + int a,b,h,i; + float best_Wo, pred, err, s,c, phi1_; + + best_Wo = refine_Wo(pexp, model, 1, model->L); + + step = (float)model->L/30; + printf("\nL: %d step: %3.2f am,bm: ", model->L, step); + for(m=(float)model->L/4; m<=model->L; m+=step) { + a = floor(m); + b = floor(m+step); + if (b > model->L) b = model->L; + h = b-a; + + printf("%d,%d,(%d) ", a, b, h); + c = s = 0.0; + if (h>1) { + for(i=a; iphi_prev[i] + N*i*best_Wo; + err = model->phi[i] - pred; + c += cos(err); s += sin(err); + } + phi1_ = atan2(s,c); + for(i=a; iphi_prev[i] + N*i*best_Wo; + printf("%d: %4.3f -> ", i, model->phi[i]); + model->phi[i] = pred + phi1_; + model->phi[i] = atan2(sin(model->phi[i]),cos(model->phi[i])); + printf("%4.3f ", model->phi[i]); + } + } + } +} + + +#define MAX_BINS 40 +//static float bins[] = {2600.0, 2800.0, 3000.0, 3200.0, 3400.0, 3600.0, 3800.0, 4000.0}; +static float bins[] = {/* + + 1000.0, 1100.0, 1200.0, 1300.0, 1400.0, + 1500.0, 1600.0, 1700.0, 1800.0, 1900.0,*/ + + 2000.0, 2400.0, 2800.0, + 3000.0, 3400.0, 3600.0, 4000.0}; + +void smooth_phase3(struct PEXP *pexp, MODEL *model) { + int m, i; + int nbins; + int b; + float f, best_Wo, pred, err; + COMP av[MAX_BINS]; + + nbins = sizeof(bins)/sizeof(float); + best_Wo = refine_Wo(pexp, model, 1, model->L); + + /* clear all bins */ + + for(i=0; iL; m++) { + f = m*model->Wo*FS/TWO_PI; + if (f > bins[0]) { + + /* find bin */ + + for(i=0; i bins[i]) && (f <= bins[i+1])) + b = i; + assert(b < MAX_BINS); + + /* est predicted phase from average */ + + pred = pexp->phi_prev[m] + N*m*best_Wo; + err = model->phi[m] - pred; + av[b].real += cos(err); av[b].imag += sin(err); + } + + } + + /* use averages to est phases */ + + for(m=1; m<=model->L; m++) { + f = m*model->Wo*FS/TWO_PI; + if (f > bins[0]) { + + /* find bin */ + + for(i=0; i bins[i]) && (f <= bins[i+1])) + b = i; + assert(b < MAX_BINS); + + /* add predicted phase error to this bin */ + + printf("L %d m %d f %4.f b %d\n", model->L, m, f, b); + + pred = pexp->phi_prev[m] + N*m*best_Wo; + err = atan2(av[b].imag, av[b].real); + printf(" %d: %4.3f -> ", m, model->phi[m]); + model->phi[m] = pred + err; + model->phi[m] = atan2(sin(model->phi[m]),cos(model->phi[m])); + printf("%4.3f\n", model->phi[m]); + } + } + printf("\n"); +} + + +/* + Try to code the phase of the largest amplitude in each band. Randomise the + phase of the other harmonics. The theory is that only the largest harmonic + will be audible. +*/ + +void cb_phase1(struct PEXP *pexp, MODEL *model) { + int m, i; + int nbins; + int b; + float f, best_Wo; + float max_val[MAX_BINS]; + int max_ind[MAX_BINS]; + + nbins = sizeof(bins)/sizeof(float); + best_Wo = refine_Wo(pexp, model, 1, model->L); + + for(i=0; iL; m++) { + f = m*model->Wo*FS/TWO_PI; + if (f > bins[0]) { + + /* find bin */ + + for(i=0; i bins[i]) && (f <= bins[i+1])) + b = i; + assert(b < MAX_BINS); + + if (model->A[m] > max_val[b]) { + max_val[b] = model->A[m]; + max_ind[b] = m; + } + } + + } + + /* randomise phase of other harmonics */ + + for(m=1; m<=model->L; m++) { + f = m*model->Wo*FS/TWO_PI; + if (f > bins[0]) { + + /* find bin */ + + for(i=0; i bins[i]) && (f <= bins[i+1])) + b = i; + assert(b < MAX_BINS); + + if (m != max_ind[b]) + model->phi[m] = pexp->phi_prev[m] + N*m*best_Wo; + } + } +} + + +/* + Theory is only the phase of the envelope of signal matters within a + Critical Band. So we estimate the position of an impulse that + approximates the envelope of the signal. +*/ + +void cb_phase2(struct PEXP *pexp, MODEL *model) { + int st, m, i, a, b, step; + float diff,w,c,s,phi1_; + float A[MAX_AMP]; + + for(m=1; m<=model->L; m++) { + A[m] = model->A[m]; + model->A[m] = 0; + } + + st = 2*model->L/4; + step = 3; + model->phi[1] = pexp->phi_prev[1] + (pexp->Wo_prev+model->Wo)*N/2.0; + + printf("L=%d ", model->L); + for(m=st; m model->L) + b = model->L; + + c = s = 0; + for(i=a; iphi[i+1] - model->phi[i]; + //w = (model->A[i+1] + model->A[i])/2; + w = 1.0; + c += w*cos(diff); s += w*sin(diff); + } + phi1_ = atan2(s,c); + printf("replacing: "); + for(i=a; iphi[i] = i*phi1_; + //model->phi[i] = i*model->phi[1]; + //model->phi[i] = m*(pexp->Wo_prev+model->Wo)*N/2.0; + model->A[m] = A[m]; + printf("%d ", i); + } + printf(" . "); + } + printf("\n"); +} + + +static void smooth_phase4(MODEL *model) { + int m; + float phi_m, phi_m_1; + + if (model->L > 25) { + printf("\nL %d\n", model->L); + for(m=model->L/2; m<=model->L; m+=2) { + if ((m+1) <= model->L) { + phi_m = (model->phi[m] - model->phi[m+1])/2.0; + phi_m_1 = (model->phi[m+1] - model->phi[m])/2.0; + model->phi[m] = phi_m; + model->phi[m+1] = phi_m_1; + printf("%d %4.3f %4.3f ", m, phi_m, phi_m_1); + } + } + } + +} + +/* try repeating last frame, just advance phases to account for time shift */ + +static void repeat_phases(struct PEXP *pexp, MODEL *model) { + int m; + + *model = pexp->prev_model; + for(m=1; m<=model->L; m++) + model->phi[m] += N*m*model->Wo; + +} + +/*---------------------------------------------------------------------------*\ + + phase_experiment() + + Phase quantisation experiments. + +\*---------------------------------------------------------------------------*/ + +void phase_experiment(struct PEXP *pexp, MODEL *model, char *arg) { + int m; + float before[MAX_AMP], best_Wo; + + assert(pexp != NULL); + memcpy(before, &model->phi[0], sizeof(float)*MAX_AMP); + + if (strcmp(arg,"q3") == 0) { + quant_phases(model, 1, model->L, 3); + update_snr_calc(pexp, model, before); + update_variance_calc(pexp, model, before); + } + + if (strcmp(arg,"dec2") == 0) { + if ((pexp->frames % 2) != 0) { + predict_phases(pexp, model, 1, model->L); + update_snr_calc(pexp, model, before); + update_variance_calc(pexp, model, before); + } + } + + if (strcmp(arg,"repeat") == 0) { + if ((pexp->frames % 2) != 0) { + repeat_phases(pexp, model); + update_snr_calc(pexp, model, before); + update_variance_calc(pexp, model, before); + } + } + + if (strcmp(arg,"vq") == 0) { + sparse_vq_pred_error(pexp, model); + update_snr_calc(pexp, model, before); + update_variance_calc(pexp, model, before); + } + + if (strcmp(arg,"pred") == 0) + predict_phases_state(pexp, model, 1, model->L); + + if (strcmp(arg,"pred1k") == 0) + predict_phases(pexp, model, 1, model->L/4); + + if (strcmp(arg,"smooth") == 0) { + smooth_phase(pexp, model,0); + update_snr_calc(pexp, model, before); + } + if (strcmp(arg,"smoothtrain") == 0) + smooth_phase(pexp, model,1); + if (strcmp(arg,"smoothvq") == 0) { + smooth_phase(pexp, model,2); + update_snr_calc(pexp, model, before); + } + + if (strcmp(arg,"smooth2") == 0) + smooth_phase2(pexp, model); + if (strcmp(arg,"smooth3") == 0) + smooth_phase3(pexp, model); + if (strcmp(arg,"smooth4") == 0) + smooth_phase4(model); + if (strcmp(arg,"vqsmooth3") == 0) { + sparse_vq_pred_error(pexp, model); + smooth_phase3(pexp, model); + } + + if (strcmp(arg,"cb1") == 0) { + cb_phase1(pexp, model); + update_snr_calc(pexp, model, before); + } + + if (strcmp(arg,"top") == 0) { + //top_amp(pexp, model, 1, model->L/4, 4, 1); + //top_amp(pexp, model, model->L/4, model->L/3, 4, 1); + //top_amp(pexp, model, model->L/3+1, model->L/2, 4, 1); + //top_amp(pexp, model, model->L/2, model->L, 6, 1); + //rand_phases(model, model->L/2, 3*model->L/4); + //struct_phases(pexp, model, model->L/2, 3*model->L/4); + //update_snr_calc(pexp, model, before); + } + + if (strcmp(arg,"pred23") == 0) { + predict_phases2(pexp, model, model->L/2, model->L); + update_snr_calc(pexp, model, before); + } + + if (strcmp(arg,"struct23") == 0) { + struct_phases(pexp, model, model->L/2, 3*model->L/4 ); + update_snr_calc(pexp, model, before); + } + + if (strcmp(arg,"addnoise") == 0) { + int m; + float max; + + max = 0; + for(m=1; m<=model->L; m++) + if (model->A[m] > max) + max = model->A[m]; + max = 20.0*log10(max); + for(m=1; m<=model->L; m++) + if (20.0*log10(model->A[m]) < (max-20)) { + model->phi[m] += (PI/4)*(1.0 -2.0*rand()/RAND_MAX); + //printf("m %d\n", m); + } + } + + /* normalise phases */ + + for(m=1; m<=model->L; m++) + model->phi[m] = atan2(sin(model->phi[m]), cos(model->phi[m])); + + /* update states */ + + //best_Wo = refine_Wo(pexp, model, model->L/2, model->L); + pexp->phi1 += N*model->Wo; + + for(m=1; m<=model->L; m++) + pexp->phi_prev[m] = model->phi[m]; + pexp->Wo_prev = model->Wo; + pexp->frames++; + pexp->prev_model = *model; +} + +#ifdef OLD_STUFF + //quant_phases(model, 1, model->L, 3); + //update_variance_calc(pexp, model, before); + //print_sparse_pred_error(pexp, model, 1, model->L, 40.0); + + //sparse_vq_pred_error(pexp, model); + + //quant_phases(model, model->L/4+1, model->L, 3); + + //predict_phases1(pexp, model, 1, model->L/4); + //quant_phases(model, model->L/4+1, model->L, 3); + + //quant_phases(model, 1, model->L/8, 3); + + //update_snr_calc(pexp, model, before); + //update_variance_calc(pexp, model, before); + + //fixed_bits_per_frame(pexp, model, 40); + //struct_phases(pexp, model, 1, model->L/4); + //rand_phases(model, 10, model->L); + //for(m=1; m<=model->L; m++) + // model->A[m] = 0.0; + //model->A[model->L/2] = 1000; + //repeat_phases(model, 20); + //predict_phases(pexp, model, 1, model->L/4); + //quant_phases(model, 1, 10, 3); + //quant_phases(model, 10, 20, 2); + //repeat_phases(model, 20); + //rand_phases(model, 3*model->L/4, model->L); + // print_phi1_pred_error(model, 1, model->L); + //predict_phases(pexp, model, 1, model->L/4); + //first_order_band(model, model->L/4, model->L/2); + //first_order_band(model, model->L/2, 3*model->L/4); + //if (fabs(model->Wo - pexp->Wo_prev)< 0.1*model->Wo) + + //print_pred_error(pexp, model, 1, model->L, 40.0); + //print_sparse_pred_error(pexp, model, 1, model->L, 40.0); + + //phi1_est = est_phi1(model, 1, model->L/4); + //print_phi1_pred_error(model, 1, model->L/4); + + //first_order_band(model, 1, model->L/4, phi1_est); + //sub_linear(model, 1, model->L/4, phi1_est); + + //top_amp(pexp, model, 1, model->L/4, 4); + //top_amp(pexp, model, model->L/4, model->L/2, 4); + + //first_order_band(model, 1, model->L/4, phi1_est); + //first_order_band(model, model->L/4, model->L/2, phi1_est); + + //if (fabs(model->Wo - pexp->Wo_prev) > 0.2*model->Wo) + // rand_phases(model, model->L/2, model->L); + + //top_amp(pexp, model, 1, model->L/4, 4); + //top_amp(pexp, model, model->L/4, model->L/2, 8); + //top_amp(pexp, model, model->L/4+1, model->L/2, 10, 1); + //top_amp(pexp, model, 1, model->L/4, 10, 1); + //top_amp(pexp, model, model->L/4+1, 3*model->L/4, 10, 1); + //top_amp(pexp, model, 1, 3*model->L/4, 20, 1); + + #ifdef REAS_CAND1 + predict_phases(pexp, model, 1, model->L/4); + top_amp(pexp, model, model->L/4+1, 3*model->L/4, 10, 1); + rand_phases(model, 3*model->L/4+1, model->L); + #endif + + #ifdef REAS_CAND2 + if ((pexp->frames % 2) == 0) { + //printf("quant\n"); + predict_phases(pexp, model, 1, model->L/4); + //top_amp(pexp, model, model->L/4+1, 3*model->L/4, 20, 1); + top_amp(pexp, model, model->L/4+1, 7*model->L/8, 20, 1); + rand_phases(model, 7*model->L/8+1, model->L); + } + else { + //printf("predict\n"); + predict_phases(pexp, model, 1, model->L); + } + #endif + + //#define REAS_CAND3 + #ifdef REAS_CAND3 + if ((pexp->frames % 3) != 0) { + printf("pred\n"); + predict_phases(pexp, model, 1, model->L); + } + else { + predict_phases(pexp, model, 1, model->L/4); + fixed_bits_per_frame(pexp, model, model->L/4+1, 60); + } + #endif + //predict_phases(pexp, model, model->L/4, model->L); + + + //print_pred_error(pexp, model, 1, model->L); + //limit_prediction_error(pexp, model, model->L/2, model->L, PI/2); +#endif diff --git a/codec2-dev/src/phaseexp.h b/codec2-dev/src/phaseexp.h new file mode 100644 index 00000000..b43db75e --- /dev/null +++ b/codec2-dev/src/phaseexp.h @@ -0,0 +1,39 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: phaseexp.h + AUTHOR......: David Rowe + DATE CREATED: June 2012 + + Experimental functions for quantising, modelling and synthesising phase. + +\*---------------------------------------------------------------------------*/ + +/* + Copyright (C) 2012 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 . +*/ + +#ifndef __PHASEEXP__ +#define __PHASEEXP__ + +#include "kiss_fft.h" + +struct PEXP; + +struct PEXP * phase_experiment_create(); +void phase_experiment_destroy(struct PEXP *pexp); +void phase_experiment(struct PEXP *pexp, MODEL *model, char *arg); + +#endif diff --git a/codec2-dev/src/quantise.c b/codec2-dev/src/quantise.c index bc9545a1..0c3fbf39 100644 --- a/codec2-dev/src/quantise.c +++ b/codec2-dev/src/quantise.c @@ -518,54 +518,183 @@ void lspjvm_quantise(float *x, float *xq, int ndim) } } -/* 3 stage VQ LSP quantsier kindly submitted by Anssi, OH3GDD */ +#define MBEST_STAGES 3 -void lspanssi_quantise(float *x, float *xq, int ndim) +struct MBEST_LIST { + int index[MBEST_STAGES]; /* index of each stage that lead us to this error */ + float error; +}; + +struct MBEST { + int entries; /* number of entries in mbest list */ + struct MBEST_LIST *list; +}; + + +static struct MBEST *mbest_create(int entries) { + int i,j; + struct MBEST *mbest; + + assert(entries > 0); + mbest = (struct MBEST *)malloc(sizeof(struct MBEST)); + assert(mbest != NULL); + + mbest->entries = entries; + mbest->list = (struct MBEST_LIST *)malloc(entries*sizeof(struct MBEST_LIST)); + assert(mbest->list != NULL); + + for(i=0; ientries; i++) { + for(j=0; jlist[i].index[j] = 0; + mbest->list[i].error = 1E32; + } + + return mbest; +} + + +static void mbest_destroy(struct MBEST *mbest) { + assert(mbest != NULL); + free(mbest->list); + free(mbest); +} + + +/*---------------------------------------------------------------------------*\ + + mbest_insert + + Insert the results of a vector to codebook entry comparison. the + list is in order of ascending error, so those codebook entries with the + smallest error will be first on the list. + +\*---------------------------------------------------------------------------*/ + +static void mbest_insert(struct MBEST *mbest, int index[], float error) { + int i, j, found; + struct MBEST_LIST *list = mbest->list; + int entries = mbest->entries; + + found = 0; + for(i=0; ii; j--) + list[j] = list[j-1]; + for(j=0; jentries; i++) { + for(j=0; jlist[i].index[j]); + printf(" %f\n", mbest->list[i].error); + } +} + + +/*---------------------------------------------------------------------------*\ + + mbest_search + + Searches vec[] to a codebbook of vectors, and maintains a list of the mbest + closest matches. + +\*---------------------------------------------------------------------------*/ + +static void mbest_search( + const float *cb, /* VQ codebook to search */ + float vec[], /* target vector */ + float w[], /* weighting vector */ + int k, /* dimension of vector */ + int m, /* number on entries in codebook */ + struct MBEST *mbest, /* list of closest matches */ + int index[] /* indexes that lead us here */ +) { - int i, n1, n2, n3; - float e, err1[LPC_ORD], err2[LPC_ORD], err3[LPC_ORD]; + float e; + int i,j; + float diff; + + for(j=0; jlist[j].index[0]; + for(i=0; ilist[j].index[1]; + index[1] = n2 = mbest_stage2->list[j].index[0]; + for(i=0; ilist[0].index[2]; + n2 = mbest_stage3->list[0].index[1]; + n3 = mbest_stage3->list[0].index[0]; for (i=0;i