Infrasound Report with G Weighting

Note: I have updated this and the code since solving the dB(G) estimation issues.

This is a project that has been on my list for a while, but has become a little more urgent lately with proposed new wind farms in our local government area. I started getting serious a couple of weeks ago looking for a signal processing package for Python to do the job. I didn’t find any that worked with G weighting (which is the appropriate weighting to use for infrasound) so started looking at modifying or adapting one for AB and C weighting to suit. The maths and/or coding for that task beat me, but in the last couple of days I had an idea for a workaround. This is it.

I started with some test code and generated some pure sine wave test signals to make sure that when I normalised the FFT, it was to the same units as the Boom i.e. Pascals.

From there I produced an FFT in dBL (decibels Linear).

I did some curve fitting in MS Excel to come up with equations for the G weighting, and used these to convert the dBL FFT to a dB(G) FFT.

Now, getting the G weighting back to the time domain from this point is not possible (as far as I know). In a real SPL meter, the weighting circuitry operates live between the microphone and the meter, so changes in the frequency distribution are taken care of, We can’t do that this way, but if we assume that the frequency distribution doesn’t change (i.e. is averaged over the sample period) we should be able to get a reasonable approximation - especially for short sample periods.

To get around this issue of not being able to apply the G weighting in real time, the difference between the linear FFT and the db(G) FFT is averaged across the bandpass filter range. This figure is then used as an estimate of the difference between the SPL in dBL and the SPL in dB(G).

Below are some test plots from the report, which as far as I can find at the moment, it works as it should.


Initial test report, all appears to work as expected. The estimated Peak dB(G) in the time domain is reported in the bottom right corner with the actual peak dBL.

A short duration sample, also seem consistent.

A narrow(er) band test around 10Hz looks as it should as well, and the estimated Peak dB(G) is also behaving as it should. At 10 Hz the G weighting is the same as linear and that shows perfectly.

My initial attempt to estimate the difference between dBL and dB(G) was a result of over thinking the problem. Yes, subtracting dB is the same as division when it is converted back, but it’s valid when using dB!


# -*- coding: utf-8 -*-
"""
Created on Mon Jun 24 08:19:54 2024

@author: al72
"""
from obspy.clients.fdsn import Client
from obspy import UTCDateTime
from matplotlib.ticker import AutoMinorLocator
import numpy as np
import matplotlib.pyplot as plt

# define fdsn client to get data from
client = Client('https://data.raspberryshake.org/')

# define start and end times
eventTime = UTCDateTime(2024, 6, 19, 1, 13, 51) # (YYYY, m, d, H, M, S) **** Enter data****
delay = 0               #delay from event time to start of plot
start = eventTime + delay   #calculate the plot start time
duration = 40               #duration of plot in seconds
end = start + duration                # start plus plot duration in seconds (recommend a minimum of 10s)
daylightSavings = False

# Name the Event
eventName = 'Quarry Blast at Oberon Quarries'       # Name for the Report

# get data from the FSDN server and detrend it
STATION = "R21C0"               # station name
st = client.get_waveforms("AM", STATION, "00", "HDF", starttime=start, endtime=end, attach_response=True)
st.merge(method=0, fill_value='latest')
st.detrend(type="demean")
cmax = abs(st[0].max())     # maximum amplitude of raw counts

# get Instrument Response
inv = client.get_stations(network="AM", station=STATION, level="RESP")

# set up bandpass filter
flc = 0.5      #enter bandpass filter lower corner frequency
fuc = 49        #enter bandpass filter upper corner frequency
filt = [flc-0.01, flc, fuc, fuc+0.1]   

# remove instrument response now that raw trace has been processed
y = st.remove_response(inventory=inv,pre_filt=filt,output="DEF", water_level=60)

fs = 100    # sampling frequency
tstep = 1/fs    # sample time interval
f0 = 10     #signal frequency
ns = len(y[0].data)   #number of samples

t = np.linspace(0, (ns-1)*tstep, ns)    # time steps
fstep = fs / ns     # frequency interval
f = np.linspace(0, (ns-1)*fstep, ns)     #frequency steps

# FFT
yfft = np.fft.fft(y[0].data)
y_mag = np.abs(yfft)/ns     #normalise the FFT

f_plot = f[0:int(ns/2+1)]
y_mag_plot = 2*y_mag[0:int(ns/2+1)]
y_mag_plot[0] = y_mag_plot[0]/2     # do not multiply DC value

#calculate maximum sound pressure level
pmax = abs(y[0].max())   #find the max pressure amplitude
splmax = 10*np.log10(pmax*pmax) + 93.979400087        #convert to sound pressure level

# Convert FFT to dB(L)
ydB = []
FFTsum = 0
for i in range (0, len(y_mag_plot.data)):
    ydB.append(20*np.log10(abs(y_mag_plot[i])/0.00002))

# Convert DB(L) to dB(G)
ydBG = []
for i in range (0, len(ydB)):
    if f_plot[i] == 0:
        ydBG.append(ydB[i+1]+34.142*np.log(f_plot[i+1])-41.332)
    elif f_plot[i] < 1:
        ydBG.append(ydB[i]+34.142*np.log(f_plot[i])-41.332)
    elif f_plot[i] <= 3.15:
        ydBG.append(ydB[i]-1.7406*f_plot[i]**4+15.961*f_plot[i]**3-55.565*f_plot[i]**2+95.938*f_plot[i]-97.475)
    elif f_plot[i] < 12.5:
        ydBG.append(ydB[i]+17.396*np.log(f_plot[i])-40.037)
    elif f_plot[i] <= 20:
        ydBG.append(ydB[i]-.0976*f_plot[i]**2+3.8393*f_plot[i]-28.738)
    elif f_plot[i] < 31.5:
        ydBG.append(ydB[i]-.0108*f_plot[i]**2-.5724*f_plot[i]+24.782)
    else:
        ydBG.append(ydB[i]+.0076*f_plot[i]**2-1.4868*f_plot[i]+35.262)

# convert db(G) back to G weighted FFT in Pa
fftG = []
l2G = 0
j = 0
for i in range (0, len(ydB)):
    fftG.append(0.00002*10**(ydBG[i]/20))
    if f_plot[i] >= flc and f_plot[i] < fuc:
        l2G += ydB[i] - ydBG[i]
        j += 1

# estimated Linear to G weighted overall correction
l2G = l2G/j

# plot charts
# set-up figure and subplots
fig = plt.figure(figsize=(18,12), dpi=150)    #18 x 12 inches
ax1 = fig.add_subplot(4,1,1)    #left top 2 high
ax2 = fig.add_subplot(4,1,2)    #left middle 2 high
ax3 = fig.add_subplot(4,1,3)    #right top 2 high
ax4 = fig.add_subplot(4,1,4)    #right middle 2 high
ax1.plot(y[0].times(reftime=start), y[0].data, lw=1, color='k')
ax1.set_ylabel('Infrasound, Pa')
ax1.xaxis.set_minor_locator(AutoMinorLocator(10))
ax1.yaxis.set_minor_locator(AutoMinorLocator(5))
ax1.margins(x=0)
ax2.plot(f_plot, y_mag_plot, lw=1, color = 'b', label='linear')
ax2.plot(f_plot, fftG, lw=1, color = 'r', label='G weighted')
ax2.set_ylabel('FFT, Pa')
ax2.legend()
ax2.xaxis.set_minor_locator(AutoMinorLocator(10))
ax2.yaxis.set_minor_locator(AutoMinorLocator(5))
#ax2.set_yscale('log')
ax2.margins(x=0)
ax3.plot(f_plot, ydB, lw=1, color = 'g', label='dBL')
ax3.set_ylabel('FFT, dB')
ax3.plot(f_plot, ydBG, lw=1, color='r', label='dB(G)')
ax3.legend()
ax3.xaxis.set_minor_locator(AutoMinorLocator(10))
ax3.yaxis.set_minor_locator(AutoMinorLocator(5))
#ax3.set_yscale('log')
ax3.margins(x=0)
ax4.plot(y[0].times(reftime=start), 20*np.log10(abs(y[0].data)/0.00002), lw=1, color='purple')
ax4.set_ylabel('Infrasound, dBL')
ax4.xaxis.set_minor_locator(AutoMinorLocator(10))
ax4.yaxis.set_minor_locator(AutoMinorLocator(5))
ax4.margins(x=0)
ax1.grid(color='dimgray', ls = '-.', lw = 0.33)
ax2.grid(color='dimgray', ls = '-.', lw = 0.33)
ax3.grid(color='dimgray', ls = '-.', lw = 0.33)
ax4.grid(color='dimgray', ls = '-.', lw = 0.33)
ax1.grid(color='dimgray', which='minor', ls = ':', lw = 0.33)
ax2.grid(color='dimgray', which='minor', ls = ':', lw = 0.33)
ax3.grid(color='dimgray', which='minor', ls = ':', lw = 0.33)
ax4.grid(color='dimgray', which='minor', ls = ':', lw = 0.33)

# add text to report
if daylightSavings:
    dsc = 11
    dst = 'AEDT'
else:
    dsc=10 
    dst = 'AEST'
    
fig.suptitle('Infrasound Report: '+eventName+' - '+eventTime.strftime('%d/%m/%Y %H:%M:%S.%f UTC'), size='xx-large',color='b')
fig.text(0.12, 0.945, 'Start: '+start.strftime('%d/%m/%Y %H:%M:%S UTC')+'   '+(start+dsc*3600).strftime('%d/%m/%Y %H:%M:%S '+dst))
fig.text(0.12, 0.93, 'End: '+end.strftime('%d/%m/%Y %H:%M:%S UTC')+'   '+(end+dsc*3600).strftime('%d/%m/%Y %H:%M:%S '+dst))
fig.text(0.12, .915, 'Station: AM.'+STATION+'.00.HDF')
fig.text(0.12, 0.9, 'Bandpass Filter: '+str(flc)+' to '+str(fuc)+' Hz.')
fig.text(0.88, 0.73, 'Peak Pressure = '+str(round(pmax,3))+'Pa. ('+str(filt[1])+" to "+str(filt[2])+"Hz)", ha='right')    #report the peak Pressure
fig.text(0.88, 0.14, 'Peak Infrasound Pressure Level ='+str(round(splmax,1))+' dBL. ('+str(filt[1])+" to "+str(filt[2])+"Hz)", ha='right')   #report peak sound pressure level
fig.text(0.88, 0.155, 'Estimated Peak Infrasound Pressure Level ='+str(round(splmax-l2G,1))+' dB(G). ('+str(filt[1])+" to "+str(filt[2])+"Hz)", ha='right')   #report peak sound pressure level
fig.text(0.9, 0.945, 'Raspberry Shake and Boom', size='small', ha='right', color='r')
fig.text(0.9, 0.935, '@AlanSheehan18', size='small', ha='right')
fig.text(0.9, 0.925, '@raspishake', size='small', ha='right')
fig.text(0.9, 0.915, 'Oberon Citizen Science Network (OCSN)', size='small', ha='right')
fig.text(0.9, 0.905, '#CitizenScience', size='small', ha='right')

Al.

2 Likes

I did some testing with octave filters at different (e.g. 1 octave centred on 2 and 10 Hz respectively in the examples posted. These were checked against the expected difference in dB due to the G weighting and they are in good agreement.

Al.

2 Likes

As always, this is an amazing endeavor sheeny, with a likewise presentation!

I’m afraid that I cannot personally offer much more support in verifying your code, but I’m bumping this topic up again so it becomes more visible.

1 Like

Thanks Stormchaser,

That’s OK, I have the code sorted. It’s a work in progress…

I’m currently rewriting it as a module. That should simplify the code in the calling program.

Currently I’ve got a crude module that needs some tidying up and documenting (hiding some of my bad coding habits!), but I want to add more functions for standard octave and 1/3 Octave filter frequencies. Once complete, we should be able to produce reports compatible with any Sound Pressure Level meter reports or compliance condition.

The advantage of the Shake and Boom is of course it monitors continuously, so if conditions change to cause a non-compliance, it can be captured.

Once complete I’ll share it on Github.

Once the module is complete it shouldn’t be too hard to split a sample into sub-samples, calculate the G weighting, Peak, average and RMS amplitudes for each sub-sample, average the FFTs from the sub-samples, and reconstruct the sample trace and sound pressure level plots. That should go a long way to accounting for changes in frequency distribution within the sample (especially if they are long).

Likewise with the octave and 1/3 octave filters, those reports can be generated for direct comparison as well.

Here’s some of the latest test reports with RMS and Average Pressures (Amplitudes) added. The effect on average and RMS amplitude for a transient signal and change in sample length is well demonstrated in these plots.



1 Like