Python scipy.signal 模块，freqz() 实例源码

def freqz_(sys, w, dt=8e-9):
    """
    This function computes the frequency response of a zpk system at an
    array of frequencies.

    It loosely mimicks 'scipy.signal.frequresp'.

    Parameters
    ----------
    system: (zeros, poles, k)
        zeros and poles both in rad/s, k is the actual coefficient, not DC gain
    w: np.array
        frequencies in rad/s
    dt: sampling time

    Returns
    -------
    np.array(..., dtype=np.complex) with the response
    """
    z, p, k = sys
    b, a = sig.zpk2tf(z, p, k)
    _, h = sig.freqz(b, a, worN=w*dt)
    return h

def frequencyResponse(self, freqs=None):
        if self.bandType == 'allpass':
            return spsig.freqz(1, worN=freqs)
        if self.bandType == 'allstop':
            return spsig.freqz(0, worN=freqs)

        numCoef = self.numCoef
        denomCoef = self.denomCoef

        if self.zeroPhase:
            # http://www.mathworks.com/matlabcentral/newsreader/view_thread/245017
            numCoef = np.convolve(numCoef,numCoef[::-1])
            denomCoef = np.convolve(denomCoef,denomCoef[::-1])

        # freqz does not preserve dtype of arguments, report bug XXX - idfah
        return spsig.freqz(numCoef, denomCoef, worN=freqs)

def another():
# plot the figure of signals
  fig = plt.figure(1)
  ax = fig.add_subplot(311)
  ax.set_title('input signal')
  ax.plot(t[1:Npts],x[1:Npts])    # just plot part of the signal
  ax = fig.add_subplot(312)
  ax.set_title('expected signal')
  ax.plot(t[1:Npts],xo[1:Npts])
  ax = fig.add_subplot(313)
  ax.set_title('filtered signal')
  ax.plot(t[1:Npts],y[1:Npts])
  fig.savefig('pic1.png')

# plot the mag & phase response of the LPF
  w,h = signal.freqz(b,1)
  hdb = 20 * np.log10(np.abs(h))
  hphs = np.unwrap(np.arctan2(np.imag(h),np.real(h)))
  fig2=plt.figure(2)
  ax2 = fig2.add_subplot(211)
  ax2.set_title('frequency response')
  ax2.plot(w/max(w),hdb)
  ax2 = fig2.add_subplot(212)
  ax2.set_title('phase response')
  ax2.plot(w/max(w),hphs)
  fig2.savefig('pic2.png')

def test_absorption_filter(self):
        b = uwa.absorption_filter(200000)
        w, h = sp.freqz(b, 1, 4)
        h = 20*np.log10(np.abs(h))
        self.assertEqual(list(np.round(h)), [0.0, -2.0, -8.0, -17.0])

def freqz(b, a=1, fs=2.0, worN=None, whole=False):
    """Plot frequency response of a filter.

    This is a convenience function to plot frequency response, and internally uses
    :func:`scipy.signal.freqz` to estimate the response. For further details, see the
    documentation for :func:`scipy.signal.freqz`.

    :param b: numerator of a linear filter
    :param a: denominator of a linear filter
    :param fs: sampling rate in Hz (optional, normalized frequency if not specified)
    :param worN: see :func:`scipy.signal.freqz`
    :param whole: see :func:`scipy.signal.freqz`
    :returns: (frequency vector, frequency response vector)

    >>> import arlpy
    >>> arlpy.signal.freqz([1,1,1,1,1], fs=120000);
    >>> w, h = arlpy.signal.freqz([1,1,1,1,1], fs=120000)
    """
    import matplotlib.pyplot as plt
    w, h = _sig.freqz(b, a, worN, whole)
    f = w*fs/(2*_np.pi)
    fig = plt.figure()
    ax1 = fig.add_subplot(111)
    plt.plot(f, 20*_np.log10(abs(h)), 'b')
    plt.ylabel('Amplitude [dB]', color='b')
    plt.xlabel('Frequency [Hz]')
    plt.grid()
    ax1.twinx()
    angles = _np.unwrap(_np.angle(h))
    plt.plot(f, angles, 'g')
    plt.ylabel('Angle (radians)', color='g')
    plt.axis('tight')
    plt.show()
    return w, h

def plotfilt(b, fs:int, ofn=None):
    if fs is None:
        fs=1 #normalized freq

    L = b.size

    fg,axs = subplots(2,1,sharex=False)
    freq, response = signal.freqz(b)
    response_dB = 20*np.log10(abs(response))
    if response_dB.max()>0:
        logging.error('filter may be unstable')

    axs[0].plot(freq*fs/(2*np.pi), response_dB)
    axs[0].set_title(f'filter response  {L} taps')
    axs[0].set_ylim((-100,None))
    axs[0].set_ylabel('|H| [db]')
    axs[0].set_xlabel('frequency [Hz]')

    t = np.arange(0, L/fs, 1/fs)
    axs[1].plot(t,b)
    axs[1].set_xlabel ('time [sec]')
    axs[1].set_title('impulse response')
    axs[1].set_ylabel('amplitude')
    axs[1].autoscale(True,tight=True)

    fg.tight_layout()

    if ofn:
        ofn = Path(ofn).expanduser()
        ofn = ofn.with_suffix('.png')
        print('writing',ofn)
        fg.savefig(str(ofn),dpi=100,bbox_inches='tight')

def frequencyResponse(self, freqs=None):
        return spsig.freqz(self.impulseResponse, worN=freqs)

def sosfreqz(sos, ws = None):
    if ws is None:
        H = [1] * 512        
    else:
        H = [1] * len(ws)

    for i in range(len(sos)):
        w, h = freqz(sos[i,:3], sos[i, 3:], worN = ws)
        H *= h
    return w, H

def MOEar(self, correctionType = 'ELC'):
        """ Method to approximate middle-outer ear transfer function for linearly scaled
            frequency representations, using an FIR approximation of order 600 taps.
            As appears in :
            - A. Härmä, and K. Palomäki, ''HUTear – a free Matlab toolbox for modeling of human hearing'',
            in Proceedings of the Matlab DSP Conference, pp 96-99, Espoo, Finland 1999.
        Arguments          :
            correctionType : (string)     String which specifies the type of correction :
                                          'ELC' - Equal Loudness Curves at 60 dB (default)
                                          'MAP' - Minimum Audible Pressure at ear canal
                                          'MAF' - Minimum Audible Field
        Returns            :
            LTq            : (ndarray)    1D Array containing the transfer function, without the DC sub-band.
        """
        # Parameters
        firOrd = self.nfft
        Cr, fr, Crdb = self.OutMidCorrection(correctionType, firOrd, self.fs)
        Cr[self.nfft - 1] = 0.

        # FIR Design
        A = firwin2(firOrd, fr, Cr, nyq = self.fs/2)
        B = 1
        _, LTq = freqz(A, B, firOrd, self.fs)

        LTq = 20. * np.log10(np.abs(LTq))
        LTq -= max(LTq)
        return LTq[:self.nfft/2 + 1]

def _filter_resp(freqmin, freqmax, corners=2, zerophase=False, sr=None,
                 N=None, whole=False):
    """
    Complex frequency response of Butterworth-Bandpass Filter.

    :param freqmin: Pass band low corner frequency.
    :param freqmax: Pass band high corner frequency.
    :param corners: Filter corners
    :param zerophase: If True, apply filter once forwards and once backwards.
        This results in twice the number of corners but zero phase shift in
        the resulting filtered trace.
    :param sr: Sampling rate in Hz.
    :param N,whole: passed to scipy.signal.freqz

    :return: frequencies and complex response
    """
    df = sr
    fe = 0.5 * df
    low = freqmin / fe
    high = freqmax / fe
    # raise for some bad scenarios
    if high > 1:
        high = 1.0
        msg = "Selected high corner frequency is above Nyquist. " + \
              "Setting Nyquist as high corner."
        log.warning(msg)
    if low > 1:
        msg = "Selected low corner frequency is above Nyquist."
        raise ValueError(msg)
    [b, a] = iirfilter(corners, [low, high], btype='band',
                       ftype='butter', output='ba')
    freqs, values = freqz(b, a, N, whole=whole)
    if zerophase:
        values *= np.conjugate(values)
    return freqs, values

def _filter_attenuation(h, freq, gain):
    """Compute minimum attenuation at stop frequency"""
    from scipy.signal import freqz
    _, filt_resp = freqz(h.ravel(), worN=np.pi * freq)
    filt_resp = np.abs(filt_resp)  # use amplitude response
    filt_resp /= np.max(filt_resp)
    filt_resp[np.where(gain == 1)] = 0
    idx = np.argmax(filt_resp)
    att_db = -20 * np.log10(filt_resp[idx])
    att_freq = freq[idx]
    return att_db, att_freq

def tf_coefficients(self, frequencies=None, coefficients=None,
                                delay=False):
        """
        computes implemented transfer function - assuming no delay and
        infinite precision (actually floating-point precision)
        Returns the discrete transfer function realized by coefficients at
        frequencies.

        Parameters
        ----------
        coefficients: np.array
            coefficients as returned from iir module

        frequencies: np.array
            frequencies to compute the transfer function for

        dt: float
            discrete sampling time (seconds)

        zoh: bool
            If true, zero-order hold implementation is assumed. Otherwise,
            the delay is expected to depend on the index of biquad.

        Returns
        -------
        np.array(..., dtype=np.complex)
        """
        if frequencies is None:
            frequencies = self.frequencies
        frequencies = np.asarray(frequencies, dtype=np.float64)
        if coefficients is None:
            fcoefficients = self.coefficients
        else:
            fcoefficients = coefficients
        # discrete frequency
        w = frequencies * 2 * np.pi * self.dt * self.loops
        # the higher stages have progressively more delay to the output
        if delay:
            delay_per_cycle = np.exp(-1j * self.dt * frequencies * 2 * np.pi)
        h = np.zeros(len(w), dtype=np.complex128)
        for i in range(len(fcoefficients)):
            sos = np.asarray(fcoefficients[i], dtype=np.float64)
            # later we can use sig.sosfreqz (very recent function, dont want
            #  to update scipy now)
            ww, hh = sig.freqz(sos[:3], sos[3:], worN=np.asarray(w,
                                                           dtype=np.float64))
            if delay:
                hh *= delay_per_cycle ** i
            h += hh
        return h

def test_temporal_learning(flen=[5]):
    '''
    Function that computes CSP filters then uses those with the temporal filter MTL 
    idea, and confirms that the output has a spectral profile that is similar to expected.
    Generate y values from trace of filtered covariance to ensure that is not an issue
    '''
    def covmat(x,y):
        return (1/(x.shape[1]-1)*x.dot(y.T))
    D = scio.loadmat('/is/ei/vjayaram/code/python/pyMTL_MunichMI.mat')
    data = D['T'].ravel()
    labels = D['Y'].ravel()
    fig = plt.figure()
    fig.suptitle('Recovered frequency filters for various filter lengths')
    model = TemporalBRC(max_prior_iter=100)
    # compute cross-covariance matrices and generate X
    for find,freq in enumerate(flen):
        X = []
        for tind,d in enumerate(data):
            d = d.transpose((2,0,1))
            x = np.zeros((d.shape[0], freq))
            nsamples = d.shape[2]
            for ind, t in enumerate(d):
                for j in range(freq):
                    C = covmat(t[:,0:(nsamples-j)],t[:,j:])
                    x[ind,j] = np.trace(C + C.T)
            X.append(x)
            labels[tind] = labels[tind].ravel()

        model.fit_multi_task(X,labels)
        b = numerics.solve_fir_coef(model.prior.mu.flatten())[0]
        # look at filter properties
        w,h = freqz(b[1:],worN=100)
        w = w*500/2/np.pi # convert to Hz

        ax1 = fig.add_subplot(3,3,find+1)
        plt.plot(w, 20 * np.log10(abs(h)), 'b')
        plt.ylabel('Amplitude [dB]', color='b')
        plt.xlabel('Frequency [Hz]')
        ax2 = ax1.twinx()
        angles = np.unwrap(np.angle(h))
        plt.plot(w, angles, 'g')
        plt.ylabel('Angle (radians)', color='g')
        plt.grid()
        plt.axis('tight')
    plt.show()

def filter_freqs(lowcut, highcut, fs, plot=False, corners=4):
   '''
   The frequency band information should technically be the amplitude spectrum
   of the filtered seismograms, but the amplitude spectrum of the bandpass
   filter is sufficient. Here we use a bandpass butterworth filter.

   The function freq band takes 3 arguments:
           lowcut = low end of bandpass (in Hz)
           highcut = high end of bandpass (in Hz)
           fs      = sample rate (1.0/dt)

   It returns two vectors:
           omega = frequency axis (in rad/s)
           amp   = frequency response of the filter
   '''
   #Define bandpass parameters
   nyquist = 0.5 * fs
   fmin    = lowcut/nyquist
   fmax    = highcut/nyquist

   #Make filter shape
   b, a = iirfilter(corners, [fmin,fmax], btype='band', ftype='butter')

   #Determine freqency response
   freq_range = np.linspace(0,0.15,200)
   w, h = freqz(b,a,worN=freq_range)

   omega    = fs * w                # in rad/s
   omega_hz = (fs * w) / (2*np.pi)  # in hz
   amp   = abs(h)

   if(plot == True):
      #Checks----------------
      #plt.semilogx(omega,amp)
      #plt.axvline(1/10.0)
      #plt.axvline(1/25.0)
      #plt.axhline(np.sqrt(0.5))

      plt.plot(omega,amp)
      plt.xlabel('frequency (rad/s)')
      plt.ylabel('amplitude')
      plt.show()

   return  omega, amp

def get_filter_freqs(filter_type,freqmin,freqmax,sampling_rate,**kwargs):
   '''
   Returns frequency band information of the filter used in cross
   correlation delay time measurements.

   args-------------------------------------------------------------------------- 
   filter_type: type of filter used (e.g., bandpass, gaussian etc...)
   freqmin: minimum frequency of filter (in Hz)
   freqmax: maximum frequency of filter (in Hz)
   sampling_rate: sampling rate of seismograms (in Hz)
   kwargs-------------------------------------------------------------------------
   corners: number of corners used in filter (if bandpass). default = 2
   std_dev: standard deviation of gaussian filter (in Hz)
   plot: True or False

   returns-----------------------------------------------------------------------
   omega: frequency axis (rad/s)
   amp: frequency response of filter
   '''
   plot = kwargs.get('plot',False)
   corners = kwargs.get('corners',4)
   std_dev = kwargs.get('std_dev',0.1)
   mid_freq = kwargs.get('mid_freq',1/10.0)

   if filter_type == 'bandpass':
      nyquist = 0.5 * sampling_rate
      fmin = freqmin/nyquist
      fmax = freqmax/nyquist
      b, a = iirfilter(corners, [fmin,fmax], btype='band', ftype='butter')
      freq_range = np.linspace(0,0.15,200)
      w, h = freqz(b,a,worN=freq_range)
      omega = sampling_rate * w
      omega_hz = (sampling_rate * w) / (2*np.pi)
      amp = abs(h)

   elif filter_type == 'gaussian':
      fmin=freqmin
      fmax=freqmax
      omega_hz = np.linspace(0,0.5,200)
      omega = omega_hz*(2*np.pi)
      f_middle_hz = (fmin+fmax)/2.0
      f_middle = f_middle_hz*(2*np.pi)
      #f_middle_hz = mid_freq  #f_middle_hz = 10.0 was used in JGR manuscript
      #f_middle = f_middle_hz*(2*np.pi)
      print f_middle_hz,f_middle
      amp = np.exp(-1.0*((omega-f_middle)**2)/(2*(std_dev**2)))
      amp = amp/np.max(amp)

   if plot:
      fig,axes = plt.subplots(2)
      axes[0].plot(omega,amp)
      axes[0].set_xlabel('frequency (rad/s)')
      axes[0].set_ylabel('amplitude')
      axes[1].plot(omega_hz,amp)
      axes[1].set_xlabel('frequency (Hz)')
      axes[1].set_ylabel('amplitude')
      axes[1].axvline(fmin)
      axes[1].axvline(fmax)
      plt.show()

   return omega, amp

def get_filter_freqs(filter_type,freqmin,freqmax,sampling_rate,**kwargs):
   '''
   Returns frequency band information of the filter used in cross
   correlation delay time measurements.

   args-------------------------------------------------------------------------- 
   filter_type: type of filter used (e.g., bandpass, gaussian etc...)
   freqmin: minimum frequency of filter (in Hz)
   freqmax: maximum frequency of filter (in Hz)
   sampling_rate: sampling rate of seismograms (in Hz)
   kwargs-------------------------------------------------------------------------
   corners: number of corners used in filter (if bandpass). default = 2
   std_dev: standard deviation of gaussian filter (in Hz)
   plot: True or False

   returns-----------------------------------------------------------------------
   omega: frequency axis (rad/s)
   amp: frequency response of filter
   '''
   plot = kwargs.get('plot',False)
   corners = kwargs.get('corners',4)
   std_dev = kwargs.get('std_dev',0.1)
   mid_freq = kwargs.get('mid_freq',1/10.0)

   if filter_type == 'bandpass':
      nyquist = 0.5 * sampling_rate
      fmin = freqmin/nyquist
      fmax = freqmax/nyquist
      b, a = iirfilter(corners, [fmin,fmax], btype='band', ftype='butter')
      freq_range = np.linspace(0,0.15,200)
      w, h = freqz(b,a,worN=freq_range)
      omega = sampling_rate * w
      omega_hz = (sampling_rate * w) / (2*np.pi)
      amp = abs(h)

   elif filter_type == 'gaussian':
      fmin=freqmin
      fmax=freqmax
      omega_hz = np.linspace(0,0.5,200)
      omega = omega_hz*(2*np.pi)
      f_middle_hz = (fmin+fmax)/2.0
      f_middle = f_middle_hz*(2*np.pi)
      #f_middle_hz = mid_freq  #f_middle_hz = 10.0 was used in JGR manuscript
      #f_middle = f_middle_hz*(2*np.pi)
      print f_middle_hz,f_middle
      amp = np.exp(-1.0*((omega-f_middle)**2)/(2*(std_dev**2)))
      amp = amp/np.max(amp)

   if plot:
      fig,axes = plt.subplots(2)
      axes[0].plot(omega,amp)
      axes[0].set_xlabel('frequency (rad/s)')
      axes[0].set_ylabel('amplitude')
      axes[1].plot(omega_hz,amp)
      axes[1].set_xlabel('frequency (Hz)')
      axes[1].set_ylabel('amplitude')
      axes[1].axvline(fmin)
      axes[1].axvline(fmax)
      plt.show()

   return omega, amp