Python signal.freqs函数代码示例

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 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 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 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()