def butterworth_plot(fig=None, ax=None):
    """
    Plot of frequency response of the Butterworth filter with different orders.
    """

    if fig is None:
        fig, ax = plt.subplots(1, 2, figsize=(10, 4))
        
    b1, a1 = signal.butter(1, 10, 'low', analog=True)
    w, h1 = signal.freqs(b1, a1)
    ang1 = np.rad2deg(np.unwrap(np.angle(h1)))
    h1 = 20 * np.log10(abs(h1))
    b2, a2 = signal.butter(2, 10, 'low', analog=True)
    w, h2 = signal.freqs(b2, a2)
    ang2 = np.rad2deg(np.unwrap(np.angle(h2)))
    h2 = 20 * np.log10(abs(h2))
    b4, a4 = signal.butter(4, 10, 'low', analog=True)
    w, h4 = signal.freqs(b4, a4)
    ang4 = np.rad2deg(np.unwrap(np.angle(h4)))
    h4 = 20 * np.log10(abs(h4))
    b6, a6 = signal.butter(6, 10, 'low', analog=True)
    w, h6 = signal.freqs(b6, a6)
    ang6 = np.rad2deg(np.unwrap(np.angle(h6)))
    h6 = 20 * np.log10(abs(h6))
    w = w/10

    # PLOT
    ax[0].plot(w, h1, 'b', w, h2, 'r', w, h4, 'g', w, h6, 'y', linewidth=2)
    ax[0].axvline(1, color='black') # cutoff frequency
    ax[0].scatter(1, -3, marker='s', edgecolor='0', facecolor='1', s=400)
    #ax1.legend(('1', '2', '4', '6'), title='Filter order', loc='best')
    ax[0].set_xscale('log')
    fig.suptitle('Bode plot for low-pass Butterworth filter with different orders',
                 fontsize=16, y=1.05)
    #ax1.set_title('Magnitude', fontsize=14)
    ax[0].set_xlabel('Frequency / Critical frequency', fontsize=14)
    ax[0].set_ylabel('Magnitude [dB]', fontsize=14)
    ax[0].set_xlim(0.1, 10)
    ax[0].set_ylim(-120, 10)
    ax[0].grid(which='both', axis='both')
    ax[1].plot(w, ang1, 'b', w, ang2, 'r', w, ang4, 'g', w, ang6, 'y', linewidth=2)
    ax[1].axvline(1, color='black')  # cutoff frequency
    ax[1].legend(('1', '2', '4', '6'), title='Filter order', loc='best')
    ax[1].set_xscale('log')
    #ax2.set_title('Phase', fontsize=14)
    ax[1].set_xlabel('Frequency / Critical frequency', fontsize=14)
    ax[1].set_ylabel('Phase [degrees]', fontsize=14)
    ax[1].set_yticks(np.arange(0, -300, -45))
    ax[1].set_ylim(-300, 10)
    ax[1].grid(which='both', axis='both')
    plt.tight_layout(w_pad=1)
    axi = plt.axes([.115, .4, .15, .35])  # inset plot
    axi.plot(w, h1, 'b', w, h2, 'r', w, h4, 'g', w, h6, 'y', linewidth=2)
    #ax11.set_yticks(np.arange(0, -7, -3))
    axi.set_xticks((0.6, 1, 1.4))
    axi.set_yticks((-6, -3, 0))
    axi.set_ylim([-7, 1])
    axi.set_xlim([.5, 1.5])
    axi.grid(which='both', axis='both')

    def __init__(self, fp, fs, gpass, gstop, ftype, btype): 
        #Variables init.
        self.fp = fp
        self.fs = fs
        self.gpass =  gpass
        self.gstop = gstop
        self.ftype = ftype
        self.btype = btype

        #Filter type for plot's title.
        types_dict = {"butter":"Butterworth", "cheby1":"Chebyshev I", "cheby2":"Chebyshev II", "ellip": "Cauer"}
        self.ftype_plot = types_dict[ftype]

        self.__wsk()
        self.__filter_order()
            
        #Designing filter.
        (self.b, self.a) = signal.iirfilter(self.ord,
                                            self.wn,
                                            rp=self.gpass, 
                                            rs=self.gstop, 
                                            btype=self.btype, 
                                            analog=True, 
                                            output='ba', 
                                            ftype=ftype)

        #Frequency response of analog filter.
        (self.w, self.h) = signal.freqs(self.b, self.a, worN=1000)

        #Denormalizing variabels for ploting. Pulsation to frequency.
        self.w = (self.w * (self.sampling_w / 2)) / (2 * pi)
        self.wn = (self.wn * (self.sampling_w / 2)) / (2 * pi)

    def filterba(self, b, a, inplace=False):
        """Apply a filter to this `Spectrum` in numerator-denominator
        format.

        Parameters
        ----------
        b : :class:`~numpy.ndarray`
            Numerator of a linear filter
        a : :class:`~numpy.ndarray`
            Decnominator of a linear filter
        inplace : `bool`, optional, default: `False`
            modify this `Spectrum` in-place

        Returns
        -------
        Spectrum
            either a view of the current `Spectrum` with filtered data,
            or a new `Spectrum` with the filtered data
        """
        fresp = abs(signal.freqs(b, a, self.frequencies)[1])
        if inplace:
            self *= fresp
            return self
        else:
            new = self * fresp
            return new

  	def FilterResponse(self):
  		"""
  		Gives the response of the filter y frequency-amplitude
  		"""

  		self.w, self.h = SGN.freqs(self.ba, self.aa)
  		return self.w, self.h

 def compute_frequencies(self, N=None):
     if hasattr(self, 'sample_rate'):
         try:
             self.W, self.H = signal.freqz(self.B, self.A)
         except:
             self.W, self.H = signal.freqz(self.B)
     else:
         self.W, self.H = signal.freqs(self.B, self.A, N)

def plot_H(b, a):
    ws = np.linspace(2*np.pi*3.7, 2*np.pi*8.6, 200)
    w, h = signal.freqs(b, a, ws)
    plt.semilogx(w/(2*np.pi), 20 * np.log10(abs(h)))
    plt.title('Butterworth filter frequency response')
    plt.xlabel('Frequency [radians / second]')
    plt.ylabel('Amplitude [dB]')
    plt.show()

def make_ctle(rx_bw, peak_freq, peak_mag, w, dc_offset=0):
    """
    Generate the frequency response of a continuous time linear
    equalizer (CTLE), given the:

    - signal path bandwidth,
    - peaking specification, and
    - list of frequencies of interest.

    We use the 'invres()' function from scipy.signal, as it suggests
    itself as a natural approach, given our chosen use model of having
    the user provide the peaking frequency and degree of peaking.

    That is, we define our desired frequency response using one zero
    and two poles, where:

    - The pole locations are equal to:
       - the signal path natural bandwidth, and
       - the user specified peaking frequency.

    - The zero location is chosen, so as to provide the desired degree
      of peaking.

    Inputs:

      - rx_bw        The natural (or, unequalized) signal path bandwidth (Hz).

      - peak_freq    The location of the desired peak in the frequency
                     response (Hz).

      - peak_mag     The desired relative magnitude of the peak (dB). (mag(H(0)) = 1)

      - w            The list of frequencies of interest (rads./s).

      - dc_offset    The d.c. offset of the CTLE gain curve (dB).

    Outputs:

      - w, H         The resultant complex frequency response, at the
                     given frequencies.

    """

    p2   = -2. * pi * rx_bw
    p1   = -2. * pi * peak_freq
    z    = p1 / pow(10., peak_mag / 20.)
    if(p2 != p1):
        r1   = (z - p1) / (p2 - p1)
        r2   = 1 - r1
    else:
        r1   = -1.
        r2   = z - p1
    b, a = invres([r1, r2], [p1, p2], [])

    w, H = freqs(b, a, w)
    H   *= pow(10., dc_offset / 20.) / abs(H[0])  # Enforce d.c. offset.

    return (w, H)

def grp_delay_ana(b, a, w):
    """
    Calculate the group delay of an anlog filter.
    """
    w, H = sig.freqs(b, a, w)
    H_angle = np.unwrap(np.angle(H))
#    tau_g = np.zeros(len(w)-1)
    tau_g = (H_angle[1:]-H_angle[:-1])/(w[0]-w[1])
    return tau_g, w[:-1]

    def test_filter(self, array):
        a2 = array.filter([100], [1], 1e-2)
        assert isinstance(a2, type(array))
        utils.assert_quantity_equal(a2.frequencies, array.frequencies)

        # manually rebuild the filter to test it works
        b, a, = signal.zpk2tf([100], [1], 1e-2)
        fresp = abs(signal.freqs(b, a, array.frequencies.value)[1])
        utils.assert_array_equal(a2.value, fresp * array.value)

 def band_stop(self, low_f, high_f, axis='x', order=6):
     f_nyq = self.rate / 2
     band_low = low_f / f_nyq
     band_high = high_f / f_nyq
     b, a = sig.butter(order, [band_low, band_high], btype='bandstop')
     w, h = sig.freqs(b, a)
     vdata = self.ts[axis]
     filtered = sig.lfilter(b,a,vdata)
     return filtered

def butter_filter():

        N = 8 #order
        wn = 0.5 #frequency in the transition band when gain drops below -3dB
        type = 'low' #filter type

        b, a = signal.butter(N,wn,'low',analog=False)
        w, h = signal.freqs(b, a)
        return (b,a,w,h)

 def test_filter(self):
     array = self.create()
     a2 = array.filter([100], [1], 1e-2)
     self.assertIsInstance(a2, type(array))
     self.assertArraysEqual(a2.frequencies, array.frequencies)
     # manually rebuild the filter to test it works
     b, a, = signal.zpk2tf([100], [1], 1e-2)
     fresp = abs(signal.freqs(b, a, array.frequencies.value)[1])
     nptest.assert_array_equal(a2.value, fresp * array.value)

def test_sos_phys2filter():
    b, a = sos_phys2filter(S0, delta, f0)
    H = freqs(b, a, 2 * np.pi * fe5)[1]
    indmax = np.abs(H).argmax()
    assert np.round(np.abs(f0 - fe5[indmax])) <= 0.01 * f0
    assert np.abs(H[0]) == approx(S0)
    bmulti, amulti = sos_phys2filter(S0 * np.ones(K), delta * np.ones(K), f0 * np.ones(K))
    assert len(bmulti[0]) == K
    assert amulti.shape == (K, 3)

def crossCalib(monitor_trace, response_trace, **kwargs):
	
	m_trace=monitor_trace.copy()
	r_trace=response_trace.copy()

	if 'demean' in kwargs and kwargs['demean']:
		m_trace.detrend('demean')
		r_trace.detrend('demean')

	if 'taper' in kwargs and kwargs['taper']:
		m_trace.taper(0.05)
		r_trace.taper(0.05)

	#Paramètres des PSD
	if 'nfft' in kwargs:
		n_fft=kwargs['nfft']
	else:
		n_fft=1024

	if 'npad' in kwargs:
		n_pad=kwargs['npad']
	else:
		n_pad=n_fft*4

	if 'noverlap' in kwargs:
		n_overlap=kwargs['noverlap']
	else:
		n_overlap=int(n_fft*0.90)

	#paz par défaut: chaine générique
	if 'paz' in kwargs:
		paz=kwargs['paz']
	else:
		paz=dict()
		paz['zeros']=np.array([])
		paz['poles']=np.array([])
		paz['gain']=1
		paz['seismometer_gain']=1
		paz['datalogger_gain']=1
		paz['sensitivity']=paz['seismometer_gain']*paz['datalogger_gain']*paz['gain']

	
	fs=m_trace.stats.sampling_rate
	(P00,f)=mlab.psd(m_trace.data,Fs=fs,NFFT=n_fft,noverlap=n_overlap,pad_to=n_pad,detrend=mlab.detrend_mean,window=mlab.window_hanning)
	(P01,f)=mlab.csd(m_trace.data,r_trace.data,Fs=fs,NFFT=n_fft,noverlap=n_overlap,pad_to=n_pad,detrend=mlab.detrend_mean,window=mlab.window_hanning)
	(C,f)=mlab.cohere(m_trace.data,r_trace.data,Fs=fs,NFFT=n_fft,noverlap=n_overlap,pad_to=n_pad,detrend=mlab.detrend_mean,window=mlab.window_hanning)
	
	(b,a)=sp.zpk2tf(paz['zeros'],paz['poles'],paz['sensitivity'])
	(_w,H0)=sp.freqs(b,a,f*2*np.pi)

	H1=P01*H0/P00
	H1=H1[1:]
	C=C[1:]
	f=f[1:]

	return (H1,C,f)

def plot_filter_response(b, a):
    w, h = signal.freqs(b, a)
    plt.plot(w, 20 * np.log10(abs(h)))
    plt.xscale('log')
    plt.title('Chebyl1 filter frequency response')
    plt.xlabel('Frequency [radians / second]')
    plt.ylabel('Amplitude [dB]')
    plt.margins(0, 0.1)
    plt.grid(which='both', axis='both')
    plt.axvline(100, color='green')   # cutoff frequency
    plt.show()

def FDLS_fromfilt(B, A, N, D, shift = 0.0, phasemult = 1.0, analog = False):

    if analog:
        w, h = sps.freqs(B, A)
    else:
        w, h = sps.freqz(B, A, worN = 1024)

    Am = np.absolute(h)
    Th = np.angle(h) * phasemult 

    return FDLS(N, D, w, Am = Am, Th = Th, shift = shift) 

def Rule_7(w_start, w_end):
    # Rule 7 determining the phase of GGm at -180deg
    # this is solved visually from a plot

    w = np.logspace(w_start, w_end, 1000)

    Pz = np.polymul(G()[0], Gm()[0])
    Pp = np.polymul(G()[1], Gm()[1])
    [w, h] = scs.freqs(Pz, Pp, w)

    plt.semilogx(w, (180 / np.pi) * (phase(h) + w * Time_Delay()))
    plt.show()

def elliptical_filter(): #Design an elliptical filter with some predefined parameters

	N = 8 #order
	rp = 1 #max ripple below unity in passband [dB]
	rs = 10 #max attenuation in the stop band [dB]

	wn = 0.5 #frequency in the transition band when gain drops below rp (for ellip)
	type = 'low' #filter type

	b, a = signal.ellip(N,rp,rs,wn,'low',analog=False)
	w, h = signal.freqs(b, a)
	return (b,a,w,h)

 def __method(self):
     # TODO Check CPU time
     '''
     low pass filtering for mode estimation function 
     array must be declared and passed through the function 
     '''
     # step 1 vedran, design filter
     ''' Wn= length-2 sequence giving the critical frequencies - Fp, Fst parameters from matlab fdesign.lowpass(Fp,Fst,Ap,Ast)
     rp: maximum ripple in the pass band. (dB) - Ap parameter from matlab fdesign.lowpass(Fp,Fst,Ap,Ast)
     rs:  minimum attenuation in the stop band. (dB) - Ast parameter from matlab fdesign.lowpass(Fp,Fst,Ap,Ast)
     '''
     b, a= signal.iirfilter(self.__order, Wn=[2/25,2.5/25], rp=0.1, rs=50, btype='lowpass', 
                            analog=True, ftype='cheby2')
     # step 2 vedran, apply filter
     simsignalFilter= signal.lfilter(b, a, self.__simulationsSignal)
     # step 3 vedran, downsample the signal
     simsignalsampled = signal.decimate(simsignalFilter, 10, ftype='iir')
     # step 4 with signal.freqs(b,a,y) we obtain frequency response of the signal
     freqHz, amplitude = signal.freqs(b, a, simsignalsampled)
     print 'angular frequency ', freqHz
     print 'amplitude response ', amplitude
     for vfreq, vdamp in zip(freqHz, amplitude):
         mode= EigenValue(vfreq,vdamp)
         self.__simulationModes.append(mode)
     # step 2 vedran, apply filter
     maessignalFilter= signal.lfilter(b, a, self.__measurementSignal)
     # step 3 vedran, downsample the signal
     maessignalsampled = signal.decimate(maessignalFilter, 10, ftype='iir')
     # step 4 with signal.freqs(b,a,y) we obtain frequency response of the signal
     freqHz, amplitude = signal.freqs(b, a, maessignalsampled)
     print 'angular frequency ', freqHz
     print 'amplitude response ', amplitude
     for vfreq, vdamp in zip(freqHz, amplitude):
         mode= EigenValue(vfreq,vdamp)
         self.__simulationModes.append(mode)
     # step 4 vedran, armax, _signal.real or signal.magnitude and signal.sampling data
     sys_ident= smapi.tsa.ARMA(simsignalsampled, order=(self.__order,self.__order)).fit()
     print sys_ident
     sys_ident= smapi.tsa.ARMA(maessignalsampled, order=(self.__order,self.__order)).fit()
     print sys_ident

def draw_frequency_response(b,a):
    from scipy.signal import freqs
    w, h = freqs(b, a)
    plt.figure(figsize=(8, 6), dpi=100)
    plt.plot(w, 20 * np.log10(abs(h)))
    plt.xscale('log')
    plt.title('Butterworth filter frequency response')
    plt.xlabel('Frequency [radians / second]')
    plt.ylabel('Amplitude [dB]')
    plt.margins(0, 0.1)
    plt.grid(which='both', axis='both')
    plt.axvline(100, color='green') # cutoff frequency
    plt.savefig('frequency_response.png',dpi=500)