0
$\begingroup$

I am very new to signal processing and coding which is why my questions might be really basic. I have a signal of the acoustic pressure p'(t) and I would like to use a third octave bandpass filter in Python. I came across these two approaches to filtering with scipy:

Bandpass filters with python for low frequencies
https://scipy-cookbook.readthedocs.io/items/ButterworthBandpass.html

At the moment I am able to use one bandpass filter with the following code for a simple example (mostly taken from link 2):

import scipy
import numpy as np
from matplotlib import pyplot as plt
from scipy.signal import butter, lfilter
from scipy.signal import freqz

def butter_bandpass(lowcut, highcut, fs, order=5):
   nyq = 0.5 * fs
   low = lowcut / nyq
   high = highcut / nyq
   b, a = butter(order, [low, high], btype='band')
   return b, a

def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
   b, a = butter_bandpass(lowcut, highcut, fs, order=order)
   y = lfilter(b, a, data)
   return y

t = np.linspace(0,1,1000)    
s = np.sin(2*np.pi*100*t) + np.sin(2*np.pi*10*t)

T = t[1] - t[0]
fs = 1/T
print(fs)
y = butter_bandpass_filter(s, 90,110,fs,order=5)


plt.figure(1)    
plt.plot(t,s,'b',t,y,'r')    

plt.figure(2)
N = s.size
f = np.linspace(0, 1/T, N)
s -= np.mean(s)
fft = np.fft.fft(s)
plt.bar(f[:N // 2], np.abs(fft)[:N // 2] * 1 / N, width=5)

plt.figure(3)
N = y.size
f = np.linspace(0, 1/T, N)
y -= np.mean(y)
fft = np.fft.fft(y)
plt.bar(f[:N // 2], np.abs(fft)[:N // 2] * 1 / N, width=5)

plt.show()  

The dft shows the results I expected. So my first question would be whether this procedure for a single filter is reasonable?

My second question: How could I use this code to build a third octave bandpass filter with the following center frequencies (for example):

f_one_third = [10,12.5,16,20,25,31.5,40,50,63,80,100,125,160,200,250,315,400,500,630,800,1000,1250,1600,2000,2500,3150,4000,5000,6300,8000,10000,12500,16000,20000]

Simply using one filter after the other seems to be very inefficient and I guess that there is a much more efficient. The first link gives an approach for a third octave filter, but I am honestly not sure how to "use". I am basically missing the equivalent of the def butter_bandpass_filter function from the code above.

I would appreciate any help immensely!

$\endgroup$
0
$\begingroup$

The dft shows the results I expected. So my first question would be whether this procedure for a single filter is reasonable?

Sort of. For very low frequencies and/or high orders, using transfer function representation (b,a) runs into numerical problems. Use second order section representation and sosfilt() instead.

My second question: How could I use this code to build a third octave bandpass filter with the following center frequencies (for example):

Depends a bit on what you are planning to do with the resulting signals. Just doing lots of parallel filters is not an unreasonable approach. The alternative would be a frequency domain method, but this is unlikely to be cheaper if you need all results in the time domain. However, if you just want to build a 3rd octave spectrum analyzer, you just need the time averaged band energy as the function of time, which can be done in the frequency domain pretty easily.

$\endgroup$
  • $\begingroup$ Thanks! Alright, I will use sosfilt() then. My desired final ouput is a frequency spectrum in which the amplitudes are assigned to the third octave center frequencies, I think that this is a 3rd octave spectrum analyzer. Does that mean that I don't even need a bandpass filter like the one in the code above? $\endgroup$ – Larzeb Mar 5 at 19:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.