--- /dev/null
+/*---------------------------------------------------------------------------*\
+
+ FILE........: codec2_internal.h
+ AUTHOR......: David Rowe
+ DATE CREATED: 22 March 2011
+
+ Some internal structures and states broken out here as they are useful for
+ testing and development.
+
+\*---------------------------------------------------------------------------*/
+
+/*
+ Copyright (C) 2011 David Rowe
+
+ All rights reserved.
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License version 2.1, as
+ published by the Free Software Foundation. This program is
+ distributed in the hope that it will be useful, but WITHOUT ANY
+ WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
+ License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef __CODEC2_INTERNAL__
+#define __CODEC2_INTERNAL__
+
+/*---------------------------------------------------------------------------*\
+
+ STATES
+
+\*---------------------------------------------------------------------------*/
+
+typedef struct {
+ float w[M]; /* time domain hamming window */
+ COMP W[FFT_ENC]; /* DFT of w[] */
+ float Pn[2*N]; /* trapezoidal synthesis window */
+ float Sn[M]; /* input speech */
+ float hpf_states[2]; /* high pass filter states */
+ void *nlp; /* pitch predictor states */
+ float Sn_[2*N]; /* synthesised output speech */
+ float ex_phase; /* excitation model phase track */
+ float bg_est; /* background noise estimate for post filter */
+ float prev_Wo; /* previous frame's pitch estimate */
+ MODEL prev_model; /* previous frame's model parameters */
+ float prev_lsps[LPC_ORD]; /* previous frame's LSPs */
+ float prev_energy; /* previous frame's LPC energy */
+} CODEC2;
+
+/*---------------------------------------------------------------------------*\
+
+ FUNCTION HEADERS
+
+\*---------------------------------------------------------------------------*/
+
+void analyse_one_frame(CODEC2 *c2, MODEL *model, short speech[]);
+void synthesise_one_frame(CODEC2 *c2, short speech[], MODEL *model,float ak[]);
+
+#endif
#include <assert.h>
#include <math.h>
#include <string.h>
+#include <stdio.h>
+#include <stdlib.h>
#include "defines.h"
#include "interp.h"
#include "lsp.h"
#include "quantise.h"
+#include "dump.h"
float sample_log_amp(MODEL *model, float w);
assert(w > 0.0); assert (w <= PI);
- m = floor(w/model->Wo + 0.5);
- f = (w - m*model->Wo)/w;
+ m = 0;
+ while ((m+1)*model->Wo < w) m++;
+ f = (w - m*model->Wo)/model->Wo;
assert(f <= 1.0);
if (m < 1) {
else {
log_amp = (1.0-f)*log10(model->A[m] + 1E-6) +
f*log10(model->A[m+1] + 1E-6);
+ //printf("m=%d A[m] %f A[m+1] %f x %f %f %f\n", m, model->A[m],
+ // model->A[m+1], pow(10.0, log_amp),
+ // (1-f), f);
}
return log_amp;
}
+/*---------------------------------------------------------------------------*\
+
+ FUNCTION....: sample_log_amp_quad()
+ AUTHOR......: David Rowe
+ DATE CREATED: 9 March 2011
+
+ Samples the amplitude envelope at an arbitrary frequency w. Uses
+ quadratic interpolation in the log domain to sample between harmonic
+ amplitudes.
+
+ y(x) = ax*x + bx + c
+
+ We assume three points are x=-1, x=0, x=1, which we map to m-1,m,m+1
+
+ c = y(0)
+ b = (y(1) - y(-1))/2
+ a = y(-1) + b - y(0)
+
+\*---------------------------------------------------------------------------*/
+
+float sample_log_amp_quad(MODEL *model, float w)
+{
+ int m;
+ float a,b,c,x, log_amp;
+
+ assert(w > 0.0); assert (w <= PI);
+
+ m = floor(w/model->Wo + 0.5);
+ if (m < 2) m = 2;
+ if (m > (model->L-1)) m = model->L-1;
+ c = log10(model->A[m]+1E-6);
+ b = (log10(model->A[m+1]+1E-6) - log10(model->A[m-1]+1E-6))/2.0;
+ a = log10(model->A[m-1]+1E-6) + b - c;
+ x = (w - m*model->Wo)/model->Wo;
+
+ log_amp = a*x*x + b*x + c;
+ //printf("m=%d A[m-1] %f A[m] %f A[m+1] %f w %f x %f log_amp %f\n", m,
+ // model->A[m-1],
+ // model->A[m], model->A[m+1], w, x, pow(10.0, log_amp));
+ return log_amp;
+}
+
+/*---------------------------------------------------------------------------*\
+
+ FUNCTION....: sample_log_amp_quad_nl()
+ AUTHOR......: David Rowe
+ DATE CREATED: 10 March 2011
+
+ Samples the amplitude envelope at an arbitrary frequency w. Uses
+ quadratic interpolation in the log domain to sample between harmonic
+ amplitudes. This version can handle non-linear steps along a freq
+ axis defined by arbitrary steps.
+
+ y(x) = ax*x + bx + c
+
+ We assume three points are (x_1,y_1), (0,y0) and (x1,y1).
+
+\*---------------------------------------------------------------------------*/
+
+float sample_log_amp_quad_nl(
+ float w[], /* frequency points */
+ float A[], /* for these amplitude samples */
+ int np, /* number of frequency points */
+ float w_sample /* frequency of new samples */
+)
+{
+ int m,i;
+ float a,b,c,x, log_amp, best_dist;
+ float x_1, x1;
+ float y_1, y0, y1;
+
+ //printf("w_sample %f\n", w_sample);
+ assert(w_sample >= 0.0); assert (w_sample <= 1.1*PI);
+
+ /* find closest point to centre quadratic interpolator */
+
+ best_dist = 1E32;
+ for (i=0; i<np; i++)
+ if (fabs(w[i] - w_sample) < best_dist) {
+ best_dist = fabs(w[i] - w_sample);
+ m = i;
+ }
+
+ /* stay one point away from edge of array */
+
+ if (m < 1) m = 1;
+ if (m > (np-2)) m = np - 2;
+
+ /* find polynomial coeffs */
+
+ x_1 = w[m-1]- w[m]; x1 = w[m+1] - w[m];
+ y_1 = log10(A[m-1]+1E-6);
+ y0 = log10(A[m]+1E-6);
+ y1 = log10(A[m+1]+1E-6);
+
+ c = y0;
+ a = (y_1*x1 - y1*x_1 + c*x_1 - c*x1)/(x_1*x_1*x1 - x1*x1*x_1);
+ b = (y1 -a*x1*x1 - c)/x1;
+ x = w_sample - w[m];
+
+ //printf("%f %f %f\n", w[0], w[1], w[2]);
+ //printf("%f %f %f %f %f %f\n", x_1, y_1, 0.0, y0, x1, y1);
+ log_amp = a*x*x + b*x + c;
+ //printf("a %f b %f c %f\n", a, b, c);
+ //printf("m=%d A[m-1] %f A[m] %f A[m+1] %f w_sample %f w[m] %f x %f log_amp %f\n", m,
+ // A[m-1],
+ // A[m], A[m+1], w_sample, w[m], x, log_amp);
+ //exit(0);
+ return log_amp;
+}
+
+#define M_MAX 40
+
+float fres[] = {100, 200, 300, 400, 500, 600, 700, 800, 900, 1000,
+ 1200, 1400, 1600, 1850, 2100, 2350, 2600, 2900, 3400, 3800};
+
+/*---------------------------------------------------------------------------*\
+
+ FUNCTION....: resample_amp_nl()
+ AUTHOR......: David Rowe
+ DATE CREATED: 7 March 2011
+
+ Converts the current model with L {Am} samples spaced Wo apart to
+ RES_POINTS samples spaced Wo/RES_POINTS apart. Then subtracts
+ from the previous frames samples to get the delta.
+
+\*---------------------------------------------------------------------------*/
+
+void resample_amp_fixed(MODEL *model,
+ float w[], float A[],
+ float wres[], float Ares[],
+ float AresdB_prev[],
+ float AresdB[],
+ float deltat[])
+{
+ int i;
+
+ for(i=1; i<=model->L; i++) {
+ w[i-1] = i*model->Wo;
+ A[i-1] = model->A[i];
+ }
+
+ for(i=0; i<RES_POINTS; i++) {
+ wres[i] = fres[i]*PI/4000.0;
+ }
+
+ for(i=0; i<RES_POINTS; i++) {
+ Ares[i] = pow(10.0,sample_log_amp_quad_nl(w, A, model->L, wres[i]));
+ }
+
+ /* work out delta T vector for this frame */
+
+ for(i=0; i<RES_POINTS; i++) {
+ AresdB[i] = 20.0*log10(Ares[i]);
+ deltat[i] = AresdB[i] - AresdB_prev[i];
+ }
+
+}
+
+/*---------------------------------------------------------------------------*\
+
+ FUNCTION....: resample_amp_nl()
+ AUTHOR......: David Rowe
+ DATE CREATED: 7 March 2011
+
+ Converts the current model with L {Am} samples spaced Wo apart to M
+ samples spaced Wo/M apart. Then converts back to L {Am} samples.
+ used to prototype constant rate Amplitude encoding ideas.
+
+ Returns the SNR in dB.
+
+\*---------------------------------------------------------------------------*/
+
+float resample_amp_nl(MODEL *model, int m, float AresdB_prev[])
+{
+ int i;
+ float w[MAX_AMP], A[MAX_AMP];
+ float wres[MAX_AMP], Ares[MAX_AMP], AresdB[MAX_AMP];
+ float signal, noise, snr;
+ float new_A;
+ float deltat[MAX_AMP], deltat_q[MAX_AMP], AresdB_q[MAX_AMP];
+
+ resample_amp_fixed(model, w, A, wres, Ares, AresdB_prev, AresdB, deltat);
+
+ /* quantise delta T vector */
+
+ for(i=0; i<RES_POINTS; i++) {
+ noise = 3.0*(1.0 - 2.0*rand()/RAND_MAX);
+ //noise = 0.0;
+ deltat_q[i] = deltat[i] + noise;
+ }
+
+ /* recover Ares vector */
+
+ for(i=0; i<RES_POINTS; i++) {
+ AresdB_q[i] = AresdB_prev[i] + deltat_q[i];
+ Ares[i] = pow(10.0, AresdB_q[i]/20.0);
+ //printf("%d %f %f\n", i, AresdB[i], AresdB_q[i]);
+ }
+
+ /* update memory based on version at decoder */
+
+ for(i=0; i<RES_POINTS; i++) {
+ AresdB_prev[i] = AresdB_q[i];
+ }
+
+
+ dump_resample(wres,Ares,M_MAX);
+
+ signal = noise = 0.0;
+
+ for(i=1; i<model->L; i++) {
+ new_A = pow(10.0,sample_log_amp_quad_nl(wres, Ares, RES_POINTS, model->Wo*i));
+ signal += pow(model->A[i], 2.0);
+ noise += pow(model->A[i] - new_A, 2.0);
+ //printf("%f %f\n", model->A[i], new_A);
+ model->A[i] = new_A;
+ }
+
+ snr = 10.0*log10(signal/noise);
+ printf("snr = %3.2f\n", snr);
+ //exit(0);
+ return snr;
+}
+
+/*---------------------------------------------------------------------------*\
+
+ FUNCTION....: resample_amp()
+ AUTHOR......: David Rowe
+ DATE CREATED: 10 March 2011
+
+ Converts the current model with L {Am} samples spaced Wo apart to M
+ samples with a non-linear spacing. Then converts back to L {Am}
+ samples. used to prototype constant rate Amplitude encoding ideas.
+
+ Returns the SNR in dB.
+
+\*---------------------------------------------------------------------------*/
+
+float resample_amp(MODEL *model, int m)
+{
+ int i;
+ MODEL model_m;
+ float new_A, signal, noise, snr, log_amp_dB;
+ float n_db = 0.0;
+
+ model_m.Wo = PI/(float)m;
+ model_m.L = PI/model_m.Wo;
+
+ for(i=1; i<=model_m.L; i++) {
+ log_amp_dB = 20.0*sample_log_amp_quad(model, i*model_m.Wo);
+ log_amp_dB += n_db*(1.0 - 2.0*rand()/RAND_MAX);
+ model_m.A[i] = pow(10,log_amp_dB/20.0);
+ }
+
+ //dump_resample(&model_m);
+
+ signal = noise = 0.0;
+
+ for(i=1; i<model->L/4; i++) {
+ new_A = pow(10,sample_log_amp_quad(&model_m, i*model->Wo));
+ signal += pow(model->A[i], 2.0);
+ noise += pow(model->A[i] - new_A, 2.0);
+ //printf("%f %f\n", model->A[i], new_A);
+ model->A[i] = new_A;
+ }
+
+ snr = 10.0*log10(signal/noise);
+ //printf("snr = %3.2f\n", snr);
+ //exit(0);
+ return snr;
+}
+
/*---------------------------------------------------------------------------*\
FUNCTION....: interp_lsp()
r = TWO_PI/GLOTTAL_FFT_SIZE;
for(m=1; m<=model->L; m++) {
-
- /* generate excitation */
-
- if (model->voiced) {
- /* I think adding a little jitter helps improve low pitch
- males like hts1a. This moves the onset of each harmonic
- over at +/- 0.25 of a sample.
- */
- jitter = 0.25*(1.0 - 2.0*rand()/RAND_MAX);
- b = floor(m*model->Wo/r + 0.5);
- if (b > ((GLOTTAL_FFT_SIZE/2)-1)) {
- b = (GLOTTAL_FFT_SIZE/2)-1;
- }
- Ex[m].real = cos(ex_phase[0]*m - jitter*model->Wo*m + glottal[b]);
- Ex[m].imag = sin(ex_phase[0]*m - jitter*model->Wo*m + glottal[b]);
- }
- else {
-
- /* When a few samples were tested I found that LPC filter
- phase is not needed in the unvoiced case, but no harm in
- keeping it.
- */
- float phi = TWO_PI*(float)rand()/RAND_MAX;
- Ex[m].real = cos(phi);
- Ex[m].imag = sin(phi);
- }
-
- /* filter using LPC filter */
-
- A_[m].real = H[m].real*Ex[m].real - H[m].imag*Ex[m].imag;
- A_[m].imag = H[m].imag*Ex[m].real + H[m].real*Ex[m].imag;
-
- /* modify sinusoidal phase */
+
+ /* generate excitation */
+
+ if (model->voiced) {
+ /* I think adding a little jitter helps improve low pitch
+ males like hts1a. This moves the onset of each harmonic
+ over at +/- 0.25 of a sample.
+ */
+ jitter = 0.25*(1.0 - 2.0*rand()/RAND_MAX);
+ b = floor(m*model->Wo/r + 0.5);
+ if (b > ((GLOTTAL_FFT_SIZE/2)-1)) {
+ b = (GLOTTAL_FFT_SIZE/2)-1;
+ }
+ Ex[m].real = cos(ex_phase[0]*m - jitter*model->Wo*m + glottal[b]);
+ Ex[m].imag = sin(ex_phase[0]*m - jitter*model->Wo*m + glottal[b]);
+ }
+ else {
+
+ /* When a few samples were tested I found that LPC filter
+ phase is not needed in the unvoiced case, but no harm in
+ keeping it.
+ */
+ float phi = TWO_PI*(float)rand()/RAND_MAX;
+ Ex[m].real = cos(phi);
+ Ex[m].imag = sin(phi);
+ }
+
+ /* filter using LPC filter */
+
+ A_[m].real = H[m].real*Ex[m].real - H[m].imag*Ex[m].imag;
+ A_[m].imag = H[m].imag*Ex[m].real + H[m].real*Ex[m].imag;
+
+ /* modify sinusoidal phase */
- new_phi = atan2(A_[m].imag, A_[m].real+1E-12);
- model->phi[m] = new_phi;
+ new_phi = atan2(A_[m].imag, A_[m].real+1E-12);
+ model->phi[m] = new_phi;
}
}
projects / freetel-svn-tracking.git / commitdiff