def _init_sub(self):
        """further init function"""
        sys = signal.lti(self._num, self._den)
        
        if self._fmin is not None and self._fmax is not None :
            wnp = np.logspace(np.log10(self._fmin*2*np.pi),
                              np.log10(self._fmax*2*np.pi),
                              num=1000) # consider frequencies
            w, mag, ph = signal.bode(sys, w=wnp)
            f = w/(2.*np.pi) # consider frequencies
        else:
            w, mag, ph = signal.bode(sys, n=1000)
            f = w/(2.*np.pi) # consider frequencies
            self._fmin = f[0]
            self._fmax = f[-1]

        #self.aM.tick_params(labelsize='x-small')
        #self.aP.tick_params(labelsize='x-small')
        
        self._aM_bl, = self.aM.semilogx(f, mag)
        self.aM.set_xlim([self._fmin, self._fmax])
        [ymin, ymax] = self.aM.get_ylim()
        dy = (ymax - ymin) * 0.05
        self.aM.set_ylim([ymin - dy, ymax + dy])
        self.aM.set_ylabel('magnitude (dB)', size='small')
        
        self._aP_bl, = self.aP.semilogx(f, ph)
        self.aP.set_xlim([self._fmin, self._fmax])
        [ymin, ymax] = self.aP.get_ylim()
        dy = (ymax - ymin) * 0.05
        self.aP.set_ylim([ymin - dy, ymax + dy])
        self.aP.set_ylabel('phase (degrees)', size='small')
        
        self._aM_pl, = self.aM.plot([], [], 'ro')
        self._aP_pl, = self.aP.plot([], [], 'ro')

        viniz = (np.log10(self._fmax) - np.log10(self._fmin))/4. + \
              np.log10(self._fmin)
        self._slider = Slider(self.aS, r'freq', np.log10(self._fmin),
                              np.log10(self._fmax), valinit=viniz)
        self._slider.vline.set_alpha(0.)
        self._slider.label.set_fontsize('small')
        self._slider.on_changed(self._onchange)
        txts = filter(lambda x: type(x)==mpl.text.Text, self.aS.get_children())
        txts[1].set_visible(False)

        self.aT.set_axis_bgcolor('0.95')
        self.aT.xaxis.set_ticks_position('none')
        self.aT.yaxis.set_ticks_position('none')
        self.aT.spines['right'].set_visible(False)
        self.aT.spines['left'].set_visible(False)
        self.aT.spines['bottom'].set_visible(False)
        self.aT.spines['top'].set_visible(False)
        self.aT.set_xticks([])
        self.aT.set_yticks([])
        self._txt['f'] = self.aT.text(0, 0, '')
        self._txt['m'] = self.aT.text(0.34, 0, '')
        self._txt['p'] = self.aT.text(0.73, 0, '')

def second_ordre(f0, Q, filename='defaut.png', type='PBs', f=None):
    '''Petite fonction pour faciliter l'utilisation de la fonction "bode" du 
    module "signal" quand on s'intéresse à des filtres du 2e ordre. Il suffit 
    de donner la fréquence propre f0 et le facteur de qualité Q pour obtenir 
    ce que l'on veut. Autres paramètres:
    * filename: le nom du fichier ('defaut.png' par défaut)
    * type: le type du filtre, à choisir parmi 'PBs' (passe-bas), 
    'PBd' (passe-bande) et 'PHt' (passe-haut). On peut aussi définir soi-même 
    le numérateur sous forme d'une liste de plusieurs éléments, le degré le 
    plus haut donné en premier. NB: le '1.01' des définitions est juste là 
    pour améliorer sans effort le rendu graphique.
    * f: les fréquences à échantillonner (si None, la fonction choisit 
    d'elle-même un intervalle adéquat).
    '''
    den = [1. / f0**2, 1. /
           (Q * f0), 1]              # Le dénominateur de la fonction de transfert
    if type == 'PBs':
        num = [1.01]          # Le numérateur pour un passe-bas
    elif type == 'PBd':
        num = [1.01 / (Q * f0), 0]  # pour un passe-bande
    elif type == 'PHt':
        num = [1.01 / f0**2, 0, 0]  # pour un passe-haut
    else:
        # sinon, c'est l'utilisateur qui le définit.
        num = type
    # Définition de la fonction de transfert
    s1 = signal.lti(num, den)
    # Obtention des valeurs adéquates
    f, GdB, phase = signal.bode(s1, f)
    # Dessin du diagramme proprement dit
    diag_bode(f, GdB, phase, filename)

def tfc (num, den, freq, mf='std', pf='rad'):
    """Calculate magnitude and phase of a transfer function for a given
    frequency.

    Parameters:
        num, den: 1D-arrays
            see scipy.signal.lti for documentation
        freq: scalar
            frequency (Hz)
        mf: string, optional, default: 'std'
            accepted values: 'std', 'dB'
        pf: string, optional, default: 'rad'
            accepted values: 'rad', 'deg'
            
    Returns: mag, ph
        mag: scalar if mf valid else None
            its unity measure depends on mf
            NB: mag_dB = 20 * log10(mag_std)
        ph: scalar if pf valid else None
            its unity measure depends on pf
    """
    mag, ph = None, None
    sys = signal.lti(num, den)
    w, magl, phl = signal.bode(sys, w=[2*np.pi*freq])
    mag_dB, ph_deg = magl[0], phl[0]
    mag_std, ph_rad = 10**(mag_dB/20.), ph_deg/180.*np.pi
    if mf == 'std':
        mag = mag_std
    elif mf == 'dB':
        mag = mag_dB
    if pf == 'rad':
        ph = ph_rad
    elif pf == 'deg':
        ph = ph_deg
    return mag, ph

def bode_plot():
            # パラメータ設定
    m = 1
    c = 1
    k = 400

    A = np.array([[0, 1], [-k/m, -c/m]])
    B = np.array([[0], [k/m]])
    C = np.array([1, 0])
    D = np.array([0])

    s1 = signal.lti(A, B,   C, D)
    w, mag, phase = signal.bode(s1, np.arange(1, 500, 1)) 

    # プロット
    plt.figure(1)
    plt.subplot(211)
    plt.loglog(w, 10**(mag/20))
    plt.ylabel("Amplitude")
    plt.axis("tight")
    plt.subplot(212)
    plt.semilogx(w, phase)
    plt.xlabel("Frequency[Hz]")
    plt.ylabel("Phase[deg]")
    plt.axis("tight")
    plt.ylim(-180, 180)
    plt.savefig('../files/150613ABCD01.svg')

    def plot_bode(self,b,c):
	s1 = signal.lti(b,c)
	w, mag, phase = signal.bode(s1)
	plt.figure()
	plt.grid()
	plt.semilog(w, mag)
	plt.semilog(w, phase)
	plt.show()

 def plot(self):
     '''
     Here goes the matplotlib magic: this function generates the Bode diagram
     '''
     w, mag, phase = signal.bode(self.sig)
     self.graphMg.clear()
     self.graphPh.clear()
     self.graphMg.plot(w, mag)
     self.graphPh.plot(w, phase, 'r')

 def test_05(self):
     # Test that bode() finds a reasonable frequency range.
     # 1st order low-pass filter: H(s) = 1 / (s + 1)
     system = lti([1], [1, 1])
     n = 10
     # Expected range is from 0.01 to 10.
     expected_w = np.logspace(-2, 1, n)
     w, mag, phase = bode(system, n=n)
     assert_almost_equal(w, expected_w)

 def test_04(self):
     # Test bode() phase calculation.
     # 1st order low-pass filter: H(s) = 1 / (s + 1)
     system = lti([1], [1, 1])
     w = [0.1, 1, 10, 100]
     w, mag, phase = bode(system, w=w)
     jw = w * 1j
     y = np.polyval(system.num, jw) / np.polyval(system.den, jw)
     expected_phase = np.arctan2(y.imag, y.real) * 180.0 / np.pi
     assert_almost_equal(phase, expected_phase)

 def test_03(self):
     # Test bode() magnitude calculation.
     # 1st order low-pass filter: H(s) = 1 / (s + 1)
     system = lti([1], [1, 1])
     w = [0.1, 1, 10, 100]
     w, mag, phase = bode(system, w=w)
     jw = w * 1j
     y = np.polyval(system.num, jw) / np.polyval(system.den, jw)
     expected_mag = 20.0 * np.log10(abs(y))
     assert_almost_equal(mag, expected_mag)

 def test_02(self):
     # Test bode() phase calculation (manual sanity check).
     # 1st order low-pass filter: H(s) = 1 / (s + 1),
     #   angle(H(s=0.1)) ~= -5.7 deg
     #   angle(H(s=1)) ~= -45 deg
     #   angle(H(s=10)) ~= -84.3 deg
     system = lti([1], [1, 1])
     w = [0.1, 1, 10]
     w, mag, phase = bode(system, w=w)
     expected_phase = [-5.7, -45, -84.3]
     assert_almost_equal(phase, expected_phase, decimal=1)

 def test_01(self):
     # Test bode() magnitude calculation (manual sanity check).
     # 1st order low-pass filter: H(s) = 1 / (s + 1),
     # cutoff: 1 rad/s, slope: -20 dB/decade
     #   H(s=0.1) ~= 0 dB
     #   H(s=1) ~= -3 dB
     #   H(s=10) ~= -20 dB
     #   H(s=100) ~= -40 dB
     system = lti([1], [1, 1])
     w = [0.1, 1, 10, 100]
     w, mag, phase = bode(system, w=w)
     expected_mag = [0, -3, -20, -40]
     assert_almost_equal(mag, expected_mag, decimal=1)

    def update(**kwargs):
        lti.update(**kwargs)
        g_impulse.data_source.data['x'], g_impulse.data_source.data['y'] = (
            impulse(lti.sys, T=t))

        g_step.data_source.data['x'], g_step.data_source.data['y'] = (
            step(lti.sys, T=t))

        g_poles.data_source.data['x'] = lti.sys.poles.real
        g_poles.data_source.data['y'] = lti.sys.poles.imag

        w, mag, phase, = bode(lti.sys)
        g_bode_mag.data_source.data['x'] = w
        g_bode_mag.data_source.data['y'] = mag
        g_bode_phase.data_source.data['x'] = w
        g_bode_phase.data_source.data['y'] = phase

        push_notebook()

    def test_from_state_space(self):
        # Ensure that bode works with a system that was created from the
        # state space representation matrices A, B, C, D.  In this case,
        # system.num will be a 2-D array with shape (1, n+1), where (n,n)
        # is the shape of A.
        # A Butterworth lowpass filter is used, so we know the exact
        # frequency response.
        a = np.array([1.0, 2.0, 2.0, 1.0])
        A = linalg.companion(a).T
        B = np.array([[0.0], [0.0], [1.0]])
        C = np.array([[1.0, 0.0, 0.0]])
        D = np.array([[0.0]])
        with suppress_warnings() as sup:
            sup.filter(BadCoefficients)
            system = lti(A, B, C, D)
            w, mag, phase = bode(system, n=100)

        expected_magnitude = 20 * np.log10(np.sqrt(1.0 / (1.0 + w**6)))
        assert_almost_equal(mag, expected_magnitude)

 def test_07(self):
     # bode() should not fail on a system with pure imaginary poles.
     # The test passes if bode doesn't raise an exception.
     system = lti([1], [1, 0, 100])
     w, mag, phase = bode(system, n=2)

 def test_06(self):
     # Test that bode() doesn't fail on a system with a pole at 0.
     # integrator, pole at zero: H(s) = 1 / s
     system = lti([1], [1, 0])
     w, mag, phase = bode(system, n=2)
     assert_equal(w[0], 0.01)  # a fail would give not-a-number

from scipy import signal
import matplotlib.pyplot as plt

s1 = signal.lti([1], [1, 1])
w, mag, phase = signal.bode(s1)

plt.figure()
plt.semilogx(w, mag)    # bode magnitude plot
plt.figure()
plt.semilogx(w, phase)  # bode phase plot
plt.show()

Pour un exercice de reconnaissance d'un filtre qui puisse jouer le rôle 
d'intégrateur dans un certain intervalle de fréquences. On en génère un du 
premier ordre (passe-bas) et un du second ordre (passe-bande).
'''

from scipy import signal        # Pour les fonctions 'lti' et 'bode'
import numpy as np              # Pour np.logspace

# Ceci est un filtre passe-bas donc potentiellement intégrateur à HF

num = [1]                       # Polynôme au numérateur   (->        1       )
den = [0.01,0.999]              # Polynôme au dénominateur (-> 0.01*jw + 0.999)

f = np.logspace(-1,4,num=200)   # L'intervalle de fréquences considéré (échelle log)
s1 = signal.lti(num,den)        # La fonction de transfert
f,GdB,phase = signal.bode(s1,f) # Fabrication automatique des données

from bode import diag_bode      # Pour générer le diagramme de Bode

# Appel effectif à la fonction dédiée.
diag_bode(f,GdB,phase,'PNG/S11_integrateur.png') 

# Ceci est un filtre passe-bande du second ordre (intégrateur à HF)

num2 = [0.01,0]                 # Numérateur   (->           0.01*jw         )
den2 = [10**-2,0.01,1]          # Dénominateur (-> 0.01*(jw)**2 + 0.01*jw + 1)

f = np.logspace(-1,4,num=200)   # Intervalle de fréquences en échelle log 
s2 = signal.lti(num2,den2)      # Fonction de transfert
f,GdB,phase = signal.bode(s2,f) # Fabrication des données
diag_bode(f,GdB,phase,'PNG/S11_integrateur2.png') # et du diagramme

num = [1]
# Coefficients in denominator of transfer function
# High order to low order, eg 1*s^2 + 0.1*s + 1
den = [1, 0.1, 1]

# Scan over zeta, a parameter for a second-order system
zetaRange = np.arange(0.1,1.1,0.1)

f1 = plt.figure()
print "no prob"
# for i in range(0,2):
den = [1, 2*zetaRange[1], 1]
print den
s1 = signal.lti(num, den)
# Specify our own frequency range: np.arange(0.1, 5, 0.01)
w, mag, phase = signal.bode(s1, np.arange(0.1, 5, 0.01).tolist())
plt.semilogx (w, mag, color="blue", linewidth="1")

plt.xlabel ("Frequency")
plt.ylabel ("Magnitude")
# plt.savefig ("c:\\mag.png", dpi=300, format="png")

plt.figure()

for i in range(0,9):
    den = [1, zetaRange[i], 1]
    s1 = signal.lti(num, den)
    w, mag, phase = signal.bode(s1, np.arange(0.1, 10, 0.02).tolist())
    plt.semilogx (w, phase, color="red", linewidth="1.1")
plt.xlabel ("Frequency")
plt.ylabel ("Amplitude")

"""Tests for `gwpy.plot.filter`
"""

import numpy
from numpy import testing as nptest

from scipy import signal

from ...frequencyseries import FrequencySeries
from .. import BodePlot
from .utils import FigureTestBase

# design ZPK for BodePlot test
ZPK = [100], [1], 1e-2
FREQUENCIES, MAGNITUDE, PHASE = signal.bode(ZPK, n=100)


class TestBodePlot(FigureTestBase):
    FIGURE_CLASS = BodePlot

    def test_init(self, fig):
        assert len(fig.axes) == 2
        maxes, paxes = fig.axes
        # test magnigtude axes
        assert maxes.get_xscale() == 'log'
        assert maxes.get_xlabel() == ''
        assert maxes.get_yscale() == 'linear'
        assert maxes.get_ylabel() == 'Magnitude [dB]'
        # test phase axes
        assert paxes.get_xscale() == 'log'

        indx = np.searchsorted(self.x, [x])[0]
        x = self.x[indx]
        y = self.y[indx]
        # update the line positions
        self.lx.set_ydata(y)
        self.ly.set_xdata(x)

        self.txt.set_text('x=%1.2f, y=%1.2f' % (x, y))
        # print('x=%1.2f, y=%1.2f' % (x, y))
        plot.draw()

# define start/end freq and freq step for bode plots
startFreq = omega_3/10
endFreq = omega_4*10
freqStep = np.ceil((endFreq-startFreq)/1000)

# bode plot for magnitude and phase
tf = signal.lti(tf_num, tf_den)
w, mag, phase = signal.bode(tf, np.arange(startFreq, endFreq, freqStep).tolist())
fig, ax = plot.subplots()

plot.semilogx(w/(2*np.pi), mag, color="blue", linewidth="2")
#plot.autoscale(enable=False, axis='both', tight=False)
plot.xlabel("Frequency (Hz)")
plot.ylabel("Magnitude (dB)")
plot.axis([omega_3/(2*np.pi)/10, omega_4/(2*np.pi)*10, -100, 10])
cursor = SnaptoCursor(ax, w/(2*np.pi), mag)
plot.connect('motion_notify_event', cursor.mouse_move)

plot.show()