Avoiding "RuntimeWarning: invalid value encountered in log10" when converting FFT data to dBFS

Question

I've attached some code below that I'm having some issues with. Basically, I need to write a function to get the dBFS value of a signal at a specified Hz value, at a given channel in a .wav file. For a number of reasons, I need this to work on both 16 bit integer and 32 bit integer formats. Since SciPy's .wav io module supports 16 bit and 32 bit integer, I'm hoping that I could just convert either to 32bit float first, then perform the FFT and dBFS conversions. For example, if you have a stereo test file that has a 440 Hz tone on channel zero, the function below would be called by db_val = check_db_at_freq_channel(test_file, 440, 0) where test_file is just a file path.

I have a test file that was exported as both 16 bit int, and a copy at 32 bit int, when I run this function on the 16 bit file, it's fine, but when I run it on the 32 bit file, I get RuntimeWarning: divide by zero encountered in log10. Obviously, there are 0's being passed to the log10 function to convert to dBFS. I'm wondering if there is a way I can avoid this? Why does this only happen to 32 bit int converted to 32 bit float? Is my method of converting to 32 bit float incorrect?

import numpy
import scipy.io.wavfil
import scipy.fftpack

def check_db_at_freq_channel(input_pcm_file, hz, channel):
    data_sample_rate, data_pcm = scipy.io.wavfile.read(input_pcm_file)

    data_len_rounded = round(data_pcm.shape[0] / float(data_sample_rate))

    if data_pcm.dtype == 'int16':
        max = float(2**15)
    elif data_pcm.dtype == 'int32':
        max = float(2**31)
    else:
        raise TypeError('data_pcm.dtype {} != {} or != {}'.format(data_pcm.dtype, 'int16', 'int32'))

    # Convert Integer data to floating point.
    data_pcm = data_pcm.astype(numpy.float32, order='C') / max

    # Check for a Mono Input
    if len(data_pcm.shape) == 1:
        data_pcm = data_pcm[:, None]

    data_fft = scipy.fftpack.fft(data_pcm[:, channel], n=int(data_len_rounded * data_sample_rate))

    # Convert to dBFS
    data_db = 20 * numpy.log10(2 * numpy.abs(data_fft) / float(len(data_pcm)))

    # Get Frequency Bins
    data_freqs = scipy.fftpack.fftfreq(data_fft.size, d=(1 / float(data_sample_rate)))

    # Get Index of requested Hz value.
    data_index = numpy.where(data_freqs == data_freqs[(numpy.abs(data_freqs - hz)).argmin()])[0][0]
    return data_db[data_index]

if __name__ == '__main__':
    test_file = '/path/to/some/file.wav'
    val = check_db_at_freq_channel(test_file, 440, 0)
    print val

I'm sure there is more information I could provide, so please let me know if I'm missing some crucial information?

Answer 1

I'm a newbie at Python, so there may be a slick way to do it I don't know about, but I would say your problem is really platform independent. Which of your test cases fail is not important, what matters is your code has the potential to fail. I would recommend that you put your dB calculation in a loop with a threshold if statement to determine whether you should take the log or assign a predefined value.

I am also wondering why you have "2 *" inside the log argument. This simply adds a log(2) to the log value. You already have "20 *" which is "2 * 10 *" to account for the squaring of the magnitude to calculate power.

Ced

Answer 2

Mathematically speaking, log(0) is undefined. The standard behavior for most computer languages, such as matlab, is to return -inf. Some languages, such as C++ will return -inf and raise a div by zero exception. Given that it's just a warning you could opt to ignore or suppress the warning and move on. Another option is to massage all of the 0 values to a very small positive number.

fft_mag = math.abs(fft_data)
fft_mag[fft_mag == 0] = np.finfo(float).eps