Fréquence des filtres de Bandpass butterworth en principe

je conçois un filtre bandpass en principe suivant le cookbook . Cependant, si je désapprouve trop les fréquences de filtrage, je me retrouve avec des déchets à des filtres de haut niveau. Ce que je fais mal? 

from scipy.signal import butter, lfilter

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

if __name__ == "__main__":
    import numpy as np
    import matplotlib.pyplot as plt
    from scipy.signal import freqz  
    # Sample rate and desired cutoff frequencies (in Hz).
    fs = 25
    # Plot the frequency response for a few different orders.
    plt.figure(1)
    plt.clf()
    for order in [1, 3, 5, 6, 9]:
        b, a = butter_bandpass(0.5, 4, fs, order=order)
        w, h = freqz(b, a, worN=2000)#np.logspace(-4, 3, 2000))
        plt.semilogx((fs * 0.5 / np.pi) * w, abs(h), label="order = %d" % order)  
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Gain')
    plt.grid(True)
    plt.legend(loc='best')

    plt.figure(2)
    plt.clf()
    for order in [1, 3, 5, 6, 9]:
        b, a = butter_bandpass(0.05, 0.4, fs, order=order)
        w, h = freqz(b, a, worN=2000)#np.logspace(-4, 3, 2000))
        plt.semilogx((fs * 0.5 / np.pi) * w, abs(h), label="order = %d" % order)  
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Gain')
    plt.grid(True)
    plt.legend(loc='best')

    plt.show()

fs = 25, low = 0.5, high = 4 fs = 25, low = 0.05, high = 0.4

mise à jour: la question a été examinée et apparemment résolue sur Scipy 0.14. Cependant, le tracé semble toujours très mauvais après la mise à jour de Scipy. Ce qui est de mal?

After Scipy 0.14, even worse

6
demandé sur Fra 2014-02-18 22:51:07

4 réponses

  1. N'utilisez pas b, a = butter pour les filtres haut de gamme, que ce soit dans Matlab ou SciPy ou Octave. Le format de la fonction de transfert a des problèmes de stabilité numérique , parce que certains coefficients sont très grands tandis que d'autres sont très petits. C'est pourquoi nous avons modifié les fonctions de conception du filtre pour utiliser le format ZPK en interne . Pour voir les avantages de ceci, vous devez utiliser z, p, k = butter(output='zpk') et puis travailler avec des pôles et des zéros au lieu de le numérateur et le dénominateur.
  2. ne faites pas de filtres numériques haut de gamme en une seule étape. C'est une mauvaise idée, quel que soit le Logiciel ou le matériel sur lequel vous les implémentez. Généralement, il est préférable de les diviser en sections de second ordre . Dans Matlab, vous pouvez utiliser zp2sos pour les générer automatiquement. cela n'a pas encore été implémenté en SciPy , bien que j'ai le code GPL pour le faire ici: https://gist.github.com/endolith/4525003 (licence GPL incompatible avec scipy, elle doit donc être réécrite à partir de zéro.)
3
répondu endolith 2014-05-27 03:36:21

apparemment, le problème est un bug connu:

Github

1
répondu Fra 2014-02-18 20:02:32

C'est un problème courant dans les filtres numériques. Les filtres d'ordre élevé avec des fréquences de coupure bien inférieures à la fréquence de nyquist ont tendance à avoir des coefficients instables en raison de la précision limitée des nombres à virgule flottante. La dernière fois que J'ai vérifié (il y a quelques années, il est vrai), Matlab a fait un bien meilleur travail de conservation de la précision que scipy, bien qu'il y ait encore des problèmes avec des filtres suffisamment extrêmes.

il y a quelques options si vous ne pouvez pas utiliser matlab. La première consiste à fractionner votre filtre en sections de second ordre en cascade. En gros, vous calculez vos pôles et zéros désirés, vous les divisez en paires conjuguées complexes et vous calculez la fonction de transfert pour chaque paire.

la deuxième option consiste à rééchantillonner à une fréquence d'échantillonnage plus proche de la fréquence du filtre. Par exemple, dans votre deuxième exemple, votre taux d'échantillonnage est de 25 et votre fréquence de coupure la plus élevée est .4. Vous pouvez utiliser un filtre passe-bas anti-aliasing et puis décimer par un facteur de 10 pour un taux d'échantillonnage de 2,5. Avec un taux d'échantillonnage plus faible, les coefficients de votre filtre bandpass seront moins sensibles aux erreurs d'arrondissement. Si vous faites cela, vous devez vous assurer que le filtre anti-repliement n'a pas le même problème.

0
répondu Evan 2014-02-21 18:22:52

ce qui se passe, c'est que les ordres des filtres de bandpass (BP) créés dans le script sont en fait le double de ceux affichés dans le graphique . Rappelons que l'ordre du filtre est de l'ordre du polynôme au dénominateur de la fonction de transfert. Canonic bandpass les filtres sont toujours de ordre Pair .

les nombres indiqués sont les ordres du prototype low-pass (LP) (généralement normalisé à une fréquence de coupure de 1 rad / s) qui est utilisé pour appliquer la transformation LP-TO-BP qui double l'ordre du filtre. Ainsi, par exemple, si nous commençons avec un LP d'ordre 1, nous finissons avec un second ordre bandpass:

1/(S+1) = > LP-2-BP transf. = > K. s / (S^2+a. s+b)

Où k , un et b sont des constantes. Le numérateur d'un filtre à bande standard est K. s^ (N/2) et donc l'ordre N du filtre doit être pair.

cette question d'ordre pour bandpass (qui arrive aussi aux entailles ou aux filtres de rejet de bande) n'est pas mentionnée dans SciPy documentation . En fait, si vous imprimez la longueur du dénominateur a (en utilisant print(len(a)) ) avant plt.show() , vous verrez il a 19 coefficients, ce qui correspond à un polynôme de 18ème ordre.

0
répondu user29690 2016-05-14 03:51:51