I actually work on a real time graphic equalizer on python. I'm using pyaudio module, scipy, numpy. My equalizer is based on a third octave band filter bank from 25 Hz to 20 kHz (so 30 bands). This filter bank divides an input signal into 30 filtered signals (centered on the center frequency of each third octave band). Also, the streaming is implemented block by block (using pyaudio and callback method).
I used filtfilt from scipy.signal module but I had some discontinuities between each block (some audible click). So, I've followed Continuity issue when applying an IIR filter on successive time-frames and it works well for high frequencies.
But for low frequencies I need to follow these the steps :
1) downsampling input signal (to keep a good definition's filter);
2) filtering with lfilter_zi to keep continuity between each block (for streaming);
3) upsampling the filtered signal.
My problem is the upsampling because that breaks the continuity between each block (see figure below) Discontinuities between 2 blocks when downsampling and upsampling a signal (sinus at 1000Hz here)
Also, here is my code of third octave band filter bank :
# -*- coding: utf-8 -*-
"""
Created on Tue Feb 19 16:14:09 2019
@author: William
"""
from __future__ import division
import numpy as np
import scipy.signal as sc
def oct3_dsgn(fc, fs, type_filt='cheby2_bandpass', output='ba'):
"""
Calcul les coefficients B et A d'un filtre passe-bande pour chacune des
bandes de tiers d'octave entre 25 Hz et 20 kHz.
La fonction cheb2ord permet d'optimiser l'ordre et la bande passante du
filtre en fonction des fréquences de coupures et de leurs gains respectifs
(gpass, gstop) et de la fréquence d'échantillonnage.
"""
#------- Définition des fréquences inférieures et supérieures de la bande n
fc1 = fc / 2**(1/6)
fc2 = fc * 2**(1/6)
#------- Définition des fréquences centrales des bandes n-1 et n+1
fm1 = fc1 / 2**(1/6)
fm2 = fc2 * 2**(1/6)
#------- Définition des fréquences normalisées par rapport à f_nyquist
W1 = fc1/(fs/2)
W2 = fc2/(fs/2)
Wm1 = fm1/(fs/2)
Wm2 = fm2/(fs/2)
#------- Définition des filtres passe-bande, passe-bas et passe-haut
if type_filt == 'cheby2_bandpass':
gpass = 20*np.log10(np.sqrt(2)) # Équivalent à 3dB
gstop = 45
n_filt, Wn = sc.cheb2ord([W1, W2], [Wm1, Wm2], gpass, gstop)
if output=='ba':
B, A = sc.cheby2(n_filt, gstop, Wn, btype='band', output='ba')
elif output=='sos':
sos = sc.cheby2(n_filt, gstop, Wn, btype='band', output='sos')
elif type_filt == 'butter':
gpass = 20*np.log10(np.sqrt(2)) # Équivalent à 3dB
gstop = 30
n_filt, Wn = sc.buttord([W1, W2], [Wm1, Wm2], gpass, gstop)
if output=='ba':
B, A = sc.butter(n_filt, Wn, btype='band', output='ba')
elif output=='sos':
sos = sc.cheby2(n_filt, gstop, Wn, btype='band', output='sos')
elif type_filt == 'cheby2_low':
gpass = 20*np.log10(np.sqrt(2)) # Équivalent à 3dB
gstop = 45
n_filt, Wn = sc.cheb2ord(W2, Wm2, gpass, gstop)
if output=='ba':
B, A = sc.cheby2(n_filt, gstop, Wn, btype='low', output='ba')
elif output=='sos':
sos = sc.cheby2(n_filt, gstop, Wn, btype='low', output='sos')
elif type_filt == 'cheby2_high':
gpass = 20*np.log10(np.sqrt(2)) # Équivalent à 3dB
gstop = 45
n_filt, Wn = sc.cheb2ord(W1, Wm1, gpass, gstop)
if output=='ba':
B, A = sc.cheby2(n_filt, gstop, Wn, btype='high', output='ba')
elif output=='sos':
sos = sc.cheby2(n_filt, gstop, Wn, btype='high', output='sos')
if output == 'ba':
return B, A
elif output == 'sos':
return sos
def oct3_filter(signal, fs, gain_general = 0, gain_bande = np.zeros((30, )), plot=False):
"""
Calcul le signal filtré à partir du signal initial grâce une reconstruction
parfaite bande par bande. Le signal initial est filtré pour chacune des bandes.
Lorsque les fréquences sont trop basses, un sous-échantillonnage est opéré
pour gagner en résolution fréquentielle et en bande passante.
Pour la bande à 25 Hz, un filtre passe-bas est utilisé. La bande à 25 Hz
comprend donc toutes les fréquences inférieures ou égales à 25 Hz.
De même, pour la bande à 20 kHz un filtre passe-haut est utlisé. La bande à
20 kHz comprend donc toutes les fréquences au-dessus de 18 kHz.
"""
#------- Définition des fréquences centrales exactes en base 2
fc_oct3 = (1000) * (2**(1/3))**np.arange(-16, 14)
n_bande = len(fc_oct3)
n_signal = len(signal)
#------- Définition des matrices de stockage des signaux et des gains
signal_filt = np.zeros((n_signal, n_bande))
Gain = 10**(gain_bande/20)
#------- Affichage de la réponse des filtres 1/3 d'octaves
if plot == True:
import matplotlib.pyplot as plt
plt.figure(0, figsize=(10,5))
#------- Boucle sur les bandes de 10 Hz à 20 kHz
for ii in range(0, n_bande):
if ii == n_bande-1 :
## bande à 20 kHz
sos = oct3_dsgn(fc_oct3[ii], fs, type_filt='cheby2_high', output='sos')
signal_filt[:, ii] = Gain[ii] * sc.sosfilt(sos, signal)
## affichage de la réponse du filtre
if plot == True:
w, h = sc.freqz(sos)
plt.semilogx(w/2/np.pi*fs, 20*np.log10(abs(h)), 'k')
elif ii == 0 :
## bande à 25 Hz
n_decimate = 32
x = sc.decimate(signal, n_decimate)
sos = oct3_dsgn(fc_oct3[ii], fs//n_decimate, type_filt='cheby2_low', output='sos')
x = sc.sosfilt(sos, x)
signal_filt[:, ii] = Gain[ii] * sc.resample(x, n_signal)
## affichage de la réponse du filtre
if plot == True:
w, h = sc.freqz(sos)
plt.semilogx(w/2/np.pi*(fs//n_decimate), 20*np.log10(abs(h)), 'k')
elif n_bande-5 <= ii < n_bande-1:
## de 8 kHz à 16 kHz
sos = oct3_dsgn(fc_oct3[ii], fs, output='sos')
signal_filt[:, ii] = Gain[ii] * sc.sosfilt(sos, signal)
## affichage de la réponse du filtre
if plot == True:
w, h = sc.freqz(sos)
plt.semilogx(w/2/np.pi*fs, 20*np.log10(abs(h)), 'k')
elif n_bande-10 <= ii < n_bande-5:
## de 2,5 kHz à 6,3 kHz
n_decimate = 2
x = sc.decimate(signal, n_decimate)
sos = oct3_dsgn(fc_oct3[ii], fs//n_decimate, output='sos')
x = sc.sosfilt(sos, x)
signal_filt[:, ii] = Gain[ii] * sc.resample(x, n_signal)
## affichage de la réponse du filtre
if plot == True:
w, h = sc.freqz(sos)
plt.semilogx(w/2/np.pi*(fs//n_decimate), 20*np.log10(abs(h)), 'k')
elif n_bande-15 <= ii < n_bande-10:
## de 800 Hz à 2 kHz
n_decimate = 4#round(fs/(2*fmax))-1
x = sc.decimate(signal, n_decimate)
sos = oct3_dsgn(fc_oct3[ii], fs//n_decimate, output='sos')
x = sc.sosfilt(sos, x)
signal_filt[:, ii] = Gain[ii] * sc.resample(x, n_signal)
## affichage de la réponse du filtre
if plot == True:
w, h = sc.freqz(sos)
plt.semilogx(w/2/np.pi*(fs//n_decimate), 20*np.log10(abs(h)), 'k')
elif n_bande-20 <= ii < n_bande-15:
## de 250 Hz à 630 Hz
n_decimate = 8#round(fs/(2*fmax))-1
x = sc.decimate(signal, n_decimate)
sos = oct3_dsgn(fc_oct3[ii], fs//n_decimate, output='sos')
x = sc.sosfilt(sos, x)
signal_filt[:, ii] = Gain[ii] * sc.resample(x, n_signal)
### affichage de la réponse du filtre
if plot == True:
w, h = sc.freqz(sos)
plt.semilogx(w/2/np.pi*(fs//n_decimate), 20*np.log10(abs(h)), 'k')
elif n_bande-25 <= ii < n_bande-20:
## de 80 Hz à 200 Hz
n_decimate = 16#round(fs/(2*fmax))-1
x = sc.decimate(signal, n_decimate)
sos = oct3_dsgn(fc_oct3[ii], fs//n_decimate, output='sos')
x = sc.sosfilt(sos, x)
signal_filt[:, ii] = Gain[ii] * sc.resample(x, n_signal)
## affichage de la réponse du filtre
if plot == True:
w, h = sc.freqz(sos)
plt.semilogx(w/2/np.pi*(fs//n_decimate), 20*np.log10(abs(h)), 'k')
elif n_bande-29 <= ii < n_bande-25:
## de 25 Hz à 63 Hz
n_decimate = 32#round(fs/(2*fmax))-1
x = sc.decimate(signal, n_decimate)
sos = oct3_dsgn(fc_oct3[ii], fs//n_decimate, output='sos')
x = sc.sosfilt(sos, x)
signal_filt[:, ii] = Gain[ii] * sc.resample(x, n_signal)
## affichage de la réponse du filtre
if plot == True:
w, h = sc.freqz(sos)
plt.semilogx(w/2/np.pi*(fs//n_decimate), 20*np.log10(abs(h)), 'k')
if plot == True:
plt.grid(which='both', linestyle='-', color='grey')
# plt.xticks([20, 50, 100, 200, 500, 1000, 2000, 5000, 10000, 20000],
# ["20", "50", "100", "200", "500", "1K",
# "2K", "5K", "10K", "20K"])
plt.xlabel('Fréquence [Hz]'), plt.ylabel('Gain [dB]')
plt.title('Réponse en fréquence des filtres 1/3 d\'octaves')
plt.xlim((10, 22e3)), plt.ylim((-5, 1))
plt.show()
#------- Sommation des signaux filtrés pour recomposer le signal d'origine
S = signal_filt.sum(axis=1)
S = S - np.mean(S)
## tuckey_window = sc.tukey(len(S), alpha=0.01)
## S = tuckey_window * S
G = 10**(gain_general/20)
return G * S
