Also, I think his is an unstable filter.

I won't have time to look at this in the coming two weeks, but please, run some tests. Maybe create a Notebook that shows what's going on? We should definitely use second-order sections (Scipy didn't support that when I wrote this code).

Note that I haven't really used fast weighting, instead I use averaging.

Hello,
I've noticed the same thing as A-Mendi for pole of the filter

Also, could you please confirm if it's normal to split the signal into chunks and extract the value before the last for each filtered chunks (code below from function integrate):

data = data.reshape(newshape)
return lfilter(b, a, data)[...,n-1] / integration_time

The standard suggests to apply the filter on the full time signal. Thus, one would expect to obtain a time resolution equal to the input. With the current approach, the time resolution is given by the integration time.

@FRidh Please, let me know what you think.
Thank you!

Hi. It's been a long time since I last looked at this code. I'm not working in the acoustics field anymore either and no longer have access to the standards, so I won't be much of help here I'm afraid.

I think that is eventually up to you and other that want to contribute. I'm fine with reviewing/merging pull requests and at some point giving merge rights. I was actually thinking of automating pushing to PyPI whenever a tag is added so we could have rapid releases. But I don't have the time to go into the technical details of methods and standards.

Thanks @topherbuckley for suggesting to help with that!

I am not very familiar with the implementation of code modifications through GitHub, but I can help with the technical details. I also have access to standards.

Alright. So I had a moment to look at this closer finally. I'm not really sure what was trying to be achieved by the current filter. It has a zero at s=1, and two poles; one at 1.0 and another at integration_time. I just implemented this exponential weighting function in some other code I'm working on, and I'll be the first to admit that the standards are VERY confusing on this topic. I read up on some existing implementations using this as a reference. I highly recommend it if you are going to be diving deeper into these standards and implementing them.

Anyways, in that reference book, a digital filter is already designed (eq 7.9 if you end up getting it) and so the code for the filter can be produced with:

b= np.poly1d([1, 0])
a= np.array([(1 + (sample_frequency * integration_time)), -sample_frequency * integration_time])

This would replace lines 160-161 in the current code. This could also be done using an analog filter first, then using the same line on 161 to convert it to a digital filter. Therefore line 160 could also just be replaced with the following two lines while keeping the bilinear transform directly after:

b = np.poly1d([1])
a = np.poly1d([integration_time, 1])

The above would represent the "low-pass filter with one real pole at -1/tau" as stated in the IS standard referenced by @A-Mendi. Since the 2nd implementation better matches the wording of the standards, I'd say thats the way to go, even though its a bit more computationally expensive.

As far as what @duracell013 is talking about regarding the filtered chunks, yes you are right that it should theoretically be applied to the full signal, but in practice this would lead to overflow errors for very large signals (imagine turning this on for hours/days and filling up a VERY large array). So what is done in practice is using a buffer of a particular length that would result in approximately the same value as what would be given by the full signal. In other words, even if the signal was one minute long, and fast-weighted, you'd really only need about a half a second in your buffer array to get within 0.1 dB of the value you'd get from the full one minute array of data. This is simply because the exponential weighting is doing what its supposed to do, i.e. weighting all past data near zero and all recent data on an exponential slope. I haven't seen this written anywhere in the standard or elsewhere, but from testing on steady state signals I've found this to be around 0.53 s for the fast time constant and 4.5 s for the slow time constant if you want to be within 0.1 dB.

I'm calling it quits for tonight, but I can implement the above sometime tomorrow perhaps. I have similar code I've been working on recently in another library.

A follow up question for @FRidh though; what is the purpose of all the reshaping or arrays and defining integration_time and sample_frequency as arrays? It adds several lines to the code, and makes it quite cryptic in my opinion. Is there some significant speed advantage to this or were you preparing this for multidimensional array processing for some reason? The use of the Ellipses for array slicing was completely new to me, and seems to be fairly rare according to this. Can I rewrite this, at least in this section in question, to be clearer using standard slicing removing all the excessive array modifications?

NIce work!

A follow up question for @FRidh though; what is the purpose of all the reshaping or arrays and defining integration_time and sample_frequency as arrays? It adds several lines to the code, and makes it quite cryptic in my opinion. Is there some significant speed advantage to this or were you preparing this for multidimensional array processing for some reason?

Yep, vectorization in case multiple signals are fed. This is used in multiple places and gives a significant performance boost in case of multi-channel signals. I typically use this code with the Signal class, which can also handle multiple channels. Therefore, I suggest not breaking that.

@FRidh Understood regarding the multi-channel signals. I created a pull request for just the simple filter update, but still have some followup regarding the integration time buffer implementation. As I'm not that familiar with your code yet, can you give me a better idea of the integrate() method is used elsewhere? Is it used elsewhere other than in this iec_61672_1_2013.py or the corresponsind test_iec_61672_1_2013.py?

How I have it set up in my other library is to be used as a constantly running while-loop since the main goal of the program was to act as a sound level meter. This kind of behavior obviously could cause problems if you expect it to be used simultaneously with something else. I also made a multi-threaded version of this but that would take a fair bit of time to dig back up out of old versions so I'd prefer not to. So if we can make the single while loop work for buffering, then that would be easy for me to set up, but I'd want your input as to whether this should be added as a new method, or within the integrate method or elsewhere.

My current method uses PyAudio as the data stream, which would also be a question for you. Are you using PyAudio, or is there some other streaming of data going on somewhere else that would feed the "data" input parameter for your integrate() method? We might approach this by giving instructions to start a data stream, either via PyAudio or another library, prior to starting the exponential integration loop, since PyAudio (as I'm sure others) start a separate thread to handle the audio streaming which doesn't seem to troublesome in conjunction with my exponential filtering while-loop.

To give you an idea of what it takes, I've trimmed down some of my own code and pasted below. Sorry, I haven't had the time to make sure variable names agree with your own and whatnot, but it should be close enough to understand. Also, some of the architecture of the code would conflict with your own right now, but that can be worked out. Let me know if you have any questions and I can elaborate on anything.

# System Imports
import datetime
import time
import sys

# Third-Party Imports
import pyaudio
import numpy as np
from numpy import log10, sqrt, mean
from scipy.signal import lfilter, bilinear

Fs = 44100
tau = 1.0
z_weighted_level = 0.0
stop_bool = False
FORMAT = pyaudio.paFloat32
CHANNELS = 2
# This is hardcoded at the moment to represent the USB audio card.
# Use pyaudio method get_host_api_info_by_index(i) to find appropriate index
INPUT_DEVICE_INDEX = 0


def start_exponential_averaging():
    """
    4.5s determined to be 0.1 db difference in exp averaging, so going with 5
    This should be updated if the time constant is changed.

    buffer_length_total represents the total kept buffer length,
    i.e. all data after 5 seconds is thrown out.
    buffer_length_single is for the purposes of keeping the refresh rate
    of the results some time below the total buffer length.

    """

    buffer_length_total = 5.0
    buffer_length_single = 1.0
    CHUNK = Fs * int(round(buffer_length_single))
    p = pyaudio.PyAudio()
    stream = p.open(input_device_index=INPUT_DEVICE_INDEX,
                    format=FORMAT,
                    channels=CHANNELS,
                    rate=Fs,
                    input=True,
                    frames_per_buffer=CHUNK)

    exponential_averaging(Fs, stream, buffer_length_single,
                          buffer_length_total, CHUNK)

    # Wait for above while loop averaging to be externally stopped via
    # the stop_bool then close all streams and PyAudio objects
    stream.stop_stream()
    stream.close()
    p.terminate()


def exponential_averaging(Fs, stream, buffer_length_single,
                          buffer_length_total, CHUNK):
    global stop_bool
    global z_weighted_level

    data = np.array([])
    start_time = time.time()
    delay = 0
    while not stop_bool:
        ideal_end_time = start_time + buffer_length_single
        block = stream.read(CHUNK - delay)
        if len(data) >= Fs * buffer_length_total:
            data = np.array([])
        data = np.append(data, np.fromstring(block, dtype=np.float32))

        z_weighted_level = z_weighted_level_calc(data, 1.0)
        # a_weighted and c_weighted can go here as well

        # Wait until integration time is over before starting a new calculation
        # This is only applicable when multithreading and the calculation
        # thread finishes before the start of the next stream chunk.
        comparer_start = ideal_end_time - time.time()
        comparer = comparer_start
        while comparer >= 0:
            comparer = ideal_end_time - time.time()

        calculation_time = time.time() - start_time
        delay = int(round((calculation_time - buffer_length_single) * Fs))
        start_time = ideal_end_time
        ideal_end_time = ideal_end_time + buffer_length_single


def z_weighted_level_calc(self, data, avg_time):
    # Check for overflow before squaring
    if float(np.max(data))**2 > sys.float_info.max:
        print('float32 dtype resulting in overflow at z_weighted_level_calc')
    z_weighted_data_squared = data**2
    # Exponential Averaging
    b_analog = np.poly1d([1])
    a_analog = np.poly1d([tau, 1])
    b, a = bilinear(b_analog, a_analog, fs=Fs)
    z_weighted_exp_avg_data = lfilter(
        b, a, z_weighted_data_squared)
    # Square Root and Logarithm
    z_weighted_level = 20 * \
        log10(sqrt(mean(z_weighted_exp_avg_data)))
    # print('{:.1f} dBZ'.format(z_weighted_level))
    return z_weighted_level

This is simply because the exponential weighting is doing what its supposed to do, i.e. weighting all past data near zero and all recent data on an exponential slope. I haven't seen this written anywhere in the standard or elsewhere, but from testing on steady state signals I've found this to be around 0.53 s for the fast time constant and 4.5 s for the slow time constant if you want to be within 0.1 dB.

@topherbuckley I understand the limitation of processing very long signals. However, my understanding of the current implementation means that we only process a very small segment at a time (duration = time constant), which might be too short to obtain a correct level. Am I wrong?

In my applications I need to calculate the SPL with a very small time constant (10 ms). Do you think this approach is still valid for such duration?

@duracell013 , yes you are right. A 10 ms time constant would give noisy enough results even if done properly. Averaging over only 10 ms would likely lead to extremely noisy results, and no, it would not properly capture what the standards seek to measure.

The code I posted above is set for a time constant of 1 s, and avoids the problem you speak of by taking "CHUNKS" in a size of 1 second (in this case, but this can be changed via the buffer_length_total and buffer_length_single variables) and then adds them to a 5 second buffer continuously, i.e. from 5-6 seconds, the first recorded CHUNK is dropped and the new CHUNK appended to the end. 5 seconds is chosen in order to get within 0.1 dB of the long term average for a steady-state signal. For your case, of a 10 ms time constant, this should still work, but after adjusting these buffer lengths. I'm not sure how long the buffer would have to be for such a small time constant in order to get to within 0.1 dB of the full signal average for steady-state signals, but based on my results with the 0.125 s and 1 s constants, it looks like its around a factor 4.5 - 5 longer than the time constant itself, so I'd recommend buffering for at least that length (i.e. 0.010 s * 4.5).

I'm not sure of your end goal here, but if for example you wanted to have a text file with recorded sound levels each taken x seconds after the first, then you would simply set buffer_length_total = 4.5 * 0.010 , and buffer_length_single = x, noting buffer_length_total > buffer_length_single. You would then either make another array to store your values and print to file later, or print to file at the end of each loop in the "while not stop_bool:" loop. For such a small time constant, I'd recommend the former, otherwise the time spent writing to file will likely exceed your time constant, and you'll have a lot of skipped data. I've timed this code, and each loop takes anywhere from 5 to 30 ms, so careful of that overlap in your case. You may have to consider a multi-threaded approach if you want to use such a small time constant, in order to avoid skipping any data. For my pursposes, with larger time constants, especially with the 1 s, I'm only missing a small percentage of each CHUNK, you'd be missing entire CHUNKs.

Side note: I messed the buffering mentioned above now that I reviewed it. See below for update. Just changed one line data = np.array([]) --> data = data[len(block):], which would take the data buffers most recent 4 blocks only, the proceed to add the newest to the end in the next lines.

# System Imports
import datetime
import time
import sys

# Third-Party Imports
import pyaudio
import numpy as np
from numpy import log10, sqrt, mean
from scipy.signal import lfilter, bilinear

Fs = 44100
tau = 1.0
z_weighted_level = 0.0
stop_bool = False
FORMAT = pyaudio.paFloat32
CHANNELS = 2
# This is hardcoded at the moment to represent the USB audio card.
# Use pyaudio method get_host_api_info_by_index(i) to find appropriate index
INPUT_DEVICE_INDEX = 0


def start_exponential_averaging():
    """
    4.5s determined to be 0.1 db difference in exp averaging, so going with 5
    This should be updated if the time constant is changed.

    buffer_length_total represents the total kept buffer length,
    i.e. all data after 5 seconds is thrown out.
    buffer_length_single is for the purposes of keeping the refresh rate
    of the results some time below the total buffer length.

    """

    buffer_length_total = 5.0
    buffer_length_single = 1.0
    CHUNK = Fs * int(round(buffer_length_single))
    p = pyaudio.PyAudio()
    stream = p.open(input_device_index=INPUT_DEVICE_INDEX,
                    format=FORMAT,
                    channels=CHANNELS,
                    rate=Fs,
                    input=True,
                    frames_per_buffer=CHUNK)

    exponential_averaging(Fs, stream, buffer_length_single,
                          buffer_length_total, CHUNK)

    # Wait for above while loop averaging to be externally stopped via
    # the stop_bool then close all streams and PyAudio objects
    stream.stop_stream()
    stream.close()
    p.terminate()


def exponential_averaging(Fs, stream, buffer_length_single,
                          buffer_length_total, CHUNK):
    global stop_bool
    global z_weighted_level

    data = np.array([])
    start_time = time.time()
    delay = 0
    while not stop_bool:
        ideal_end_time = start_time + buffer_length_single
        block = stream.read(CHUNK - delay)
        if len(data) >= Fs * buffer_length_total:
            data = data[len(block):]
        data = np.append(data, np.fromstring(block, dtype=np.float32))

        z_weighted_level = z_weighted_level_calc(data, 1.0)
        # a_weighted and c_weighted can go here as well

        # Wait until integration time is over before starting a new calculation
        # This is only applicable when multithreading and the calculation
        # thread finishes before the start of the next stream chunk.
        comparer_start = ideal_end_time - time.time()
        comparer = comparer_start
        while comparer >= 0:
            comparer = ideal_end_time - time.time()

        calculation_time = time.time() - start_time
        delay = int(round((calculation_time - buffer_length_single) * Fs))
        start_time = ideal_end_time
        ideal_end_time = ideal_end_time + buffer_length_single


def z_weighted_level_calc(self, data, avg_time):
    # Check for overflow before squaring
    if float(np.max(data))**2 > sys.float_info.max:
        print('float32 dtype resulting in overflow at z_weighted_level_calc')
    z_weighted_data_squared = data**2
    # Exponential Averaging
    b_analog = np.poly1d([1])
    a_analog = np.poly1d([tau, 1])
    b, a = bilinear(b_analog, a_analog, fs=Fs)
    z_weighted_exp_avg_data = lfilter(
        b, a, z_weighted_data_squared)
    # Square Root and Logarithm
    z_weighted_level = 20 * \
        log10(sqrt(mean(z_weighted_exp_avg_data)))
    # print('{:.1f} dBZ'.format(z_weighted_level))
    return z_weighted_level

If on the other hand, you are not seeking to do any real-time measurements, but instead looking only to get the filtered version of a time-definite signal, you could keep the existing code and just delete the following lines

    n = np.floor(integration_time * sample_frequency).astype(int)
    data = data[..., 0:n*(samples//n)]
    newshape = list(data.shape[0:-1])
    newshape.extend([-1, n])
    data = data.reshape(newshape)
    #data = data.reshape((-1, n)) # Divide in chunks over which to perform the integration.

and modify the next line from

return lfilter(b, a, data)[...,n-1] / integration_time # Perform the integration. Select the final value of the integration.

to

return lfilter(b, a, data) # Perform the integration. Select the final value of the integration.

which would filter the full data array signal.

One thing I'd appreciate your input on, is whether you think normalizing by the integration_time, or the data length (as was done in the existing code) is necessary? I haven't gotten my code from the other library to the end goal of showing dBA and dBC quite yet, so I'm not sure if my gains/normalizations are set appropriately yet or not. Normalizing it like this had not occurred to me as necessary, but I could be wrong.