--- /dev/null
+% 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
${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
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
#
codebookjvm.c
codebookmel.c
codebooklspmelvq.c
+ codebooknewamp1.c
codebookge.c
- golay23.c\r
- freedv_api.c\r
- freedv_vhf_framing.c\r
- freedv_data_channel.c\r
- varicode.c\r
- modem_stats.c\r
-)\r
+ golay23.c
+ freedv_api.c
+ freedv_vhf_framing.c
+ freedv_data_channel.c
+ varicode.c
+ modem_stats.c
+)
set(CODEC2_PUBLIC_HEADERS
golay23.h
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
/* 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];
}
for(i=0; i<FFT_ENC/2; i++)
- fprintf(fpwb,"%f\t",Pwb[i].real);
+ fprintf(fpwb,"%f\t",Pwb[i]);
fprintf(fpwb,"\n");
}
-void dump_Pw(COMP Pw[]) {
+void dump_Pw(float Pw[]) {
int i;
char s[MAX_STR];
}
for(i=0; i<FFT_ENC/2; i++)
- fprintf(fpw,"%f\t",Pw[i].real);
+ fprintf(fpw,"%f\t",Pw[i]);
fprintf(fpw,"\n");
}
void dump_model(MODEL *m);
void dump_quantised_model(MODEL *m);
void dump_Pwn(COMP Pw[]);
-void dump_Pw(COMP Pw[]);
+void dump_Pw(float Pw[]);
void dump_Rw(float Rw[]);
void dump_lsp(float lsp[]);
void dump_weights(float w[], int ndim);
/* post filter */
void dump_bg(float e, float bg_est, float percent_uv);
-void dump_Pwb(COMP Pwb[]);
+void dump_Pwb(float Pwb[]);
#endif
*e = powf(10.0, xq[1]/10.0);
}
+/*---------------------------------------------------------------------------*\
+
+ FUNCTION....: interp_para()
+ AUTHOR......: David Rowe
+ DATE CREATED: Jan 2017
+
+ General 2nd order parabolic interpolator. Used splines orginally,
+ but this is much simpler and we don't need much accuracy. Given two
+ vectors of points xp and yp, find interpolated values y at points x.
+
+\*---------------------------------------------------------------------------*/
+
+void interp_para(float y[], float xp[], float yp[], int np, float x[], int n)
+{
+ assert(np >= 3);
+
+ int k,i;
+ float xi, x1, y1, x2, y2, x3, y3, a, b;
+
+ k = 0;
+ for (i=0; i<n; i++) {
+ xi = x[i];
+
+ /* k is index into xp of where we start 3 points used to form parabola */
+
+ while ((xp[k+1] < xi) && (k < (np-3)))
+ k++;
+
+ x1 = xp[k]; y1 = yp[k]; x2 = xp[k+1]; y2 = yp[k+1]; x3 = xp[k+2]; y3 = yp[k+2];
+
+ printf("k: %d np: %d i: %d xi: %f x1: %f y1: %f\n", k, np, i, xi, x1, y1);
+
+ a = ((y3-y2)/(x3-x2)-(y2-y1)/(x2-x1))/(x3-x1);
+ b = ((y3-y2)/(x3-x2)*(x2-x1)+(y2-y1)/(x2-x1)*(x3-x2))/(x3-x1);
+
+ y[i] = a*(xi-x2)*(xi-x2) + b*(xi-x2) + y2;
+ }
+}
+
+
+/*---------------------------------------------------------------------------*\
+
+ FUNCTION....: ftomel()
+ AUTHOR......: David Rowe
+ DATE CREATED: Jan 2017
+
+ Non linear sampling of frequency axis, reducing the "rate" is a
+ first step before VQ
+
+\*---------------------------------------------------------------------------*/
+
+float ftomel(float fHz) {
+ float mel = floorf(2595.0*log10f(1.0 + fHz/700.0)+0.5);
+ return mel;
+}
+
+void mel_sample_freqs_kHz(float rate_K_sample_freqs_kHz[], int K)
+{
+ float mel_start = ftomel(200.0); float mel_end = ftomel(3700.0);
+ float step = (mel_end-mel_start)/(K-1);
+ float mel;
+ int k;
+
+ mel = mel_start;
+ for (k=0; k<K; k++) {
+ rate_K_sample_freqs_kHz[k] = 0.7*(pow(10.0, (mel/2595.0)) - 1.0);
+ mel += step;
+ }
+}
+
+
+/*---------------------------------------------------------------------------*\
+
+ FUNCTION....: resample_const_rate_f()
+ AUTHOR......: David Rowe
+ DATE CREATED: Jan 2017
+
+ Resample Am from time-varying rate L=floor(pi/Wo) to fixed rate K.
+
+\*---------------------------------------------------------------------------*/
+
+void resample_const_rate_f(MODEL *model, float rate_K_vec[], float rate_K_sample_freqs_kHz[], int K)
+{
+ int m;
+ float AmdB[MAX_AMP+1], rate_L_sample_freqs_kHz[MAX_AMP+1], AmdB_peak;
+
+ /* convert rate L=pi/Wo amplitude samples to fixed rate K */
+
+ AmdB_peak = -100.0;
+ for(m=1; m<=model->L; 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);
+}
+
+
+
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
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)
--- /dev/null
+/*---------------------------------------------------------------------------*\
+
+ 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 <http://www.gnu.org/licenses/>.
+*/
+
+#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<M_PITCH; i++) {
+ Sn[i] = 1.0;
+ }
+
+ int K = 20;
+ float rate_K_sample_freqs_kHz[K];
+ float model_octave[FRAMES][MAX_AMP+2]; // model params in matrix format, useful for C <-> Octave
+ float rate_K_surface[FRAMES][K]; // rate K vecs for each frame, form a surface that makes pretty graphs
+
+ for(f=0; f<FRAMES; f++)
+ for(m=0; m<MAX_AMP+2; m++)
+ model_octave[f][m] = 0.0;
+
+ mel_sample_freqs_kHz(rate_K_sample_freqs_kHz, K);
+ //for(int k=0; k<K; k++)
+ // printf("k: %d sf: %f\n", k, rate_K_sample_freqs_kHz[k]);
+
+ FILE *fin = fopen(argv[1], "rb");
+ if (fin == NULL) {
+ fprintf(stderr, "Problem opening hts1.raw\n");
+ exit(1);
+ }
+
+ for(f=0; f<FRAMES; f++) {
+ assert(fread(buf,sizeof(short),N_SAMP,fin) == N_SAMP);
+
+ /* shift buffer of input samples, and insert new samples */
+
+ for(i=0; i<M_PITCH-N_SAMP; i++) {
+ Sn[i] = Sn[i+N_SAMP];
+ }
+ for(i=0; i<N_SAMP; i++) {
+ Sn[i+M_PITCH-N_SAMP] = buf[i];
+ }
+
+ /* Estimate Sinusoidal Model Parameters ----------------------*/
+
+ nlp(nlp_states, Sn, N_SAMP, P_MIN, P_MAX, &pitch, Sw, W, &prev_uq_Wo);
+ model.Wo = TWO_PI/pitch;
+
+ dft_speech(fft_fwd_cfg, Sw, Sn, w);
+ two_stage_pitch_refinement(&model, Sw);
+ estimate_amplitudes(&model, Sw, W, 1);
+
+ /* Resample at rate K ----------------------------------------*/
+
+ resample_const_rate_f(&model, &rate_K_surface[f][0], rate_K_sample_freqs_kHz, K);
+ for(int k=0; k<K; k++)
+ printf("k: %d sf: %f sv: %f\n", k, rate_K_sample_freqs_kHz[k], rate_K_surface[f][k]);
+ printf("\n");
+
+ /* log vectors */
+
+ model_octave[f][0] = model.Wo;
+ model_octave[f][1] = model.L;
+ for(m=1; m<=model.L; m++) {
+ model_octave[f][m+1] = model.A[m];
+ }
+ }
+
+ fclose(fin);
+
+ /* save vectors in Octave format */
+
+ FILE *fout = fopen("tnewamp1_out.txt","wt");
+ assert(fout != NULL);
+ fprintf(fout, "# Created by tnewamp1.c\n");
+ octave_save_float(fout, "rate_K_surface_c", (float*)rate_K_surface, FRAMES, K, K);
+ octave_save_float(fout, "model_c", (float*)model_octave, FRAMES, MAX_AMP+2, MAX_AMP+2);
+ fclose(fout);
+
+ printf("Done! Now run\n octave:1> tnewamp1(\"../build_linux/src/hts1a\")\n");
+ return 0;
+}
+
projects / freetel-svn-tracking.git / commitdiff