pyfilterbank icon indicating copy to clipboard operation
pyfilterbank copied to clipboard

Wrong dBA meassurements?

Open JonNeat opened this issue 4 years ago • 7 comments

Hi guys,

First thanks for writing this library. I'm using it for a sound project where I read the input from a microphone using the sounddevice library for python then giving the data to your library for dBa measurements.

The result is a bit off. When recording a quiet room using the microphone (samplerate 44100, dtype = int16), I get 30dBa. But when turning on a decibel-device (UNI-T UT353) it reads 40dBa.

Is my code wrong, or the UNI-T decibel-device?

Could it be that my ndarray contains int16, and not int24?

JonNeat avatar Aug 08 '19 15:08 JonNeat

Please post some example code. How did you calibrate? Did you use the data from the same microphone?

JonNeat [email protected] schrieb am Do., 8. Aug. 2019, 17:58:

Hi guys,

First thanks for writing this library. I'm using it for a sound project where I read the input from a microphone using the sounddevice library for python then giving the data to your library for dBa measurements.

The result is a bit off. The recordings from the microphone recording at a samplerate 44100, dtype = int16, I get 30dBa in a quiet room. But when turning on my UNI-T UT353 it reads 40dBa.

Is my code wrong, or the UNI-T decibel-device?

Could it be that my ndarray contains int16, and not int24?

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/SiggiGue/pyfilterbank/issues/17?email_source=notifications&email_token=AAXQ6HZZUWBTHXTV6TIBIG3QDQ7DPA5CNFSM4IKMDEL2YY3PNVWWK3TUL52HS4DFUVEXG43VMWVGG33NNVSW45C7NFSM4HEGHGDQ, or mute the thread https://github.com/notifications/unsubscribe-auth/AAXQ6HZKZDICF2NOUXJNFCDQDQ7DPANCNFSM4IKMDELQ .

SiggiGue avatar Aug 09 '19 07:08 SiggiGue

@SiggiGue

import numpy as np
import argparse
import sounddevice as sd
import datetime
import math
import sys
import queue
import splweighting as weighting #https://github.com/SiggiGue/pyfilterbank/blob/master/pyfilterbank/splweighting.py
import calendar;
import time;
np.set_printoptions(threshold=sys.maxsize)

def rms_flat(a):  # from matplotlib.mlab
    """
    Return the root mean square of all the elements of *a*, flattened out.
    """
    return np.sqrt(np.mean(np.absolute(a)**2))

q = queue.Queue()

def audio_callback(indata, frames, time, status):
    q.put(indata.copy())
    
with sd.InputStream(callback=audio_callback, blocksize=10000, samplerate=44100, dtype='int16'):
    f = open("output.txt", "a")
    itrSave = 50
    currentItr = 0
    while True:
        try:
            indata = q.get()
            weightedData = weighting.weight_signal(indata, 44100, 'A')
            dBa = 20*np.log10(rms_flat(weightedData))
            original = 20 * np.log10(rms_flat(indata))
            print('Original:   {:+.2f} dB'.format(original))
            print('A-weighted: {:+.2f} dB'.format(dBa))
            ts = calendar.timegm(time.gmtime())
            f.write(str(ts) + "," + str(dBa)+"\n")
            if currentItr == itrSave:
                f.flush()
                currentItr = 0
            
            currentItr += 1

        except KeyboardInterrupt:
            print("closing down")
            f.close()
            break

JonNeat avatar Aug 09 '19 12:08 JonNeat

The above code records from an attached microphone, then applies A-weighting and prints it to a file.

JonNeat avatar Aug 09 '19 12:08 JonNeat

If you want to measure the real sound pressure level A-weighted you must calibrate your audio recording system. The best practice is to do it using a certified sound pressure level calibrator, which emits a sine wave of 1 Pa RMS at 1 kHz at the microphone attached to it.

p_rms_calib = 1

You record the sine wave with your system gains set up to record the audio as close to the full scale amplitude without clipping, as you can, then get the RMS value of the recording:

myRecData = someRecFunc() p_rms_rec = RMS( myRecData )

The calibration factor of your recording system is

CF = p_rms_calib / p_rms_rec # CF for Calibration Factor

With CF you can adjust any of your recording data made with the exactly same hardware configuration used to calculate CF, so no more changing gains. If you do so, another CF must be obtained. The way to use it is by dividing your recorded data by it

myInterestData = someRecFunc() myInterestData_calibrated = myInterestData / CF

myInterestData_calibrated is now a representation of a real sound pressure wave, which can be weighted and used to calculate SPL.

The equation

mydBData = 20*log10(myData)

Calculates the myData level, which is also measured using decibel (dB). To calculate SPL, that is, decibel with reference to 20 uPa you need

mySPLdata = 20*log10(effective_sound_pressure / reference_sound_pressure)

The effective_sound_pressure is the RMS of myInterestData_calibrated. The reference_sound_pressure is 20 uPa, or 2e-05 Pa.

Hope it helps!

Chum4k3r avatar Aug 09 '19 15:08 Chum4k3r

@Chum4k3r thanks for the reply! Some of this is hard to understand without deep knowledge into sound. Can you apply it to my current code? If I understand you correct, my current code is OK, but just needs calibration?

I've simplified my above code a bit, so it's easier to read.

JonNeat avatar Aug 10 '19 09:08 JonNeat

Yes, your code is alright, it just needs the calibration step.

The calibration is something you only need to do once, so you can add a function that calibrate the system before the actual run.


import numpy as np
import argparse
import sounddevice as sd
import datetime
import math
import sys
import queue
import splweighting as weighting #https://github.com/SiggiGue/pyfilterbank/blob/master/pyfilterbank/splweighting.py
import calendar;
import time;
np.set_printoptions(threshold=sys.maxsize)

# Calibration function
def calibration(channelNumber, sampleRate):
    mySignal = sd.rec(int(5*sampleRate), samplerate=sampleRate,
                                 mapping=[channelNumber], dtype='float32', blocking=True)
    RMS = rms_flat(mySignal)
    return 1/RMS

def rms_flat(a):  # from matplotlib.mlab
    """
    Return the root mean square of all the elements of *a*, flattened out.
    """
    return np.sqrt(np.mean(np.absolute(a)**2))

q = queue.Queue()

def audio_callback(indata, frames, time, status):
    q.put(indata.copy())

# HERE YOU GET THE CALIBRATION VALUE
input("Place the calibrator on the microphone and press enter.")
CF = calibration(1, 44100)

input("Now take the calibrator off the microphone and press enter to start the real measurement.")
with sd.InputStream(callback=audio_callback, blocksize=10000, samplerate=44100, dtype='float32'):
    f = open("output.txt", "a")
    itrSave = 50
    currentItr = 0
    while True:
        try:
            indata = q.get()/CF  # HERE YOU APPLY IT
            weightedData = weighting.weight_signal(indata, 44100, 'A')
            dBa = 20*np.log10(rms_flat(weightedData) / 2e-5)  # referencing 2e-5 Pa
            original = 20 * np.log10(rms_flat(indata) / 2e-5)  # referencing 2e-5 Pa
            print('Original:   {:+.2f} dB'.format(original))
            print('A-weighted: {:+.2f} dB'.format(dBa))
            ts = calendar.timegm(time.gmtime())
            f.write(str(ts) + "," + str(dBa)+"\n")
            if currentItr == itrSave:
                f.flush()
                currentItr = 0
            
            currentItr += 1

        except KeyboardInterrupt:
            print("closing down")
            f.close()
            break

I've also changed the data type for a better data representation, since your output data also gets cast to float on the weighting steps, and added the level referencing on "decibel" calculation.

This is needed because "decibel" is just a scale, and anything can be scaled by it, so to have decibel that represents sound pressure level you must divide the scaled value (sound pressure) by the reference value.

Chum4k3r avatar Aug 10 '19 11:08 Chum4k3r

Thanks for the code @Chum4k3r . I'll see if I can find a calibrator :).

JonNeat avatar Aug 11 '19 10:08 JonNeat