1/3 octave from audio file with python

Question

I am a beginner in signal processing and I would like to apply third-octave band pass filters on mp3 or wav file (yields about 30 new filtered time series) center frequencies: 39 Hz, 50 Hz, 63 Hz, 79 Hz, 99 Hz, 125 Hz, 157 Hz, 198 Hz, 250 Hz, 315 Hz, 397 Hz, 500 Hz, …

First way...

After I read mp3 file, I got a stereo signal. Then I segmented 1 audio file into 4096 samples. Then I separate it to left and right channel. Now I have arrays of data for each channel. Next, I apply Fast Fourier Transform to a sample file. The problem is I need to apply third-octave band pass filters to these samples. I need suggestion on how should I do since I quite not understand acoustics library.

Another way ...

I found some website related to my expectation but he uses octave bandpass filter. I use the code from Michael's reply on https://www.dsprelated.com/thread/7036/octave-bandpass-filter-on-audio-wav-files So I would like to apply this code to third-octave.

from scipy import signal
import numpy as np
import matplotlib.pyplot as plt
import math

sampleRate = 44100.0
nyquistRate = sampleRate/2.0

#center = [39, 50, 63, 79, 99, 125, 157, 198, 250, 315, 397, 500, 630, 
794, 1000, 1260, 1588, 2000, 2520, 3176, 4000, 5040, 6352, 8000, 10080, 
12704, 16000, 20160, 2508, 32000]
centerFrequency_Hz = 480.0;
lowerCutoffFrequency_Hz=centerFrequency_Hz/math.sqrt(2);
upperCutoffFrequenc_Hz=centerFrequency_Hz*math.sqrt(2);

# Determine numerator (b) and denominator (a) coefficients of the digital 
# Infinite Impulse Response (IIR) filter.
b, a = signal.butter( N=4, Wn=np.array([ lowerCutoffFrequency_Hz, 
upperCutoffFrequenc_Hz])/nyquistRate, btype='bandpass', analog=False, 
output='ba');

# Compute frequency response of the filter.
w, h = signal.freqz(b, a)

fig = plt.figure()
plt.title('Digital filter frequency response')
ax1 = fig.add_subplot(111)

plt.plot(w, 20 * np.log10(abs(h)), 'b')
plt.ylabel('Amplitude [dB]', color='b')
plt.xlabel('Frequency [rad/sample]')

ax2 = ax1.twinx()
angles = np.unwrap(np.angle(h))
plt.plot(w, angles, 'g')
plt.ylabel('Angle (radians)', color='g')
plt.grid()
plt.axis('tight')
plt.show()

fs, speech = wavfile.read(filename='segmented/atb30.wav');
speech = speech[:, 0]
fig=plt.figure()
plt.title('Speech Signal')
plt.plot(speech)

filteredSpeech=signal.filtfilt(b, a, speech)
fig=plt.figure()
plt.title('480 Hz Octave-band Filtered Speech')
plt.plot(filteredSpeech)

Answer 1

According to equations (5) and (6) from ANSI S1.11: Specification for Octave, Half-Octave, and Third Octave Band Filter Sets, for 1/3-octave the lower and upper frequencies of each band are given by:

factor = np.power(G, 1.0/6.0)
lowerCutoffFrequency_Hz=centerFrequency_Hz/factor;
upperCutoffFrequency_Hz=centerFrequency_Hz*factor;

Where G is either 2 (when designing filters according to the specified base-2 rules), or np.power(10, 0.3) (when designing filters according to the specified base-10 rules). In your case the center frequencies you have provided are close to the values obtained using base-2 rules, so you should also G = 2 to be consistent.

Note that for your given center frequency array, a few values of the upper frequencies will be greater than the Nyquist frequency (half your sampling rate). Those would correspondingly yield invalid upper normalized frequency inputs when trying to design the filter with scipy.signal.butter. To avoid that you should limit your center frequency array to values less than ~19644Hz:

centerFrequency_Hz = np.array([39, 50, 63, 79, 99, 125, 157, 198, 250, 315, 
397, 500, 630, 794, 1000, 1260, 1588, 2000, 2520, 3176, 4000, 5040, 6352, 8000, 10080, 
12704, 16000])

Also, scipy.signal.butter can handle one set of lower and upper frequencies at a time, so you should loop over the lower and upper frequency arrays to design each bandpass filter:

for lower,upper in zip(lowerCutoffFrequency_Hz, upperCutoffFrequency_Hz):
    # Determine numerator (b) and denominator (a) coefficients of the digital 
    # Infinite Impulse Response (IIR) filter.
    b, a = signal.butter( N=4, Wn=np.array([ lower, 
    upper])/nyquistRate, btype='bandpass', analog=False, 
    output='ba');

    # Compute frequency response of the filter.
    w, h = signal.freqz(b, a)

    plt.plot(w, 20 * np.log10(abs(h)), 'b')

    # Filter signal
    filteredSpeech = signal.filtfilt(b, a, speech)

This should give you a plot similar to the following for the magnitude responses:

You may at this point notice some signs of instability in the lowest band. To avoid this problem, you could switch to the sos representation:

for lower,upper in zip(lowerCutoffFrequency_Hz, upperCutoffFrequency_Hz):
    # Design filter
    sos = signal.butter( N=4, Wn=np.array([ lower, 
    upper])/nyquistRate, btype='bandpass', analog=False, 
    output='sos');

    # Compute frequency response of the filter.
    w, h = signal.sosfreqz(sos)

    plt.plot(w, 20 * np.log10(abs(h)), 'b')

    # Filter signal
    filteredSpeech = signal.sosfiltfilt(sos, speech)