对具有初始条件的连续块进行滤波后重新采样（以避免不连续）

Question

我实际上在python上开发一个实时图形均衡器。我正在使用pyaudio模块，scipy，numpy。我的均衡器是基于第三倍频程滤波器组从25赫兹到20千赫（所以30波段）。该滤波器组将输入信号分成30个滤波信号（以每三个倍频程的中心频率为中心）。此外，流是逐块实现的（使用pyaudio和回调方法）。你知道吗

我使用了来自scipy.signal模块的filtfilt，但是我在每个块之间有一些不连续性（一些可听见的咔嗒声）。所以，我遵循了Continuity issue when applying an IIR filter on successive time-frames，它适用于高频。你知道吗

但对于低频，我需要遵循以下步骤：

1）对输入信号进行下采样（以保持良好的滤波器清晰度）

2）使用lfilter_zi进行过滤以保持每个块之间的连续性（对于流）

3）对滤波后的信号进行上采样。你知道吗

我的问题是上采样，因为这会破坏每个块之间的连续性（参见下图） Discontinuities between 2 blocks when downsampling and upsampling a signal (sinus at 1000Hz here)

另外，这是我的第三倍频程滤波器组代码：

    # -*- 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