pesqmod.c

pesq算法

/*****************************************************************************

Perceptual Evaluation of Speech Quality (PESQ)
ITU-T Draft Recommendation P.862.
Version 1.1 - 15 November 2000.

NOTICE

The Perceptual Evaluation of Speech Quality (PESQ) algorithm and the copyright
therein is the property of British Telecommunications plc and Royal KPN NV,
and is protected by UK, US and other patents.  Permission is granted to use
PESQ for the purpose of evaluation of ITU-T recommendation P.862.  Any other
use of this software or the PESQ algorithm requires a license, which may be
obtained from:

OPTICOM GmbH
Michael Keyhl, Am Weichselgarten 7, D- 91058 Erlangen, Germany
Phone: [phone]     Fax: [phone]     E-mail: [e-mail]

PsyTechnics Limited
Richard Reynolds, B54 Adastral Park, Ipswich IP5 3RE, UK
Phone: [phone] or [phone]    Fax: [phone]
E-mail: [e-mail]

Patent-only licences should be obtained from Opticom. PsyTechnics or Opticom
can provide licences, and further information, for other PESQ products.

Further information is also available from: www.pesq.org

By using this software you acknowledge that PESQ is protected by copyright
and by patents and is being made available to you for the purpose of
evaluation of ITU-T Recommendation P.862.

You must not use PESQ for any other purpose without first obtaining a written
license from British Telecommunications plc and Royal KPN NV, from their
agents listed above.

You must not disclose, reproduce or otherwise release PESQ to any third
party without the prior written permission of British Telecommunications plc
and Royal KPN NV.

Authors:
  Antony Rix (BT) <[e-mail]>
  Mike Hollier (BT)
  Andries Hekstra (KPN Research)
  John Beerends (KPN Research)

*****************************************************************************/

#include <math.h>
#include <stdio.h>
#include "pesq.h"
#include "pesqpar.h"
#include "dsp.h"

#define        CRITERIUM_FOR_SILENCE_OF_5_SAMPLES        500.

float Sl, Sp;



int *nr_of_hz_bands_per_bark_band;
double *centre_of_band_bark;
double *centre_of_band_hz;
double *width_of_band_bark;
double *width_of_band_hz;
double *pow_dens_correction_factor;
double *abs_thresh_power;


void input_filter(
    SIGNAL_INFO * ref_info, SIGNAL_INFO * deg_info, float * ftmp )
{
    DC_block( (*ref_info).data, (*ref_info).Nsamples );
    DC_block( (*deg_info).data, (*deg_info).Nsamples );

    apply_filters( (*ref_info).data, (*ref_info).Nsamples );
    apply_filters( (*deg_info).data, (*deg_info).Nsamples );
}

void calc_VAD( SIGNAL_INFO * sinfo )
{
    apply_VAD( sinfo, sinfo-> data, sinfo-> VAD, sinfo-> logVAD );
}

int id_searchwindows( SIGNAL_INFO * ref_info, SIGNAL_INFO * deg_info,
    ERROR_INFO * err_info )
{
    long  Utt_num = 0;
    long  count, VAD_length;
    long  this_start;
    int   speech_flag = 0;
    float VAD_value;
    long  del_deg_start;
    long  del_deg_end;

    VAD_length = ref_info-> Nsamples / Downsample;

    del_deg_start = MINUTTLENGTH - err_info-> Crude_DelayEst / Downsample;
    del_deg_end =
        ((*deg_info).Nsamples - err_info-> Crude_DelayEst) / Downsample -
        MINUTTLENGTH;

    for (count = 0; count < VAD_length; count++)
    {
        VAD_value = ref_info-> VAD [count];

        if( (VAD_value > 0.0f) && (speech_flag == 0) ) 
        {
            speech_flag = 1;
            this_start = count;
            err_info-> UttSearch_Start [Utt_num] = count - SEARCHBUFFER;
            if( err_info-> UttSearch_Start [Utt_num] < 0 )
                err_info-> UttSearch_Start [Utt_num] = 0;
        }

        if( ((VAD_value == 0.0f) || (count == (VAD_length-1))) &&
            (speech_flag == 1) ) 
        {
            speech_flag = 0;
            err_info-> UttSearch_End [Utt_num] = count + SEARCHBUFFER;
            if( err_info-> UttSearch_End [Utt_num] > VAD_length - 1 )
                err_info-> UttSearch_End [Utt_num] = VAD_length -1;

            if( ((count - this_start) >= MINUTTLENGTH) &&
                (this_start < del_deg_end) &&
                (count > del_deg_start) )
                Utt_num++;            
        }
    }

    err_info-> Nutterances = Utt_num;
    return Utt_num;
} 

void id_utterances( SIGNAL_INFO * ref_info, SIGNAL_INFO * deg_info,
    ERROR_INFO * err_info )
{
    long  Utt_num = 0;
    long  Largest_uttsize = 0;
    long  count, VAD_length;
    int   speech_flag = 0;
    float VAD_value;
    long  this_start;
    long  last_end;
    long  del_deg_start;
    long  del_deg_end;

    VAD_length = ref_info-> Nsamples / Downsample;

    del_deg_start = MINUTTLENGTH - err_info-> Crude_DelayEst / Downsample;
    del_deg_end =
        ((*deg_info).Nsamples - err_info-> Crude_DelayEst) / Downsample -
        MINUTTLENGTH;

    for (count = 0; count < VAD_length ; count++)
    {
        VAD_value = ref_info-> VAD [count];
        if( (VAD_value > 0.0f) && (speech_flag == 0) ) 
        {
            speech_flag = 1;
            this_start = count;
            err_info-> Utt_Start [Utt_num] = count;
        }

        if( ((VAD_value == 0.0f) || (count == (VAD_length-1))) &&
            (speech_flag == 1) ) 
        {
            speech_flag = 0;
            err_info-> Utt_End [Utt_num] = count;

            if( ((count - this_start) >= MINUTTLENGTH) &&
                (this_start < del_deg_end) &&
                (count > del_deg_start) )
                Utt_num++;            
        }
    }

    err_info-> Utt_Start [0] = SEARCHBUFFER;
    err_info-> Utt_End [err_info-> Nutterances-1] = (VAD_length - SEARCHBUFFER);
    
    for (Utt_num = 1; Utt_num < err_info-> Nutterances; Utt_num++ )
    {
        this_start = err_info-> Utt_Start [Utt_num];
        last_end = err_info-> Utt_End [Utt_num - 1];
        count = (this_start + last_end) / 2;
        err_info-> Utt_Start [Utt_num] = count;
        err_info-> Utt_End [Utt_num - 1] = count;
    }

    this_start = (err_info-> Utt_Start [0] * Downsample) + err_info-> Utt_Delay [0];
    if( this_start < (SEARCHBUFFER * Downsample) )
    {
        count = SEARCHBUFFER +
                (Downsample - 1 - err_info-> Utt_Delay [0]) / Downsample;
        err_info-> Utt_Start [0] = count;
    }
    last_end = (err_info-> Utt_End [err_info-> Nutterances-1] * Downsample) +
               err_info-> Utt_Delay [err_info-> Nutterances-1];
    if( last_end > ((*deg_info).Nsamples - SEARCHBUFFER * Downsample) )
    {
        count = ( (*deg_info).Nsamples -
                  err_info-> Utt_Delay [err_info-> Nutterances-1] ) / Downsample -
                SEARCHBUFFER;
        err_info-> Utt_End [err_info-> Nutterances-1] = count;
    }

    for (Utt_num = 1; Utt_num < err_info-> Nutterances; Utt_num++ )
    {
        this_start =
            (err_info-> Utt_Start [Utt_num] * Downsample) +
            err_info-> Utt_Delay [Utt_num];
        last_end =
            (err_info-> Utt_End [Utt_num - 1] * Downsample) +
            err_info-> Utt_Delay [Utt_num - 1];
        if( this_start < last_end )
        {
            count = (this_start + last_end) / 2;
            this_start =
                (Downsample - 1 + count - err_info-> Utt_Delay [Utt_num]) / Downsample;
            last_end =
               (count - err_info-> Utt_Delay [Utt_num - 1]) / Downsample;
            err_info-> Utt_Start [Utt_num] = this_start;
            err_info-> Utt_End [Utt_num - 1] = last_end;
        }
    }

    for (Utt_num = 0; Utt_num < err_info-> Nutterances; Utt_num++ )
        if( (err_info-> Utt_End [Utt_num] - err_info-> Utt_Start [Utt_num])
             > Largest_uttsize )
            Largest_uttsize = 
                err_info-> Utt_End [Utt_num] - err_info-> Utt_Start [Utt_num];

    err_info-> Largest_uttsize = Largest_uttsize;
}

void utterance_split( SIGNAL_INFO * ref_info, SIGNAL_INFO * deg_info,
    ERROR_INFO * err_info, float * ftmp )
{
    long Utt_id;
    long Utt_DelayEst;
    long Utt_Delay;
    float Utt_DelayConf;
    long Utt_Start;
    long Utt_End;
    long Utt_SpeechStart;
    long Utt_SpeechEnd;
    long Utt_Len;
    long step;
    long Best_ED1, Best_ED2;
    long Best_D1, Best_D2;
    float Best_DC1, Best_DC2;
    long Best_BP;
    long Largest_uttsize = 0;

    Utt_id = 0;
    while( (Utt_id < err_info-> Nutterances) &&
           (err_info-> Nutterances < MAXNUTTERANCES) )
    {
        Utt_DelayEst = err_info-> Utt_DelayEst [Utt_id];
        Utt_Delay = err_info-> Utt_Delay [Utt_id];
        Utt_DelayConf = err_info-> Utt_DelayConf [Utt_id];
        Utt_Start = err_info-> Utt_Start [Utt_id];
        Utt_End = err_info-> Utt_End [Utt_id];

        Utt_SpeechStart = Utt_Start;
        while( (Utt_SpeechStart < Utt_End) && (ref_info-> VAD [Utt_SpeechStart] <= 0.0f) )
            Utt_SpeechStart++;
        Utt_SpeechEnd = Utt_End;
        while( (Utt_SpeechEnd > Utt_Start) && (ref_info-> VAD [Utt_SpeechEnd] <= 0.0f) )
            Utt_SpeechEnd--;
        Utt_SpeechEnd++;
        Utt_Len = Utt_SpeechEnd - Utt_SpeechStart;

        if( Utt_Len >= 200 )
        {
            split_align( ref_info, deg_info, err_info, ftmp,
                Utt_Start, Utt_SpeechStart, Utt_SpeechEnd, Utt_End,
                Utt_DelayEst, Utt_DelayConf,
                &Best_ED1, &Best_D1, &Best_DC1,
                &Best_ED2, &Best_D2, &Best_DC2,
                &Best_BP );

            if( (Best_DC1 > Utt_DelayConf) && (Best_DC2 > Utt_DelayConf) )
            {
                for (step = err_info-> Nutterances-1; step > Utt_id; step-- )
                {
                    err_info-> Utt_DelayEst [step +1] = err_info-> Utt_DelayEst [step];
                    err_info-> Utt_Delay [step +1] = err_info-> Utt_Delay [step];
                    err_info-> Utt_DelayConf [step +1] = err_info-> Utt_DelayConf [step];
                    err_info-> Utt_Start [step +1] = err_info-> Utt_Start [step];
                    err_info-> Utt_End [step +1] = err_info-> Utt_End [step];
                    err_info-> UttSearch_Start [step +1] = err_info-> Utt_Start [step];
                    err_info-> UttSearch_End [step +1] = err_info-> Utt_End [step];
                }
                err_info-> Nutterances++;

                err_info-> Utt_DelayEst [Utt_id] = Best_ED1;
                err_info-> Utt_Delay [Utt_id] = Best_D1;
                err_info-> Utt_DelayConf [Utt_id] = Best_DC1;

                err_info-> Utt_DelayEst [Utt_id +1] = Best_ED2;
                err_info-> Utt_Delay [Utt_id +1] = Best_D2;
                err_info-> Utt_DelayConf [Utt_id +1] = Best_DC2;

                err_info-> UttSearch_Start [Utt_id +1] = err_info-> UttSearch_Start [Utt_id];
                err_info-> UttSearch_End [Utt_id +1] = err_info-> UttSearch_End [Utt_id];

                if( Best_D2 < Best_D1 )
                {
                    err_info-> Utt_Start [Utt_id] = Utt_Start;
                    err_info-> Utt_End [Utt_id] = Best_BP;
                    err_info-> Utt_Start [Utt_id +1] = Best_BP;
                    err_info-> Utt_End [Utt_id +1] = Utt_End;
                }
                else
                {
                    err_info-> Utt_Start [Utt_id] = Utt_Start;
                    err_info-> Utt_End [Utt_id] = Best_BP + (Best_D2 - Best_D1) / (2 * Downsample);
                    err_info-> Utt_Start [Utt_id +1] = Best_BP - (Best_D2 - Best_D1) / (2 * Downsample);
                    err_info-> Utt_End [Utt_id +1] = Utt_End;
                }

                if( (err_info-> Utt_Start [Utt_id] - SEARCHBUFFER) * Downsample + Best_D1 < 0 )
                    err_info-> Utt_Start [Utt_id] =
                        SEARCHBUFFER + (Downsample - 1 - Best_D1) / Downsample;

                if( (err_info-> Utt_End [Utt_id +1] * Downsample + Best_D2) >
                    ((*deg_info).Nsamples - SEARCHBUFFER * Downsample) )
                    err_info-> Utt_End [Utt_id +1] =
                        ((*deg_info).Nsamples - Best_D2) / Downsample - SEARCHBUFFER;

            }
            else Utt_id++;
        }
        else Utt_id++;
    }

    for (Utt_id = 0; Utt_id < err_info-> Nutterances; Utt_id++ )
        if( (err_info-> Utt_End [Utt_id] - err_info-> Utt_Start [Utt_id])
             > Largest_uttsize )
            Largest_uttsize =