pesqmod.c

pesq算法

                err_info-> Utt_End [Utt_id] - err_info-> Utt_Start [Utt_id];

    err_info-> Largest_uttsize = Largest_uttsize;
}

void utterance_locate( SIGNAL_INFO * ref_info, SIGNAL_INFO * deg_info,
    ERROR_INFO * err_info, float * ftmp )
{    
    long Utt_id;
    
    id_searchwindows( ref_info, deg_info, err_info );

    for (Utt_id = 0; Utt_id < err_info-> Nutterances; Utt_id++)
    {
        crude_align( ref_info, deg_info, err_info, Utt_id, ftmp);
        time_align(ref_info, deg_info, err_info, Utt_id, ftmp );
    }

    id_utterances( ref_info, deg_info, err_info );

    utterance_split( ref_info, deg_info, err_info, ftmp );   
}


void short_term_fft (int Nf, SIGNAL_INFO *info, float *window, long start_sample, float *hz_spectrum, float *fft_tmp) {
    int n, k;        

    for (n = 0; n < Nf; n++ )
    {
        fft_tmp [n] = info-> data [start_sample + n] * window [n];
    }
    RealFFT(fft_tmp, Nf);

    for (k = 0; k < Nf / 2; k++ ) 
    {
        hz_spectrum [k] = fft_tmp [k << 1] * fft_tmp [k << 1] + fft_tmp [1 + (k << 1)] * fft_tmp [1 + (k << 1)];
    }    

    hz_spectrum [0] = 0;
}

void freq_warping (int number_of_hz_bands, float *hz_spectrum, int Nb, float *pitch_pow_dens, long frame) {

    int        hz_band = 0;
    int        bark_band;
    double    sum;

    for (bark_band = 0; bark_band < Nb; bark_band++) {
        int n = nr_of_hz_bands_per_bark_band [bark_band];
        int i;

        sum = 0;
        for (i = 0; i < n; i++) {
            sum += hz_spectrum [hz_band++];
        }
        
        sum *= pow_dens_correction_factor [bark_band];
        sum *= Sp;
        pitch_pow_dens [frame * Nb + bark_band] = (float) sum;
    }
}

float total_audible (int frame, float *pitch_pow_dens, float factor) {
    int        band;
    float     h, threshold;
    double  result;
    
    result = 0.;
    for (band= 1; band< Nb; band++) {
        h = pitch_pow_dens [frame * Nb + band];
        threshold = (float) (factor * abs_thresh_power [band]);
        if (h > threshold) {
            result += h;
        }
    }
    return (float) result;
}

void time_avg_audible_of (int number_of_frames, int *silent, float *pitch_pow_dens, float *avg_pitch_pow_dens, int total_number_of_frames) 
{
    int    frame;
    int    band;

    for (band = 0; band < Nb; band++) {
        double result = 0;
        for (frame = 0; frame < number_of_frames; frame++) {
            if (!silent [frame]) {
                float h = pitch_pow_dens [frame * Nb + band];
                if (h > 100 * abs_thresh_power [band]) {
                    result += h;
                }
            }
        }

        avg_pitch_pow_dens [band] = (float) (result / total_number_of_frames);
    }
}            

void freq_resp_compensation (int number_of_frames, float *pitch_pow_dens_ref, float *avg_pitch_pow_dens_ref, float *avg_pitch_pow_dens_deg, float constant)
{
    int band;

    for (band = 0; band < Nb; band++) {
        float    x = (avg_pitch_pow_dens_deg [band] + constant) / (avg_pitch_pow_dens_ref [band] + constant);
        int        frame;

        if (x > (float) 100.0) {x = (float) 100.0;} 
        if (x < (float) 0.01) {x = (float) 0.01;}   

        for (frame = 0; frame < number_of_frames; frame++) {        
            pitch_pow_dens_ref [frame * Nb + band] *= x;
        }        
    }
}

#define ZWICKER_POWER       0.23 

void intensity_warping_of (float *loudness_dens, int frame, float *pitch_pow_dens)
{
    int        band;
    float    h;
    double    modified_zwicker_power;

    for (band = 0; band < Nb; band++) {
        float threshold = (float) abs_thresh_power [band];
        float input = pitch_pow_dens [frame * Nb + band];

        if (centre_of_band_bark [band] < (float) 4) {
            h =  (float) 6 / ((float) centre_of_band_bark [band] + (float) 2);
        } else {
            h = (float) 1;
        }
        if (h > (float) 2) {h = (float) 2;}
        h = (float) pow (h, (float) 0.15); 
        modified_zwicker_power = ZWICKER_POWER * h;

        if (input > threshold) {
            loudness_dens [band] = (float) (pow (threshold / 0.5, modified_zwicker_power)
                                                    * (pow (0.5 + 0.5 * input / threshold, modified_zwicker_power) - 1));
        } else {
            loudness_dens [band] = 0;
        }

        loudness_dens [band] *= (float) Sl;
    }    
}

float pseudo_Lp (int n, float *x, float p) {   
    double totalWeight = 0;
    double result = 0;
    int    band;

    for (band = 1; band < Nb; band++) {
        float h = (float) fabs (x [band]);        
        float w = (float) width_of_band_bark [band];
        float prod = h * w;

        result += pow (prod, p);
        totalWeight += w;
    }

    result /= totalWeight;
    result = pow (result, 1/p);
    result *= totalWeight;
    
    return (float) result;
}  
void multiply_with_asymmetry_factor (float      *disturbance_dens, 
                                     int         frame, 
                                     const float   * const pitch_pow_dens_ref, 
                                     const float   * const pitch_pow_dens_deg) 
{
    int   i;
    float ratio, h;

    for (i = 0; i < Nb; i++) {
        ratio = (pitch_pow_dens_deg [frame * Nb + i] + (float) 50)
                  / (pitch_pow_dens_ref [frame * Nb + i] + (float) 50);

        h = (float) pow (ratio, (float) 1.2);    
        if (h > (float) 12) {h = (float) 12;}
        if (h < (float) 3) {h = (float) 0.0;}

        disturbance_dens [i] *= h;
    }
}

double pow_of (const float * const x, long start_sample, long stop_sample, long divisor) {
    long    i;
    double  power = 0;

    if (start_sample < 0) {
        exit (1);
    }

    if (start_sample > stop_sample) {
        exit (1);
    }

    for (i = start_sample; i < stop_sample; i++) {
        float h = x [i];
        power += h * h;        
    }
    
    power /= divisor;
    return power;
}


int compute_delay (long              start_sample, 
                   long                 stop_sample, 
                   long                 search_range, 
                   float            *time_series1, 
                   float            *time_series2,
                   float            *max_correlation) {

    double            power1, power2, normalization;
    long            i;
    float           *x1, *x2, *y;
    double            h;
    long            n = stop_sample - start_sample;   
    long            power_of_2 = nextpow2 (2 * n);
    long            best_delay;

    power1 = pow_of (time_series1, start_sample, stop_sample, stop_sample - start_sample) * (double) n/(double) power_of_2;
    power2 = pow_of (time_series2, start_sample, stop_sample, stop_sample - start_sample) * (double) n/(double) power_of_2;
    normalization = sqrt (power1 * power2);

    if ((power1 <= 1E-6) || (power2 <= 1E-6)) {
        *max_correlation = 0;
        return 0;
    }

    x1 = (float *) safe_malloc ((power_of_2 + 2) * sizeof (float));;
    x2 = (float *) safe_malloc ((power_of_2 + 2) * sizeof (float));;
    y = (float *) safe_malloc ((power_of_2 + 2) * sizeof (float));;
    
    for (i = 0; i < power_of_2 + 2; i++) {
        x1 [i] = 0.;
        x2 [i] = 0.;
        y [i] = 0.;
    }

    for (i = 0; i < n; i++) {
        x1 [i] = (float) fabs (time_series1 [i + start_sample]);
        x2 [i] = (float) fabs (time_series2 [i + start_sample]);
    }

    RealFFT (x1, power_of_2);
    RealFFT (x2, power_of_2);

    for (i = 0; i <= power_of_2 / 2; i++) { 
        x1 [2 * i] /= power_of_2;
        x1 [2 * i + 1] /= power_of_2;                
    }

    for (i = 0; i <= power_of_2 / 2; i++) { 
        y [2*i] = x1 [2*i] * x2 [2*i] + x1 [2*i + 1] * x2 [2*i + 1];
        y [2*i + 1] = -x1 [2*i + 1] * x2 [2*i] + x1 [2*i] * x2 [2*i + 1];
    }    
  
    RealIFFT (y, power_of_2);

    best_delay = 0;
    *max_correlation = 0;

    for (i = -search_range; i <= -1; i++) {
        h = (float) fabs (y [(i + power_of_2)]) / normalization;
        if (fabs (h) > (double) *max_correlation) {
            *max_correlation = (float) fabs (h);
            best_delay= i;
        }
    }

    for (i = 0; i < search_range; i++) {
        h = (float) fabs (y [i]) / normalization;
        if (fabs (h) > (double) *max_correlation) {
            *max_correlation = (float) fabs (h);
            best_delay= i;
        }
    }

    safe_free (x1);
    safe_free (x2);
    safe_free (y);
    
    return best_delay;
}

    
#define NUMBER_OF_PSQM_FRAMES_PER_SYLLABE       20
 

float Lpq_weight (int         start_frame,
                  int         stop_frame,
                  float         power_syllable,
                  float      power_time,
                  float        *frame_disturbance,
                  float        *time_weight) {

    double    result_time= 0;
    double  total_time_weight_time = 0;
    int        start_frame_of_syllable;
    
    for (start_frame_of_syllable = start_frame; 
         start_frame_of_syllable <= stop_frame; 
         start_frame_of_syllable += NUMBER_OF_PSQM_FRAMES_PER_SYLLABE/2) {

        double  result_syllable = 0;
        int     count_syllable = 0;
        int     frame;

        for (frame = start_frame_of_syllable;
             frame < start_frame_of_syllable + NUMBER_OF_PSQM_FRAMES_PER_SYLLABE;
             frame++) {
            if (frame <= stop_frame) {
                float h = frame_disturbance [frame];
                result_syllable +=  pow (h, power_syllable); 
            }
            count_syllable++;                
        }

        result_syllable /= count_syllable;
        result_syllable = pow (result_syllable, (double) 1/power_syllable);        
    
        result_time+=  pow (time_weight [start_frame_of_syllable - start_frame] * result_syllable, power_time); 
        total_time_weight_time += pow (time_weight [start_frame_of_syllable - start_frame], power_time);
    }

    result_time /= total_time_weight_time;
    result_time= pow (result_time, (float) 1 / power_time);

    return (float) result_time;
}

void set_to_sine (SIGNAL_INFO *info, float amplitude, float omega) {
    long i;

    for (i = 0; i < info-> Nsamples; i++) {
        info-> data [i] = amplitude * (float) sin (omega * i);
    }
}