Overall aim: to take a short signal (likely to be noisy with some tonal components) and filter it into third-octave bands (TOBs), apply frequency-dependent spectral adjustments, then recombine the signal, apply some other time-domain processing, then refilter again back into TOBs. Phase and amplitude to be retained as faithfully as possible throughout the process.
My approach: using Python language, generate TOB Butterworth IIR filter and apply forwards-backwards to ensure zero-phase. The cutoff frequencies for each filter band are adjusted to ensure the two-way application results in the same cutoff (3dB down) frequencies as the use of a one-way equivalent filter. The filtered signals would, in my application, be adjusted with band-dependent factors, however I've omitted this stage from my example, so I can see what the effects of the filtering+recombination process are. The bandpass filtered signals are then summed over the frequencies into a 'recombined' signal. They would then also be processed in some other ways, again omitted from my example for simplicity. Finally the recombined signal is again filtered using the same TOB two-way process, and the RMS amplitudes calculated for the total and TOB components of the signals at each stage.
Problem/question: I have noticed that the RMS amplitudes of my filtered+recombined signals are slightly higher than the signal entering the filter. I have tested two different types of complex signal as inputs: i) pink noise, and ii) a sum of sinusoids (with randomised phase shifts) at each TOB mid-frequency. The pink noise signal seems to get an increase of roughly ~1.2-1.4 dB for each application of the filtering+recombination process, while the sinusoid sum signal gets a smaller increase of ~0.3-0.6 dB. Why is the process adding energy into the signals? For most filtering processes, you might expect a small loss in total energy, but here I observe an increase.
Additional notes: Due to the random numbers used to generate the signals, the script produces slightly different outputs with every run. The factor used to adjust the cutoff frequencies for the forwards-backwards filter has been derived by experimentation - the appropriate factor to use has been found to be dependent on the ratio between the filter order and the sampling frequency, as well as the type of Butterworth design (bandpass/highpass/lowpass) - here it has been set according to the desired forwards-backwards filter order of 8 (ie the one-way filter order*2). The method of padding the filter has been applied on the basis of this article: Boore, DM, 2005. Previous investigation showed that using the default pad behaviour also showed an increase in the energy, so this parameter does not seem to be the main cause. I'm aware that my example signals are quite short (1 second) for the frequency range I'm interested in (5-630 Hz), but this is unavoidable with the signals I have to process.
import numpy as np
from scipy.signal import (butter, sosfiltfilt, sosfreqz)
import matplotlib.pyplot as plt
from random import random
# pink noise generation
# credit to python-acoustics library by Python Acoustics https://github.com/python-acoustics/python-acoustics
N = 8192 # signal length (samples), and sampling frequency
state = np.random.RandomState()
uneven = N % 2
X = state.randn(N // 2 + 1 + uneven) + 1j * state.randn(N // 2 + 1 + uneven)
S = np.sqrt(np.arange(len(X)) + 1.) # +1 to avoid divide by zero
pink = (np.fft.irfft(X / S)).real
if uneven:
pink = pink[:-1]
T = N/8192 # signal duration, 1 s
t = np.linspace(0, T-(T/N), N) # signal time vector
# define filter order
order = 4
fbw_order = order*2 # fowards-backwards filter order
# generate third-octave band frequencies
b = 3
ind = np.arange(-23, -1, 1) # range of frequency indices
G10 = 10**(3/10) # octave ratio coefficient (base-ten)
OctRatio = G10**(0.5/b) # octave ratio
fm = G10**(ind/b)*1000 # output range of exact fractional mid-frequencies
f1 = fm/OctRatio # output range of exact lower band-edge frequencies
f2 = fm*OctRatio # output range of exact upper band-edge frequencies
# generate sine-based complex signal
sine = np.zeros(8192)
for f in fm:
sine += (0.003/fm.size)*np.sin(2*np.pi*f*t + random()*np.pi)
# forwards-backwards filter pre-warped frequencies
f1_fbw_HP = f1/1.12
f2_fbw_LP = f2*1.11
f1_fbw_BP = f1/1.0135
f2_fbw_BP = f2*1.0135
# generate forwards-backwards HP and LP filters
sos_HP = butter(order, f1_fbw_HP[0], btype='highpass', output='sos', fs=N)
sos_LP = butter(order, f2_fbw_LP[-1], btype='lowpass', output='sos', fs=N)
# filter pink noise with LP and HP filters
padN = np.min([N-1, int(N*1.5*order/f2_fbw_LP[-1])])
psignal = sosfiltfilt(sos_LP, pink, padtype='constant', padlen=padN)
padN = np.min([N-1, int(N*1.5*order/f1_fbw_HP[0])])
psignal = sosfiltfilt(sos_HP, psignal, padtype='constant', padlen=padN)
# calculate RMS and dB re 1e-9
psignal_RMS = np.sqrt(np.mean(np.square(psignal)))
psignal_RMS_dB = 20*np.log10(np.sqrt(np.mean(np.square(psignal)))/1e-9)
ssignal_RMS = np.sqrt(np.mean(np.square(sine)))
ssignal_RMS_dB = 20*np.log10(np.sqrt(np.mean(np.square(sine)))/1e-9)
# loop to filter signals into third-octave bands
psignal_TOB_fbw = np.zeros((psignal.size, fm.size))
ssignal_TOB_fbw = np.zeros((sine.size, fm.size))
for ii, f in enumerate(fm):
sos_fbw_BP = butter(order, [f1_fbw_BP[ii], f2_fbw_BP[ii]], btype='bandpass',
output='sos', fs=N)
padN = np.min([N-1, int(N*1.5*order/f2_fbw_BP[ii])])
psignal_TOB_fbw[:, ii] = sosfiltfilt(sos_fbw_BP, pink, padtype='constant',
padlen=padN)
ssignal_TOB_fbw[:, ii] = sosfiltfilt(sos_fbw_BP, sine, padtype='constant',
padlen=padN)
# calculate TOB root-mean-square and dB re 1e-9
psignal_TOB_fbw_RMS = np.sqrt(np.mean(np.square(psignal_TOB_fbw), axis=0))
psignal_TOB_fbw_RMS_dB = 20*np.log10(psignal_TOB_fbw_RMS/1e-9)
ssignal_TOB_fbw_RMS = np.sqrt(np.mean(np.square(ssignal_TOB_fbw), axis=0))
ssignal_TOB_fbw_RMS_dB = 20*np.log10(ssignal_TOB_fbw_RMS/1e-9)
# recombine TOBs and calculate total RMS and dB re 1e-9
psignal_fbw_re = np.sum(psignal_TOB_fbw, axis=1)
psignal_fbw_RMS = np.sqrt(np.mean(np.square(psignal_fbw_re)))
psignal_fbw_RMS_dB = 20*np.log10(psignal_fbw_RMS/1e-9)
ssignal_fbw_re = np.sum(ssignal_TOB_fbw, axis=1)
ssignal_fbw_RMS = np.sqrt(np.mean(np.square(ssignal_fbw_re)))
ssignal_fbw_RMS_dB = 20*np.log10(ssignal_fbw_RMS/1e-9)
# loop to refilter recombined signal back into TOBs
psignal_TOB_fbw2x = np.zeros((psignal_fbw_re.size, fm.size))
ssignal_TOB_fbw2x = np.zeros((ssignal_fbw_re.size, fm.size))
for ii, f in enumerate(fm):
sos_fbw_BP = butter(order, [f1_fbw_BP[ii], f2_fbw_BP[ii]], btype='bandpass',
output='sos', fs=N)
padN = np.min([N-1, int(N*1.5*order/f2_fbw_BP[ii])])
psignal_TOB_fbw2x[:, ii] = sosfiltfilt(sos_fbw_BP, psignal_fbw_re,
padtype='constant', padlen=padN)
ssignal_TOB_fbw2x[:, ii] = sosfiltfilt(sos_fbw_BP, ssignal_fbw_re,
padtype='constant', padlen=padN)
# calculate TOB root-mean-square and dB re 1e-9
psignal_TOB_fbw2x_RMS = np.sqrt(np.mean(np.square(psignal_TOB_fbw2x), axis=0))
psignal_TOB_fbw2x_RMS_dB = 20*np.log10(psignal_TOB_fbw2x_RMS/1e-9)
ssignal_TOB_fbw2x_RMS = np.sqrt(np.mean(np.square(ssignal_TOB_fbw2x), axis=0))
ssignal_TOB_fbw2x_RMS_dB = 20*np.log10(ssignal_TOB_fbw2x_RMS/1e-9)
# recombine TOBs and calculate total RMS and dB re 1e-9
psignal_fbw2x_re = np.sum(psignal_TOB_fbw2x, axis=1)
psignal_fbw2x_RMS = np.sqrt(np.mean(np.square(psignal_fbw2x_re)))
psignal_fbw2x_RMS_dB = 20*np.log10(psignal_fbw2x_RMS/1e-9)
ssignal_fbw2x_re = np.sum(ssignal_TOB_fbw2x, axis=1)
ssignal_fbw2x_RMS = np.sqrt(np.mean(np.square(ssignal_fbw2x_re)))
ssignal_fbw2x_RMS_dB = 20*np.log10(ssignal_fbw2x_RMS/1e-9)
# calculate and display level differences
print(psignal_fbw_RMS_dB - psignal_RMS_dB)
print(psignal_fbw2x_RMS_dB - psignal_fbw_RMS_dB)
print(psignal_TOB_fbw2x_RMS_dB - psignal_TOB_fbw_RMS_dB)
print(ssignal_fbw_RMS_dB - ssignal_RMS_dB)
print(ssignal_fbw2x_RMS_dB - ssignal_fbw_RMS_dB)
print(ssignal_TOB_fbw2x_RMS_dB - ssignal_TOB_fbw_RMS_dB)
