Python signal.freqz方法代码示例

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import freqz [as 别名]
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 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import freqz [as 别名]
def lsf2mfbe(lsf, mel_filters):

    NFFT = 512
    M = get_filterbank(n_filters=mel_filters, NFFT=NFFT, normalize=False, htk=True)

    mfbe = np.zeros(( len(lsf), mel_filters), dtype=np.float64)
    spec = np.zeros((len(lsf), NFFT/2+1,), dtype=np.float64)
   
    x = np.zeros((NFFT,), dtype=np.float64)
    x[0] = 1.0
    b = np.ones((1,), dtype=np.float64)
 
    for i, lsf_vec in enumerate(lsf):     
        #convert lsf to filter polynomial
        a_poly = lsf2poly(lsf_vec)
        # compute power spectrum
        w, H = freqz(b=1.0, a=a_poly, worN=NFFT, whole=True)
        spec_vec = np.abs(H[:(NFFT/2+1)])
        #spec_vec = np.square(spec_vec)
        # apply filterbank matrix
        mfbe[i,:] = np.log10( np.dot(M,spec_vec) )
        spec[i,:] = spec_vec
         
    return mfbe, spec 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import freqz [as 别名]
def test_ticket1441(self):
        """Regression test for ticket 1441."""
        # Because freqz previously used arange instead of linspace,
        # when N was large, it would return one more point than
        # requested.
        N = 100000
        w, h = freqz([1.0], worN=N)
        assert_equal(w.shape, (N,)) 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import freqz [as 别名]
def test_basic(self):
        w, h = freqz([1.0], worN=8)
        assert_array_almost_equal(w, np.pi * np.arange(8.0) / 8)
        assert_array_almost_equal(h, np.ones(8)) 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import freqz [as 别名]
def test_basic_whole(self):
        w, h = freqz([1.0], worN=8, whole=True)
        assert_array_almost_equal(w, 2 * np.pi * np.arange(8.0) / 8)
        assert_array_almost_equal(h, np.ones(8)) 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import freqz [as 别名]
def test_plot(self):

        def plot(w, h):
            assert_array_almost_equal(w, np.pi * np.arange(8.0) / 8)
            assert_array_almost_equal(h, np.ones(8))

        assert_raises(ZeroDivisionError,
                      freqz, [1.0], worN=8, plot=lambda w, h: 1 / 0)
        freqz([1.0], worN=8, plot=plot) 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import freqz [as 别名]
def test_absorption_filter(self):
        b = uwa.absorption_filter(200000)
        w, h = sp.freqz(b, 1, 4)
        h = utils.mag2db(np.abs(h))
        self.assertEqual(list(np.round(h)), [0.0, -3.0, -11.0, -22.0]) 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import freqz [as 别名]
def filter2spectrum(self, param):
        return np.abs(freqz([1], param)[1]) 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import freqz [as 别名]
def freqz_cas(sos,w):
    """
    Cascade frequency response
    
    Mark Wickert October 2016
    """
    Ns,Mcol = sos.shape
    w,Hcas = signal.freqz(sos[0,:3],sos[0,3:],w)
    for k in range(1,Ns):
        w,Htemp = signal.freqz(sos[k,:3],sos[k,3:],w)
        Hcas *= Htemp
    return w, Hcas 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import freqz [as 别名]
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] 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import freqz [as 别名]
def arma_periodogram(ar, ma, worN=None, whole=0):
    '''periodogram for ARMA process given by lag-polynomials ar and ma

    Parameters
    ----------
    ar : array_like
        autoregressive lag-polynomial with leading 1 and lhs sign
    ma : array_like
        moving average lag-polynomial with leading 1
    worN : {None, int}, optional
        option for scipy.signal.freqz   (read "w or N")
        If None, then compute at 512 frequencies around the unit circle.
        If a single integer, the compute at that many frequencies.
        Otherwise, compute the response at frequencies given in worN
    whole : {0,1}, optional
        options for scipy.signal.freqz
        Normally, frequencies are computed from 0 to pi (upper-half of
        unit-circle.  If whole is non-zero compute frequencies from 0 to 2*pi.

    Returns
    -------
    w : array
        frequencies
    sd : array
        periodogram, spectral density

    Notes
    -----
    Normalization ?

    This uses signal.freqz, which does not use fft. There is a fft version
    somewhere.

    '''
    w, h = signal.freqz(ma, ar, worN=worN, whole=whole)
    sd = np.abs(h)**2/np.sqrt(2*np.pi)
    if np.sum(np.isnan(h)) > 0:
        # this happens with unit root or seasonal unit root'
        print('Warning: nan in frequency response h, maybe a unit root')
    return w, sd 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import freqz [as 别名]
def mfbe2lsf(mfbe, lsf_order):

    NFFT = 512
    M = get_filterbank(n_filters=mfbe.shape[1], NFFT=NFFT, normalize=False, htk=True)

    M_inv = pinv(M)
    p = lsf_order

    lsf = np.zeros(( len(mfbe), lsf_order), dtype=np.float64)
    spec = np.zeros((len(mfbe), NFFT/2+1), dtype=np.float64)

    for i, mfbe_vec in enumerate(mfbe):
    
        # invert mel filterbank
        spec_vec = np.dot(M_inv, np.power(10, mfbe_vec))

        # floor reconstructed spectrum
        spec_vec = np.maximum(spec_vec, 1e-9)
 
        # squared magnitude 2-sided spectrum
        twoside = np.r_[spec_vec, np.flipud(spec_vec[1:-1])]
        twoside = np.square(twoside) 
        r = np.fft.ifft(twoside)
        r = r.real

        # reference from talkbox
        # a,_,_ = TB.levinson(r, order=p)
  
        # levinson-durbin
        a = LA.solve_toeplitz(r[0:p],r[1:p+1])
        a = np.r_[1.0, -1.0*a]
   
        lsf[i,:] = poly2lsf(a)
   
        # reconstructed all-pole spectrum
        w, H = freqz(b=1.0, a=a, worN=NFFT, whole=True)
        spec[i,:] = np.abs(H[:(NFFT/2+1)])
            
    return lsf, spec 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import freqz [as 别名]
def spec2lsf(spec, lsf_order=30):

    NFFT = 2*(spec.shape[0]-1)
    n_frames = spec.shape[1]

    p = lsf_order

    lsf = np.zeros(( n_frames, lsf_order), dtype=np.float64)
    spec_rec = np.zeros(spec.shape)

    for i, spec_vec in enumerate(spec.T):
    
        # floor reconstructed spectrum
        spec_vec = np.maximum(spec_vec, 1e-9)
 
        # squared magnitude 2-sided spectrum
        twoside = np.r_[spec_vec, np.flipud(spec_vec[1:-1])]
        twoside = np.square(twoside) 
        r = np.fft.ifft(twoside)
        r = r.real
  
        # levinson-durbin
        a = LA.solve_toeplitz(r[0:p],r[1:p+1])
        a = np.r_[1.0, -1.0*a]
   
        lsf[i,:] = poly2lsf(a)
   
        # reconstructed all-pole spectrum
        w, H = freqz(b=1.0, a=a, worN=NFFT, whole=True)
        spec_rec[:,i] = np.abs(H[:(NFFT/2+1)])
            
    return lsf, spec_rec 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import freqz [as 别名]
def arma_periodogram(ar, ma, worN=None, whole=0):
    """
    Periodogram for ARMA process given by lag-polynomials ar and ma

    Parameters
    ----------
    ar : array_like
        autoregressive lag-polynomial with leading 1 and lhs sign
    ma : array_like
        moving average lag-polynomial with leading 1
    worN : {None, int}, optional
        option for scipy.signal.freqz (read "w or N")
        If None, then compute at 512 frequencies around the unit circle.
        If a single integer, the compute at that many frequencies.
        Otherwise, compute the response at frequencies given in worN
    whole : {0,1}, optional
        options for scipy.signal.freqz
        Normally, frequencies are computed from 0 to pi (upper-half of
        unit-circle.  If whole is non-zero compute frequencies from 0 to 2*pi.

    Returns
    -------
    w : array
        frequencies
    sd : array
        periodogram, spectral density

    Notes
    -----
    Normalization ?

    This uses signal.freqz, which does not use fft. There is a fft version
    somewhere.
    """
    w, h = signal.freqz(ma, ar, worN=worN, whole=whole)
    sd = np.abs(h) ** 2 / np.sqrt(2 * np.pi)
    if np.any(np.isnan(h)):
        # this happens with unit root or seasonal unit root'
        import warnings
        warnings.warn('Warning: nan in frequency response h, maybe a unit '
                      'root', RuntimeWarning)
    return w, sd 

# 需要导入模块: from scipy import signal [as 别名]
# 或者: from scipy.signal import freqz [as 别名]
def compute_frequency_response(filter_coefs, a_vals, fs):
    """Compute the frequency response of a filter.

    Parameters
    ----------
    filter_coefs : 1d or 2d array
        If 1d, interpreted as the B-value filter coefficients.
        If 2d, interpreted as the second-order (sos) filter coefficients.
    a_vals : 1d array or None
        The A-value filter coefficients for a filter.
        If second-order filter coefficients are provided in `filter_coefs`, must be None.
    fs : float
        Sampling rate, in Hz.

    Returns
    -------
    f_db : 1d array
        Frequency vector corresponding to attenuation decibels, in Hz.
    db : 1d array
        Degree of attenuation for each frequency specified in `f_db`, in dB.

    Examples
    --------
    Compute the frequency response for an FIR filter:

    >>> from neurodsp.filt.fir import design_fir_filter
    >>> filter_coefs = design_fir_filter(fs=500, pass_type='bandpass', f_range=(8, 12))
    >>> f_db, db = compute_frequency_response(filter_coefs, 1, fs=500)

    Compute the frequency response for an IIR filter, which uses SOS coefficients:

    >>> from neurodsp.filt.iir import design_iir_filter
    >>> sos_coefs = design_iir_filter(fs=500, pass_type='bandpass',
    ...                               f_range=(8, 12), butterworth_order=3)
    >>> f_db, db = compute_frequency_response(sos_coefs, None, fs=500)
    """

    if filter_coefs.ndim == 1 and a_vals is not None:
        # Compute response for B & A value filter coefficient inputs
        w_vals, h_vals = freqz(filter_coefs, a_vals, worN=int(fs * 2))
    elif filter_coefs.ndim == 2 and a_vals is None:
        # Compute response for sos filter coefficient inputs
        w_vals, h_vals = sosfreqz(filter_coefs, worN=int(fs * 2))
    else:
        raise ValueError("The organization of the filter coefficient inputs is not understood.")

    f_db = w_vals * fs / (2. * np.pi)
    db = 20 * np.log10(abs(h_vals))

    return f_db, db