An even faster, parallel, MATLAB executable (MEX) compilation of the PESQ measure

In a previous post of mine (which you should read now if you haven’t already) I explained how to create a MATLAB executable for the widely used PESQ algorithm. The main reason for wanting to do this was to save time when running a large amount of speech quality tests and the speed increases obtained from using a PESQ MEX function were amazing! At the time of writing that post, the MEX function was approximately 8 times faster than anything else online that I could find. In fact, it is still the fastest implementation I can find, however, I think there is room for improvement and I finally found some time to get it working. In this post, I will show step-by-step how you can compile the PESQ MEX function to accept audio vectors directly from MATLAB and, which, should give great speed increases when run on parallel cores.

First things first, what was wrong with the previous method of compiling the PESQ binaries into a MEX function? Well, for one, it would only accept character arrays describing where the audio files were on disk. The main problem with this is that if you have already processed some audio files in the MATLAB workspace then you would have to write them to disk first before you could pass the file paths to the PESQ MEX function (and then likely delete them after). For a Hard Disk Drive (HDD) that is an unreasonable amount of disk Input/Ouput (I/O) and the HDD’s I/O speed limit then became the bottleneck of the method. For a Solid State Drive (SDD) this is perhaps less of a problem (speed-wise, although the lifetime of a drive is still dependent on the read/write counts), however, it would be good to see the performance increase that an SSD can obtain on a HDD as well.

So what we have is a PESQ MEX function that now needs to be given 1D floating point arrays for the reference signal and degraded signal instead the character array paths used previously. Let’s get stuck in to how to get it working.

From here, I will assume you already have a working PESQ MEX function.

First, lets go to our pesqmain.c file and navigate to the entry point, mexFunction(). We want to declare some new variables for the reference and degraded signals,

float *ref_, *deg_;
size_t refLen, degLen;

after this we want to read all but the last two arguments as the usual regular string arguments.

argv[0] = mxMalloc( sizeof("pesq") );
strcpy (argv[0], "pesq");	
for (i=1; i<(argc-2); i++) {
    argv[i] = mxMalloc( mxGetN(prhs[i-1]) + 1 );
    argv[i] = mxArrayToString(prhs[i-1]);
}

The next part in mexFunction() is where we obtain the last two arguments as audio vectors directly from MATLAB (instead of character array paths),

i = (argc-2);
argv[i] = mxMalloc( 1 );
argv[i] = "v";
ref_ = (float *)mxGetData(prhs[i-1]);
refLen = mxGetNumberOfElements(prhs[i-1]);    
i = (argc-1);
argv[i] = mxMalloc( 1 );
argv[i] = "v";
deg_ = (float *)mxGetData(prhs[i-1]);
degLen = mxGetNumberOfElements(prhs[i-1]);

(note the dummy path "v") and then we will finish off the mexFunction() entry point with:

res = main(argc,argv,ref_,refLen,deg_,degLen);
nlhs = 1;
plhs[0] = mxCreateDoubleScalar(res);
for (i=0; i<(argc-2); i++) {
    mxFree(argv[i]);
}

(notice the extra arguments to main())

Ok, so now that we’ve correctly received the audio vectors into the mex function entry point, we can look at passing these vectors (arrays) on to the PESQ main() function. First we will need to change the function definition and its declaration with:

double main(int argc, char **argv, float *ref_, int refLen, float *deg_, int degLen)

Now find the only call to pesq_measure() from within main() and immediately before the call add the lines:

ref_info.data = ref_;
ref_info.Nsamples = refLen;
deg_info.data = deg_;
deg_info.Nsamples = degLen;

This is where we pass the audio data arrays and their lengths over to the pesq_measure() function using the already defined data structures, ref_info and deg_info.

So that pesq_measure() doesn’t overwrite the arrays with NULL, comment out the first two lines in pesq_measure() which assign NULL values to data:

//ref_info-> data = NULL;
...
//deg_info-> data = NULL;

Awesome, now we are done with the pesqmain.c file!

What happens next is a little more tricky, we need to let the PESQ algorithm know that we already have “read/parsed” the audio data and that it doesn’t need to read from a file anymore. This can be done by removing any file I/O from the C file called pesqio.c.

Open up pesqio.c and change the declaration for input_data and comment out all file I/O lines:

float *input_data;
...
//FILE *Src_file = fopen( sinfo-> path_name, "rb" );
...
/* Comment this big chunk out
if( Src_file == NULL ){
    *Error_Flag = 1;
    *Error_Type = "Could not open source file";
    printf ("%s!\n", *Error_Type);
    safe_free( input_data );
    return;
}
if( fseek( Src_file, 0L, SEEK_END ) != 0 ){
    *Error_Flag = 1;
    *Error_Type = "Could not reach end of source file";
    safe_free( input_data );
    printf ("%s!\n", *Error_Type);
    fclose( Src_file );
    return;
}
file_size = ftell( Src_file );
if( file_size < 0L ){ 
    *Error_Flag = 1;
    *Error_Type = "Could not measure length of source file";
    safe_free( input_data );
    printf ("%s!\n", *Error_Type);
    fclose( Src_file );
    return;
}
if( fseek( Src_file, 0L, SEEK_SET ) != 0 ){
    *Error_Flag = 1;
    *Error_Type = "Could not reach start of source file";
    safe_free( input_data );
    printf ("%s!\n", *Error_Type);
    fclose( Src_file );
    return;
}
name_len = strlen( sinfo-> path_name );
if( name_len > 4 ){
    if( strcmp( sinfo-> path_name + name_len - 4, ".wav" ) == 0 )
        header_size = 22;
    if( strcmp( sinfo-> path_name + name_len - 4, ".WAV" ) == 0 )
        header_size = 22;
    if( strcmp( sinfo-> path_name + name_len - 4, ".raw" ) == 0 )
        header_size = 0;
    if( strcmp( sinfo-> path_name + name_len - 4, ".src" ) == 0 )
        header_size = 0;
}
if( name_len > 2 ){
    if( strcmp( sinfo-> path_name + name_len - 2, ".s" ) == 0 )
        header_size = 0;
}
if( header_size > 0 )
    fread( input_data, 2, header_size, Src_file );
*/
...
//fclose( Src_file ); <- comment all of these out, there's a few here and there

Good, now lets replace this line:

input_data = (short *) safe_malloc( 16384 * sizeof(short) );

with these:

input_data = (float *) safe_malloc( sinfo-> Nsamples * sizeof(float) );    
memcpy( input_data, sinfo->data, sinfo-> Nsamples * sizeof(float) );

so that we’ve copied the audio data to the same place we would read our file data from, input_data.

We also want to overwrite the Nsamples variable with the already known number of samples, so replace this line:

Nsamples = (file_size / 2) - header_size;

with this line:

Nsamples = sinfo-> Nsamples;

Now because we already have the data and we don’t actually need to read from a file, we should comment out (or replace) all of this:

/*
to_read = Nsamples;
    while( to_read > 16384 ){
        read_count = fread( input_data, sizeof(short), 16384, Src_file );
        if( read_count < 16384 ) {
            *Error_Flag = 1;
            *Error_Type = "Error reading source file.";
            printf ("%s!\n", *Error_Type);
            safe_free( input_data );
            safe_free( sinfo-> data );
            sinfo-> data = NULL;
            fclose( Src_file );
            return;
        }
        if( sinfo-> apply_swap ){
            p_byte = (char *)input_data;
            for( count = 0L; count < read_count; count++ ){
                s = p_byte[count << 1];
                p_byte[count << 1] = p_byte[(count << 1)+1];
                p_byte[(count << 1)+1] = s;
            }
        }
        to_read -= read_count;
        p_input = input_data;
        while( read_count > 0 ){
            read_count--;
            *(read_ptr++) = (float)(*(p_input++));
        }
    }
    read_count = fread( input_data, sizeof(short), to_read, Src_file );
    if( read_count < to_read ) {
        *Error_Flag = 1;
        *Error_Type = "Error reading source file";
        printf ("%s!\n", *Error_Type);
        safe_free( input_data );
        safe_free( sinfo-> data );
        sinfo-> data = NULL;
        fclose( Src_file );
        return;
    }
    if( sinfo-> apply_swap ){
        p_byte = (char *)input_data;
        for( count = 0L; count < read_count; count++ ){
            s = p_byte[count << 1];
            p_byte[count << 1] = p_byte[(count << 1)+1];
            p_byte[(count << 1)+1] = s;
        }
    }
    p_input = input_data;
    while( read_count > 0 ){
        read_count--;
        *(read_ptr++) = (float)(*(p_input++));
    }
*/

with this:

read_count = Nsamples;
while( read_count > 0 ){
    read_count--;
    *(read_ptr++) = (float)(*(input_data++));
}

and one last part, we don’t want to free any of our MATLAB data so let’s comment this out:

//safe_free( input_data );

That’s it! Let’s go ahead and compile the new MEX function with the following:

mex('pesqmain.c', 'pesqmod.c', 'pesqio.c', 'pesqdsp.c', 'dsp.c')

Now you can call this MEX function from MATLAB and instead of using character arrays you can send single type MATLAB vectors directly over to the MEX function as the last two arguments!

If you’re unsure of how to do this I have written a nice little helper function which you can download, a description is included in the source code.

My tip would be to try running this new MEX function (or helper function) from within a parfor loop in MATLAB! Or if you know how to run algorithms on the GPU in MATLAB, give that a go! and let me know what sort of performance enhancements you get in the comments.  🙂

Enjoy!