From 30325d24d107dbf133da39f7c96d1510fd1c9449 Mon Sep 17 00:00:00 2001 From: erdgeist Date: Fri, 15 Aug 2025 12:42:40 +0200 Subject: Bump to codec2 version 1.2.0 --- sine.c | 591 +++++++++++++++++++++++++++++++---------------------------------- 1 file changed, 286 insertions(+), 305 deletions(-) (limited to 'sine.c') diff --git a/sine.c b/sine.c index d912cd7..814c269 100644 --- a/sine.c +++ b/sine.c @@ -27,64 +27,66 @@ /*---------------------------------------------------------------------------*\ - INCLUDES + INCLUDES \*---------------------------------------------------------------------------*/ -#include -#include +#include "sine.h" + #include +#include +#include #include "defines.h" -#include "sine.h" #include "kiss_fft.h" #define HPF_BETA 0.125 /*---------------------------------------------------------------------------*\ - HEADERS + HEADERS \*---------------------------------------------------------------------------*/ void hs_pitch_refinement(MODEL *model, COMP Sw[], float pmin, float pmax, - float pstep); + float pstep); /*---------------------------------------------------------------------------*\ - FUNCTIONS + FUNCTIONS \*---------------------------------------------------------------------------*/ C2CONST c2const_create(int Fs, float framelength_s) { - C2CONST c2const; - - assert((Fs == 8000) || (Fs = 16000)); - c2const.Fs = Fs; - c2const.n_samp = round(Fs*framelength_s); - c2const.max_amp = floor(Fs*P_MIN_S/2); - c2const.p_min = floor(Fs*P_MIN_S); - c2const.p_max = floor(Fs*P_MAX_S); - c2const.m_pitch = floor(Fs*M_PITCH_S); - c2const.Wo_min = TWO_PI/c2const.p_max; - c2const.Wo_max = TWO_PI/c2const.p_min; - - if (Fs == 8000) { - c2const.nw = 279; - } else { - c2const.nw = 511; /* actually a bit shorter in time but lets us maintain constant FFT size */ - } + C2CONST c2const; + + assert((Fs == 8000) || (Fs == 16000)); + c2const.Fs = Fs; + c2const.n_samp = round(Fs * framelength_s); + c2const.max_amp = floor(Fs * P_MAX_S / 2); + c2const.p_min = floor(Fs * P_MIN_S); + c2const.p_max = floor(Fs * P_MAX_S); + c2const.m_pitch = floor(Fs * M_PITCH_S); + c2const.Wo_min = TWO_PI / c2const.p_max; + c2const.Wo_max = TWO_PI / c2const.p_min; + + if (Fs == 8000) { + c2const.nw = 279; + } else { + c2const.nw = 511; /* actually a bit shorter in time but lets us maintain + constant FFT size */ + } - c2const.tw = Fs*TW_S; + c2const.tw = Fs * TW_S; - /* - fprintf(stderr, "max_amp: %d m_pitch: %d\n", c2const.n_samp, c2const.m_pitch); - fprintf(stderr, "p_min: %d p_max: %d\n", c2const.p_min, c2const.p_max); - fprintf(stderr, "Wo_min: %f Wo_max: %f\n", c2const.Wo_min, c2const.Wo_max); - fprintf(stderr, "nw: %d tw: %d\n", c2const.nw, c2const.tw); - */ + /* + fprintf(stderr, "max_amp: %d m_pitch: %d\n", c2const.n_samp, c2const.m_pitch); + fprintf(stderr, "p_min: %d p_max: %d\n", c2const.p_min, c2const.p_max); + fprintf(stderr, "Wo_min: %f Wo_max: %f\n", c2const.Wo_min, c2const.Wo_max); + fprintf(stderr, "nw: %d tw: %d\n", c2const.nw, c2const.tw); + */ - return c2const; + return c2const; } /*---------------------------------------------------------------------------*\ @@ -97,13 +99,13 @@ C2CONST c2const_create(int Fs, float framelength_s) { \*---------------------------------------------------------------------------*/ -void make_analysis_window(C2CONST *c2const, codec2_fft_cfg fft_fwd_cfg, float w[], float W[]) -{ +void make_analysis_window(C2CONST *c2const, codec2_fft_cfg fft_fwd_cfg, + float w[], float W[]) { float m; - COMP wshift[FFT_ENC]; - int i,j; - int m_pitch = c2const->m_pitch; - int nw = c2const->nw; + COMP wshift[FFT_ENC]; + int i, j; + int m_pitch = c2const->m_pitch; + int nw = c2const->nw; /* Generate Hamming window centered on M-sample pitch analysis window @@ -117,20 +119,18 @@ void make_analysis_window(C2CONST *c2const, codec2_fft_cfg fft_fwd_cfg, float w[ */ m = 0.0; - for(i=0; im_pitch; - int nw = c2const->nw; - - for(i=0; im_pitch; + int nw = c2const->nw; + + for (i = 0; i < FFT_ENC; i++) { + Sw[i].real = 0.0; + Sw[i].imag = 0.0; + } - /* Centre analysis window on time axis, we need to arrange input - to FFT this way to make FFT phases correct */ + /* Centre analysis window on time axis, we need to arrange input + to FFT this way to make FFT phases correct */ - /* move 2nd half to start of FFT input vector */ + /* move 2nd half to start of FFT input vector */ - for(i=0; iWo + 5; - pmin = TWO_PI/model->Wo - 5; + pmax = TWO_PI / model->Wo + 5; + pmin = TWO_PI / model->Wo - 5; pstep = 1.0; - hs_pitch_refinement(model,Sw,pmin,pmax,pstep); + hs_pitch_refinement(model, Sw, pmin, pmax, pstep); /* Fine refinement */ - pmax = TWO_PI/model->Wo + 1; - pmin = TWO_PI/model->Wo - 1; + pmax = TWO_PI / model->Wo + 1; + pmin = TWO_PI / model->Wo - 1; pstep = 0.25; - hs_pitch_refinement(model,Sw,pmin,pmax,pstep); + hs_pitch_refinement(model, Sw, pmin, pmax, pstep); /* Limit range */ - if (model->Wo < TWO_PI/c2const->p_max) - model->Wo = TWO_PI/c2const->p_max; - if (model->Wo > TWO_PI/c2const->p_min) - model->Wo = TWO_PI/c2const->p_min; + if (model->Wo < TWO_PI / c2const->p_max) model->Wo = TWO_PI / c2const->p_max; + if (model->Wo > TWO_PI / c2const->p_min) model->Wo = TWO_PI / c2const->p_min; - model->L = floorf(PI/model->Wo); + model->L = floorf(PI / model->Wo); /* trap occasional round off issues with floorf() */ - if (model->Wo*model->L >= 0.95*PI) { - model->L--; + if (model->Wo * model->L >= 0.95 * PI) { + model->L--; } - assert(model->Wo*model->L < PI); + assert(model->Wo * model->L < PI); } /*---------------------------------------------------------------------------*\ @@ -348,35 +342,39 @@ void two_stage_pitch_refinement(C2CONST *c2const, MODEL *model, COMP Sw[]) \*---------------------------------------------------------------------------*/ -void hs_pitch_refinement(MODEL *model, COMP Sw[], float pmin, float pmax, float pstep) -{ - int m; /* loop variable */ - int b; /* bin for current harmonic centre */ - float E; /* energy for current pitch*/ - float Wo; /* current "test" fundamental freq. */ - float Wom; /* Wo that maximises E */ - float Em; /* mamimum energy */ - float r, one_on_r; /* number of rads/bin */ - float p; /* current pitch */ +void hs_pitch_refinement(MODEL *model, COMP Sw[], float pmin, float pmax, + float pstep) { + int m; /* loop variable */ + int b; /* bin for current harmonic centre */ + float E; /* energy for current pitch*/ + float Wo; /* current "test" fundamental freq. */ + float Wom; /* Wo that maximises E */ + float Em; /* mamimum energy */ + float r, one_on_r; /* number of rads/bin */ + float p; /* current pitch */ /* Initialisation */ - model->L = PI/model->Wo; /* use initial pitch est. for L */ + model->L = PI / model->Wo; /* use initial pitch est. for L */ Wom = model->Wo; Em = 0.0; - r = TWO_PI/FFT_ENC; - one_on_r = 1.0/r; + r = TWO_PI / FFT_ENC; + one_on_r = 1.0 / r; /* Determine harmonic sum for a range of Wo values */ - for(p=pmin; p<=pmax; p+=pstep) { + for (p = pmin; p <= pmax; p += pstep) { E = 0.0; - Wo = TWO_PI/p; + Wo = TWO_PI / p; + + float bFloat = Wo * one_on_r; + float currentBFloat = bFloat; /* Sum harmonic magnitudes */ - for(m=1; m<=model->L; m++) { - b = (int)(m*Wo*one_on_r + 0.5); - E += Sw[b].real*Sw[b].real + Sw[b].imag*Sw[b].imag; + for (m = 1; m <= model->L; m++) { + b = (int)(currentBFloat + 0.5); + E += Sw[b].real * Sw[b].real + Sw[b].imag * Sw[b].imag; + currentBFloat += bFloat; } /* Compare to see if this is a maximum */ @@ -399,35 +397,35 @@ void hs_pitch_refinement(MODEL *model, COMP Sw[], float pmin, float pmax, float \*---------------------------------------------------------------------------*/ -void estimate_amplitudes(MODEL *model, COMP Sw[], float W[], int est_phase) -{ - int i,m; /* loop variables */ - int am,bm; /* bounds of current harmonic */ - float den; /* denominator of amplitude expression */ +void estimate_amplitudes(MODEL *model, COMP Sw[], float W[], int est_phase) { + int i, m; /* loop variables */ + int am, bm; /* bounds of current harmonic */ + float den; /* denominator of amplitude expression */ - float r = TWO_PI/FFT_ENC; - float one_on_r = 1.0/r; + float r = TWO_PI / FFT_ENC; + float one_on_r = 1.0 / r; - for(m=1; m<=model->L; m++) { + for (m = 1; m <= model->L; m++) { /* Estimate ampltude of harmonic */ den = 0.0; - am = (int)((m - 0.5)*model->Wo*one_on_r + 0.5); - bm = (int)((m + 0.5)*model->Wo*one_on_r + 0.5); + am = (int)((m - 0.5) * model->Wo * one_on_r + 0.5); + bm = (int)((m + 0.5) * model->Wo * one_on_r + 0.5); - for(i=am; iA[m] = sqrtf(den); if (est_phase) { - int b = (int)(m*model->Wo/r + 0.5); /* DFT bin of centre of current harmonic */ + int b = (int)(m * model->Wo / r + + 0.5); /* DFT bin of centre of current harmonic */ - /* Estimate phase of harmonic, this is expensive in CPU for - embedded devicesso we make it an option */ + /* Estimate phase of harmonic, this is expensive in CPU for + embedded devicesso we make it an option */ - model->phi[m] = atan2f(Sw[b].imag,Sw[b].real); + model->phi[m] = atan2f(Sw[b].imag, Sw[b].real); } } } @@ -443,119 +441,110 @@ void estimate_amplitudes(MODEL *model, COMP Sw[], float W[], int est_phase) \*---------------------------------------------------------------------------*/ -float est_voicing_mbe( - C2CONST *c2const, - MODEL *model, - COMP Sw[], - float W[] - ) -{ - int l,al,bl,m; /* loop variables */ - COMP Am; /* amplitude sample for this band */ - int offset; /* centers Hw[] about current harmonic */ - float den; /* denominator of Am expression */ - float error; /* accumulated error between original and synthesised */ - float Wo; - float sig, snr; - float elow, ehigh, eratio; - float sixty; - COMP Ew; - Ew.real = 0; - Ew.imag = 0; - - int l_1000hz = model->L*1000.0/(c2const->Fs/2); - sig = 1E-4; - for(l=1; l<=l_1000hz; l++) { - sig += model->A[l]*model->A[l]; - } +float est_voicing_mbe(C2CONST *c2const, MODEL *model, COMP Sw[], float W[]) { + int l, al, bl, m; /* loop variables */ + COMP Am; /* amplitude sample for this band */ + int offset; /* centers Hw[] about current harmonic */ + float den; /* denominator of Am expression */ + float error; /* accumulated error between original and synthesised */ + float Wo; + float sig, snr; + float elow, ehigh, eratio; + float sixty; + COMP Ew; + Ew.real = 0; + Ew.imag = 0; + + int l_1000hz = model->L * 1000.0 / (c2const->Fs / 2); + sig = 1E-4; + for (l = 1; l <= l_1000hz; l++) { + sig += model->A[l] * model->A[l]; + } - Wo = model->Wo; - error = 1E-4; + Wo = model->Wo; + error = 1E-4; - /* Just test across the harmonics in the first 1000 Hz */ + /* Just test across the harmonics in the first 1000 Hz */ - for(l=1; l<=l_1000hz; l++) { - Am.real = 0.0; - Am.imag = 0.0; - den = 0.0; - al = ceilf((l - 0.5)*Wo*FFT_ENC/TWO_PI); - bl = ceilf((l + 0.5)*Wo*FFT_ENC/TWO_PI); + for (l = 1; l <= l_1000hz; l++) { + Am.real = 0.0; + Am.imag = 0.0; + den = 0.0; + al = ceilf((l - 0.5) * Wo * FFT_ENC / TWO_PI); + bl = ceilf((l + 0.5) * Wo * FFT_ENC / TWO_PI); - /* Estimate amplitude of harmonic assuming harmonic is totally voiced */ + /* Estimate amplitude of harmonic assuming harmonic is totally voiced */ - offset = FFT_ENC/2 - l*Wo*FFT_ENC/TWO_PI + 0.5; - for(m=al; m V_THRESH) - model->voiced = 1; - else - model->voiced = 0; - - /* post processing, helps clean up some voicing errors ------------------*/ - - /* - Determine the ratio of low freqency to high frequency energy, - voiced speech tends to be dominated by low frequency energy, - unvoiced by high frequency. This measure can be used to - determine if we have made any gross errors. - */ - - int l_2000hz = model->L*2000.0/(c2const->Fs/2); - int l_4000hz = model->L*4000.0/(c2const->Fs/2); - elow = ehigh = 1E-4; - for(l=1; l<=l_2000hz; l++) { - elow += model->A[l]*model->A[l]; - } - for(l=l_2000hz; l<=l_4000hz; l++) { - ehigh += model->A[l]*model->A[l]; - } - eratio = 10.0*log10f(elow/ehigh); + snr = 10.0 * log10f(sig / error); + if (snr > V_THRESH) + model->voiced = 1; + else + model->voiced = 0; + + /* post processing, helps clean up some voicing errors ------------------*/ - /* Look for Type 1 errors, strongly V speech that has been - accidentally declared UV */ + /* + Determine the ratio of low frequency to high frequency energy, + voiced speech tends to be dominated by low frequency energy, + unvoiced by high frequency. This measure can be used to + determine if we have made any gross errors. + */ + + int l_2000hz = model->L * 2000.0 / (c2const->Fs / 2); + int l_4000hz = model->L * 4000.0 / (c2const->Fs / 2); + elow = ehigh = 1E-4; + for (l = 1; l <= l_2000hz; l++) { + elow += model->A[l] * model->A[l]; + } + for (l = l_2000hz; l <= l_4000hz; l++) { + ehigh += model->A[l] * model->A[l]; + } + eratio = 10.0 * log10f(elow / ehigh); - if (model->voiced == 0) - if (eratio > 10.0) - model->voiced = 1; + /* Look for Type 1 errors, strongly V speech that has been + accidentally declared UV */ - /* Look for Type 2 errors, strongly UV speech that has been - accidentally declared V */ + if (model->voiced == 0) + if (eratio > 10.0) model->voiced = 1; - if (model->voiced == 1) { - if (eratio < -10.0) - model->voiced = 0; + /* Look for Type 2 errors, strongly UV speech that has been + accidentally declared V */ - /* A common source of Type 2 errors is the pitch estimator - gives a low (50Hz) estimate for UV speech, which gives a - good match with noise due to the close harmoonic spacing. - These errors are much more common than people with 50Hz3 - pitch, so we have just a small eratio threshold. */ + if (model->voiced == 1) { + if (eratio < -10.0) model->voiced = 0; - sixty = 60.0*TWO_PI/c2const->Fs; - if ((eratio < -4.0) && (model->Wo <= sixty)) - model->voiced = 0; - } - //printf(" v: %d snr: %f eratio: %3.2f %f\n",model->voiced,snr,eratio,dF0); + /* A common source of Type 2 errors is the pitch estimator + gives a low (50Hz) estimate for UV speech, which gives a + good match with noise due to the close harmoonic spacing. + These errors are much more common than people with 50Hz3 + pitch, so we have just a small eratio threshold. */ - return snr; + sixty = 60.0 * TWO_PI / c2const->Fs; + if ((eratio < -4.0) && (model->Wo <= sixty)) model->voiced = 0; + } + // printf(" v: %d snr: %f eratio: %3.2f %f\n",model->voiced,snr,eratio,dF0); + + return snr; } /*---------------------------------------------------------------------------*\ @@ -564,32 +553,29 @@ float est_voicing_mbe( AUTHOR......: David Rowe DATE CREATED: 11/5/94 - Init function that generates the trapezoidal (Parzen) sythesis window. + Init function that generates the trapezoidal (Parzen) synthesis window. \*---------------------------------------------------------------------------*/ -void make_synthesis_window(C2CONST *c2const, float Pn[]) -{ - int i; +void make_synthesis_window(C2CONST *c2const, float Pn[]) { + int i; float win; - int n_samp = c2const->n_samp; - int tw = c2const->tw; + int n_samp = c2const->n_samp; + int tw = c2const->tw; /* Generate Parzen window in time domain */ win = 0.0; - for(i=0; iL; l++) { - b = (int)(l*model->Wo*FFT_DEC/TWO_PI + 0.5); - if (b > ((FFT_DEC/2)-1)) { - b = (FFT_DEC/2)-1; - } - Sw_[b].real = model->A[l]*cosf(model->phi[l]); - Sw_[b].imag = model->A[l]*sinf(model->phi[l]); + for (l = 1; l <= model->L; l++) { + b = (int)(l * model->Wo * FFT_DEC / TWO_PI + 0.5); + if (b > ((FFT_DEC / 2) - 1)) { + b = (FFT_DEC / 2) - 1; } + Sw_[b].real = model->A[l] * cosf(model->phi[l]); + Sw_[b].imag = model->A[l] * sinf(model->phi[l]); + } - /* Perform inverse DFT */ + /* Perform inverse DFT */ - codec2_fftri(fftr_inv_cfg, Sw_,sw_); + codec2_fftri(fftr_inv_cfg, Sw_, sw_); - /* Overlap add to previous samples */ + /* Overlap add to previous samples */ - #ifdef USE_KISS_FFT - #define FFTI_FACTOR ((float)1.0) - #else - #define FFTI_FACTOR ((float32_t)FFT_DEC) - #endif +#ifdef USE_KISS_FFT +#define FFTI_FACTOR ((float)1.0) +#else +#define FFTI_FACTOR ((float32_t)FFT_DEC) +#endif - for(i=0; i