From: drowe67 Date: Thu, 5 Jan 2017 02:51:06 +0000 (+0000) Subject: first stage C<->Octave test framework for newamp1 passing X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=d38896a4428461c0487595b6392f4f59aca1873f;p=freetel-svn-tracking.git first stage C<->Octave test framework for newamp1 passing git-svn-id: https://svn.code.sf.net/p/freetel/code@2951 01035d8c-6547-0410-b346-abe4f91aad63 --- diff --git a/codec2-dev/octave/tnewamp1.m b/codec2-dev/octave/tnewamp1.m new file mode 100644 index 00000000..59eb4304 --- /dev/null +++ b/codec2-dev/octave/tnewamp1.m @@ -0,0 +1,272 @@ +% tnewamp1.m +% +% Copyright David Rowe 2017 +% This program is distributed under the terms of the GNU General Public License +% Version 2 + +#{ + + Octave script to compare Octave and C versions of newamp1 processing, in order to test C port. + + c2sim -> dump files -> $ ../build_linux/unittest/tnewamp1 -> octave:1> tnewamp1 + Usage: + + 1/ build codec2 with -DDUMP - see codec2-dev/README + + 2/ Generate dump files using c2sim (just need to do this once) + $ cd codec2-dev/build_linux/src + $ ./c2sim ../../raw/hts1a.raw --phase0 --postfilter --dump hts1a --lpc 10 --dump_pitch_e hts1a_pitche.txt + + 3/ Run C version which generates a file of Octave test vectors as ouput: + + $ cd codec2-dev/build_linux/unittest + $ ./tnewamp1 ../../raw/hts1a + + 4/ Run Octave script to generate Octave test vectors and compare with C. + + octave:1> tnewamp1("../build_linux/src/hts1a") + +#} + + +% In general, this function processes a bunch of amplitudes, we then +% use c2sim to hear the results. Bunch of different experiments below + +function tnewamp1(input_prefix) + newamp; + more off; + + max_amp = 80; + postfilter = 0; % optional postfiler that runs on Am, not used atm + synth_phase = 1; + + if nargin == 1 + output_prefix = input_prefix; + end + model_name = strcat(input_prefix,"_model.txt"); + model = load(model_name); + [frames nc] = size(model); + + voicing_name = strcat(input_prefix,"_pitche.txt"); + voicing = zeros(1,frames); + + if exist(voicing_name, "file") == 2 + pitche = load(voicing_name); + voicing = pitche(:, 3); + end + + % Load in C vectors and compare ----------------------------------------- + + load("../build_linux/unittest/tnewamp1_out.txt"); + K = 20; + [frames tmp] = size(rate_K_surface_c); + [rate_K_surface sample_freqs_kHz] = resample_const_rate_f_mel(model(1:frames,:), K); +rate_K_surface + figure(1); + mesh(rate_K_surface); + figure(2); + mesh(rate_K_surface_c); + figure(3); + mesh(rate_K_surface - rate_K_surface_c); + + #{ + + for f=1:frames + printf("%d ", f); + Wo = model_(f,1); + L = min([model_(f,2) max_amp-1]); + Am = model_(f,3:(L+2)); + + Am_ = zeros(1,max_amp); + Am_(2:L) = Am(1:L-1); + + % optional post filter on linear {Am}, boosts higher amplitudes more than lower, + % improving shape of formants and reducing muffling. Note energy + % normalisation + + if postfilter + e1 = sum(Am_(2:L).^2); + Am_(2:L) = Am_(2:L) .^ 1.5; + e2 = sum(Am_(2:L).^2); + Am_(2:L) *= sqrt(e1/e2); + end + + fwrite(fam, Am_, "float32"); + fwrite(fWo, Wo, "float32"); + + if synth_phase + + % synthesis phase spectra from magnitiude spectra using minimum phase techniques + + fft_enc = 128; + phase = determine_phase(model_, f, fft_enc); + assert(length(phase) == fft_enc); + %Aw = zeros(1, fft_enc*2); + %Aw(1:2:fft_enc*2) = cos(phase); + %Aw(2:2:fft_enc*2) = -sin(phase); + + % sample phase at centre of each harmonic, not 1st entry Hm[1] in octave Hm[0] in C + % is not used + + Hm = zeros(1, 2*max_amp); + for m=1:L + b = round(m*Wo*fft_enc/(2*pi)); + Hm(2*m) = cos(phase(b)); + Hm(2*m+1) = -sin(phase(b)); + end + fwrite(fhm, Hm, "float32"); + end + end + + #} + + printf("\n") + +endfunction + + +% ----------------------------------------------------------------------------------------- +% Linear decimator/interpolator that operates at rate K, includes VQ, post filter, and Wo/E +% quantisation. Evolved from abys decimator below. Simulates the entire encoder/decoder. + +function [model_ voicing_ indexes] = experiment_rate_K_dec(model, voicing) + max_amp = 80; + [frames nc] = size(model); + model_ = zeros(frames, max_amp+3); + indexes = zeros(frames,4); + + M = 4; + + % create frames x K surface. TODO make all of this operate frame by + % frame, or at least M/2=4 frames rather than one big chunk + + K = 20; + [surface sample_freqs_kHz] = resample_const_rate_f_mel(model, K); + target_surface = surface; + + figure(1); + mesh(surface); + + % VQ rate K surface. TODO: If we are decimating by M/2=4 we really + % only need to do this every 4th frame. + + melvq; + load train_120_vq; m=5; + + for f=1:frames + mean_f(f) = mean(surface(f,:)); + surface_no_mean(f,:) = surface(f,:) - mean_f(f); + end + figure(2); + mesh(surface_no_mean); + + [res surface_no_mean_ ind] = mbest(train_120_vq, surface_no_mean, m); + indexes(:,1:2) = ind; + + for f=1:frames + surface_no_mean_(f,:) = post_filter(surface_no_mean_(f,:), sample_freqs_kHz, 1.5); + end + figure(3); + mesh(surface_no_mean_); + + surface_ = zeros(frames, K); + energy_q = 10 + 40/16*(0:15); + for f=1:frames + [mean_f_ indx] = quantise(energy_q, mean_f(f)); + indexes(f,3) = indx - 1; + %mean_f_ = mean_f(f); + surface_(f,:) = surface_no_mean_(f,:) + mean_f_; + end + + figure(); + mesh(surface_); + + % break into segments of M frames. We have 3 samples in M frame + % segment spaced M/2 apart and interpolate the rest. This evolved + % from AbyS scheme below but could be simplified to simple linear + % interpolation, or using 3 or 4 points but shift of M/2=4 frames. + + interpolated_surface_ = zeros(frames, K); + for f=1:M:frames-M + left_vec = surface_(f,:); + right_vec = surface_(f+M,:); + sample_points = [f f+M]; + resample_points = f:f+M-1; + for k=1:K + interpolated_surface_(resample_points,k) = interp_linear(sample_points, [left_vec(k) right_vec(k)], resample_points); + end + end + + % break into segments for purposes of Wo interpolation + + for f=1:M:frames + % quantise Wo + + % UV/V flag is coded using a zero index for Wo, this means we need to + % adjust Wo index slightly for the lowest Wo V frames + + if voicing(f) + index = encode_log_Wo(model(f,1), 6); + if index == 0 + index = 1; + end + indexes(f,4) = index; + voicing(f) = 1; + model_(f,1) = decode_log_Wo(indexes(f,4), 6); + else + indexes(f,4) = 0; + voicing(f) = 0; + model_(f,1) = 2*pi/100; + end + end + + + voicing_ = zeros(1, frames); + for f=1:M:frames-M + + Wo1_ = model_(f,1); + Wo2_ = model_(f+M,1); + + % uncomment to use unquantised values + %Wo1_ = model(f,1); + %Wo2_ = model(f+M,1); + + if !voicing(f) && !voicing(f+M) + model_(f:f+M-1,1) = 2*pi/100; + end + + if voicing(f) && !voicing(f+M) + model_(f:f+M/2-1,1) = Wo1_; + model_(f+M/2:f+M-1,1) = 2*pi/100; + voicing_(f:f+M/2-1) = 1; + end + + if !voicing(f) && voicing(f+M) + model_(f:f+M/2-1,1) = 2*pi/100; + model_(f+M/2:f+M-1,1) = Wo2_; + voicing_(f+M/2:f+M-1) = 1; + end + + if voicing(f) && voicing(f+M) + Wo_samples = [Wo1_ Wo2_]; + model_(f:f+M-1,1) = interp1([f f+M], Wo_samples, f:f+M-1, "linear", 0); + voicing_(f:f+M-1) = 1; + end + + #{ + printf("f: %d f+M/2: %d Wo: %f %f (%f %%) v: %d %d \n", f, f+M/2, model(f,1), model(f+M/2,1), 100*abs(model(f,1) - model(f+M/2,1))/model(f,1), voicing(f), voicing(f+M/2)); + for i=f:f+M/2-1 + printf(" f: %d v: %d v_: %d Wo: %f Wo_: %f\n", i, voicing(i), voicing_(i), model(i,1), model_(i,1)); + end + #} + end + model_(frames-M:frames,1) = pi/100; % set end frames to something sensible + + % enable these to use original (non interpolated) voicing and Wo + %voicing_ = voicing; + %model_(:,1) = model(:,1); + + model_(:,2) = floor(pi ./ model_(:,1)); % calculate L for each interpolated Wo + model_ = resample_rate_L(model_, interpolated_surface_, sample_freqs_kHz); + +endfunction diff --git a/codec2-dev/src/CMakeLists.txt b/codec2-dev/src/CMakeLists.txt index 658a8dbf..5abc3e4f 100644 --- a/codec2-dev/src/CMakeLists.txt +++ b/codec2-dev/src/CMakeLists.txt @@ -55,6 +55,11 @@ set(CODEBOOKSLSPMELVQ ${D}/lspmelvq3.txt ) +set(CODEBOOKSNEWAMP1 + ${D}/train_120_1.txt + ${D}/train_120_2.txt +) + set(CODEBOOKSGE ${D}/gecb.txt) # when crosscompiling we need a native executable @@ -126,6 +131,13 @@ add_custom_command( DEPENDS generate_codebook ${CODEBOOKSGE} ) +# codebooknewamp1.c +add_custom_command( + OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/codebooknewamp1.c + COMMAND generate_codebook newamp1vq_cb ${CODEBOOKSNEWAMP1} > ${CMAKE_CURRENT_BINARY_DIR}/codebooknewamp1.c + DEPENDS generate_codebook ${CODEBOOKSNEWAMP1} +) + # # codec2 library sources # @@ -157,14 +169,15 @@ set(CODEC2_SRCS codebookjvm.c codebookmel.c codebooklspmelvq.c + codebooknewamp1.c codebookge.c - golay23.c - freedv_api.c - freedv_vhf_framing.c - freedv_data_channel.c - varicode.c - modem_stats.c -) + golay23.c + freedv_api.c + freedv_vhf_framing.c + freedv_data_channel.c + varicode.c + modem_stats.c +) set(CODEC2_PUBLIC_HEADERS golay23.h diff --git a/codec2-dev/src/defines.h b/codec2-dev/src/defines.h index 9ec7235d..c1ebcdd9 100644 --- a/codec2-dev/src/defines.h +++ b/codec2-dev/src/defines.h @@ -94,5 +94,6 @@ extern const struct lsp_codebook lsp_cbvqanssi[]; extern const struct lsp_codebook mel_cb[]; extern const struct lsp_codebook ge_cb[]; extern const struct lsp_codebook lspmelvq_cb[]; +extern const struct lsp_codebook newamp1vq_cb[]; #endif diff --git a/codec2-dev/src/dump.c b/codec2-dev/src/dump.c index 9c29bfb7..c993d83b 100644 --- a/codec2-dev/src/dump.c +++ b/codec2-dev/src/dump.c @@ -378,7 +378,7 @@ void dump_lpc_snr(float snr) { /* Pw "before" post filter so we can plot before and after */ -void dump_Pwb(COMP Pwb[]) { +void dump_Pwb(float Pwb[]) { int i; char s[MAX_STR]; @@ -391,11 +391,11 @@ void dump_Pwb(COMP Pwb[]) { } for(i=0; i= 3); + + int k,i; + float xi, x1, y1, x2, y2, x3, y3, a, b; + + k = 0; + for (i=0; iL; m++) { + AmdB[m] = 20.0*log10(model->A[m]+1E-16); + if (AmdB[m] > AmdB_peak) { + AmdB_peak = AmdB[m]; + } + rate_L_sample_freqs_kHz[m] = m*model->Wo*4.0/M_PI; + printf("m: %d AmdB: %f AmdB_peak: %f sf: %f\n", m, AmdB[m], AmdB_peak, rate_L_sample_freqs_kHz[m]); + } + + /* clip between peak and peak -50dB, to reduce dynamic range */ + + for(m=1; m<=model->L; m++) { + if (AmdB[m] < (AmdB_peak-50.0)) { + AmdB[m] = AmdB_peak-50.0; + } + } + + interp_para(rate_K_vec, &rate_L_sample_freqs_kHz[1], &AmdB[1], model->L, rate_K_sample_freqs_kHz, K); +} + + + diff --git a/codec2-dev/src/quantise.h b/codec2-dev/src/quantise.h index bc8c056a..e6df842c 100644 --- a/codec2-dev/src/quantise.h +++ b/codec2-dev/src/quantise.h @@ -135,4 +135,9 @@ float decode_amplitudes(MODEL *model, float lsps[], float *e); +void interp_para(float y[], float xp[], float yp[], int np, float x[], int n); +float ftomel(float fHz); +void mel_sample_freqs_kHz(float rate_K_sample_freqs_kHz[], int K); +void resample_const_rate_f(MODEL *model, float rate_K_vec[], float rate_K_sample_freqs_kHz[], int K); + #endif diff --git a/codec2-dev/unittest/CMakeLists.txt b/codec2-dev/unittest/CMakeLists.txt index a3713d56..5c328f9d 100644 --- a/codec2-dev/unittest/CMakeLists.txt +++ b/codec2-dev/unittest/CMakeLists.txt @@ -113,3 +113,6 @@ target_link_libraries(raw2h codec2) add_definitions(-D__UNITTEST__) add_executable(c2validate c2validate.c) target_link_libraries(c2validate codec2) + +add_executable(tnewamp1 tnewamp1.c ../src/quantise.c ../src/kiss_fft.c ../src/sine.c ../src/nlp.c ../src/dump.c ../src/octave.c ${CODEBOOKS}) +target_link_libraries(tnewamp1 codec2) diff --git a/codec2-dev/unittest/tnewamp1.c b/codec2-dev/unittest/tnewamp1.c new file mode 100644 index 00000000..586ba9a0 --- /dev/null +++ b/codec2-dev/unittest/tnewamp1.c @@ -0,0 +1,136 @@ +/*---------------------------------------------------------------------------*\ + + FILE........: tnewamp1.c + AUTHOR......: David Rowe + DATE CREATED: Jan 2017 + + Tests for the C version of the newamp1 amplitude modelling used for + 700c. This program outputs a file of Octave vectors that are loaded + and automatically tested against the Octave version of the modem by + the Octave script tnewamp1.m + +\*---------------------------------------------------------------------------*/ + +/* + Copyright (C) 2017 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, 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 "codec2_fft.h" +#include "sine.h" +#include "nlp.h" +#include "dump.h" +#include "octave.h" +#include "quantise.h" + +#define FRAMES 100 + +int main(int argc, char *argv[]) { + short buf[N_SAMP]; /* input/output buffer */ + float Sn[M_PITCH]; /* float input speech samples */ + COMP Sw[FFT_ENC]; /* DFT of Sn[] */ + codec2_fft_cfg fft_fwd_cfg; /* fwd FFT states */ + float w[M_PITCH]; /* time domain hamming window */ + COMP W[FFT_ENC]; /* DFT of w[] */ + MODEL model; + void *nlp_states; + float pitch, prev_uq_Wo; + int i,m,f; + + if (argc != 2) { + printf("usage: ./tnewamp1 RawFile\n"); + exit(1); + } + nlp_states = nlp_create(M_PITCH); + prev_uq_Wo = TWO_PI/P_MAX; + fft_fwd_cfg = codec2_fft_alloc(FFT_ENC, 0, NULL, NULL); + make_analysis_window(fft_fwd_cfg, w, W); + + for(i=0; i Octave + float rate_K_surface[FRAMES][K]; // rate K vecs for each frame, form a surface that makes pretty graphs + + for(f=0; f tnewamp1(\"../build_linux/src/hts1a\")\n"); + return 0; +} +