sk_dsp_comm/iir_design_helper.py
def freqz_resp_list(b,a=np.array([1]),mode = 'dB',fs=1.0,Npts = 1024,fsize=(6,4)):
    """
    A method for displaying digital filter frequency response magnitude,
    phase, and group delay. A plot is produced using matplotlib

    freq_resp(self,mode = 'dB',Npts = 1024)

    A method for displaying the filter frequency response magnitude,
    phase, and group delay. A plot is produced using matplotlib

    freqz_resp(b,a=[1],mode = 'dB',Npts = 1024,fsize=(6,4))

        b = ndarray of numerator coefficients
        a = ndarray of denominator coefficents
     mode = display mode: 'dB' magnitude, 'phase' in radians, or 
            'groupdelay_s' in samples and 'groupdelay_t' in sec, 
            all versus frequency in Hz
     Npts = number of points to plot; default is 1024
    fsize = figure size; defult is (6,4) inches

    Mark Wickert, January 2015
    """
    if type(b) == list:
        # We have a list of filters
        N_filt = len(b)
    f = np.arange(0,Npts)/(2.0*Npts)
    for n in range(N_filt):
        w,H = signal.freqz(b[n],a[n],2*np.pi*f)
        if n == 0:
            plt.figure(figsize=fsize)
        if mode.lower() == 'db':
            plt.plot(f*fs,20*np.log10(np.abs(H)))
            if n == N_filt-1:
                plt.xlabel('Frequency (Hz)')
                plt.ylabel('Gain (dB)')
                plt.title('Frequency Response - Magnitude')

        elif mode.lower() == 'phase':
            plt.plot(f*fs,np.angle(H))
            if n == N_filt-1:
                plt.xlabel('Frequency (Hz)')
                plt.ylabel('Phase (rad)')
                plt.title('Frequency Response - Phase')

        elif (mode.lower() == 'groupdelay_s') or (mode.lower() == 'groupdelay_t'):
            """
            Notes
            -----

            Since this calculation involves finding the derivative of the
            phase response, care must be taken at phase wrapping points 
            and when the phase jumps by +/-pi, which occurs when the 
            amplitude response changes sign. Since the amplitude response
            is zero when the sign changes, the jumps do not alter the group 
            delay results.
            """
            theta = np.unwrap(np.angle(H))
            # Since theta for an FIR filter is likely to have many pi phase
            # jumps too, we unwrap a second time 2*theta and divide by 2
            theta2 = np.unwrap(2*theta)/2.
            theta_dif = np.diff(theta2)
            f_diff = np.diff(f)
            Tg = -np.diff(theta2)/np.diff(w)
            # For gain almost zero set groupdelay = 0
            idx = np.nonzero(np.ravel(20*np.log10(H[:-1]) < -400))[0]
            Tg[idx] = np.zeros(len(idx))
            max_Tg = np.max(Tg)
            #print(max_Tg)
            if mode.lower() == 'groupdelay_t':
                max_Tg /= fs
                plt.plot(f[:-1]*fs,Tg/fs)
                plt.ylim([0,1.2*max_Tg])
            else:
                plt.plot(f[:-1]*fs,Tg)
                plt.ylim([0,1.2*max_Tg])
            if n == N_filt-1:
                plt.xlabel('Frequency (Hz)')
                if mode.lower() == 'groupdelay_t':
                    plt.ylabel('Group Delay (s)')
                else:
                    plt.ylabel('Group Delay (samples)')
                plt.title('Frequency Response - Group Delay')
        else:
            s1 = 'Error, mode must be "dB", "phase, '
            s2 = '"groupdelay_s", or "groupdelay_t"'
            print(s1 + s2)
def freqz_resp_list(b,a=np.array([1]),mode = 'dB',fs=1.0,Npts = 1024,fsize=(6,4)):
    """
    A method for displaying digital filter frequency response magnitude,
    phase, and group delay. A plot is produced using matplotlib

    freq_resp(self,mode = 'dB',Npts = 1024)

    A method for displaying the filter frequency response magnitude,
    phase, and group delay. A plot is produced using matplotlib

    freqz_resp(b,a=[1],mode = 'dB',Npts = 1024,fsize=(6,4))

        b = ndarray of numerator coefficients
        a = ndarray of denominator coefficents
     mode = display mode: 'dB' magnitude, 'phase' in radians, or 
            'groupdelay_s' in samples and 'groupdelay_t' in sec, 
            all versus frequency in Hz
     Npts = number of points to plot; default is 1024
    fsize = figure size; defult is (6,4) inches

    Mark Wickert, January 2015
    """
    if type(b) == list:
        # We have a list of filters
        N_filt = len(b)
    f = np.arange(0,Npts)/(2.0*Npts)
    for n in range(N_filt):
        w,H = signal.freqz(b[n],a[n],2*np.pi*f)
        if n == 0:
            plt.figure(figsize=fsize)
        if mode.lower() == 'db':
            plt.plot(f*fs,20*np.log10(np.abs(H)))
            if n == N_filt-1:
                plt.xlabel('Frequency (Hz)')
                plt.ylabel('Gain (dB)')
                plt.title('Frequency Response - Magnitude')

        elif mode.lower() == 'phase':
            plt.plot(f*fs,np.angle(H))
            if n == N_filt-1:
                plt.xlabel('Frequency (Hz)')
                plt.ylabel('Phase (rad)')
                plt.title('Frequency Response - Phase')

        elif (mode.lower() == 'groupdelay_s') or (mode.lower() == 'groupdelay_t'):
            """
            Notes
            -----

            Since this calculation involves finding the derivative of the
            phase response, care must be taken at phase wrapping points 
            and when the phase jumps by +/-pi, which occurs when the 
            amplitude response changes sign. Since the amplitude response
            is zero when the sign changes, the jumps do not alter the group 
            delay results.
            """
            theta = np.unwrap(np.angle(H))
            # Since theta for an FIR filter is likely to have many pi phase
            # jumps too, we unwrap a second time 2*theta and divide by 2
            theta2 = np.unwrap(2*theta)/2.
            theta_dif = np.diff(theta2)
            f_diff = np.diff(f)
            Tg = -np.diff(theta2)/np.diff(w)
            # For gain almost zero set groupdelay = 0
            idx = np.nonzero(np.ravel(20*np.log10(H[:-1]) < -400))[0]
            Tg[idx] = np.zeros(len(idx))
            max_Tg = np.max(Tg)
            #print(max_Tg)
            if mode.lower() == 'groupdelay_t':
                max_Tg /= fs
                plt.plot(f[:-1]*fs,Tg/fs)
                plt.ylim([0,1.2*max_Tg])
            else:
                plt.plot(f[:-1]*fs,Tg)
                plt.ylim([0,1.2*max_Tg])
            if n == N_filt-1:
                plt.xlabel('Frequency (Hz)')
                if mode.lower() == 'groupdelay_t':
                    plt.ylabel('Group Delay (s)')
                else:
                    plt.ylabel('Group Delay (samples)')
                plt.title('Frequency Response - Group Delay')
        else:
            s1 = 'Error, mode must be "dB", "phase, '
            s2 = '"groupdelay_s", or "groupdelay_t"'
            print(s1 + s2)
def freqz_resp_list ( b , a = np . array ( [ 1 ] ) , mode = 'dB' , fs = 1.0 , Npts = 1024 , fsize = ( 6 , 4 ) ) : if type ( b ) == list : # We have a list of filters
 N_filt = len ( b ) f = np . arange ( 0 , Npts ) / ( 2.0 * Npts ) for n in range ( N_filt ) : w , H = signal . freqz ( b [ n ] , a [ n ] , 2 * np . pi * f ) if n == 0 : plt . figure ( figsize = fsize ) if mode . lower ( ) == 'db' : plt . plot ( f * fs , 20 * np . log10 ( np . abs ( H ) ) ) if n == N_filt - 1 : plt . xlabel ( 'Frequency (Hz)' ) plt . ylabel ( 'Gain (dB)' ) plt . title ( 'Frequency Response - Magnitude' ) elif mode . lower ( ) == 'phase' : plt . plot ( f * fs , np . angle ( H ) ) if n == N_filt - 1 : plt . xlabel ( 'Frequency (Hz)' ) plt . ylabel ( 'Phase (rad)' ) plt . title ( 'Frequency Response - Phase' ) elif ( mode . lower ( ) == 'groupdelay_s' ) or ( mode . lower ( ) == 'groupdelay_t' ) : """
            Notes
            -----

            Since this calculation involves finding the derivative of the
            phase response, care must be taken at phase wrapping points 
            and when the phase jumps by +/-pi, which occurs when the 
            amplitude response changes sign. Since the amplitude response
            is zero when the sign changes, the jumps do not alter the group 
            delay results.
            """ theta = np . unwrap ( np . angle ( H ) ) # Since theta for an FIR filter is likely to have many pi phase
 # jumps too, we unwrap a second time 2*theta and divide by 2
 theta2 = np . unwrap ( 2 * theta ) / 2. theta_dif = np . diff ( theta2 ) f_diff = np . diff ( f ) Tg = - np . diff ( theta2 ) / np . diff ( w ) # For gain almost zero set groupdelay = 0
 idx = np . nonzero ( np . ravel ( 20 * np . log10 ( H [ : - 1 ] ) < - 400 ) ) [ 0 ] Tg [ idx ] = np . zeros ( len ( idx ) ) max_Tg = np . max ( Tg ) #print(max_Tg)
 if mode . lower ( ) == 'groupdelay_t' : max_Tg /= fs plt . plot ( f [ : - 1 ] * fs , Tg / fs ) plt . ylim ( [ 0 , 1.2 * max_Tg ] ) else : plt . plot ( f [ : - 1 ] * fs , Tg ) plt . ylim ( [ 0 , 1.2 * max_Tg ] ) if n == N_filt - 1 : plt . xlabel ( 'Frequency (Hz)' ) if mode . lower ( ) == 'groupdelay_t' : plt . ylabel ( 'Group Delay (s)' ) else : plt . ylabel ( 'Group Delay (samples)' ) plt . title ( 'Frequency Response - Group Delay' ) else : s1 = 'Error, mode must be "dB", "phase, ' s2 = '"groupdelay_s", or "groupdelay_t"' print ( s1 + s2 )
A method for displaying digital filter frequency response magnitude,
    phase, and group delay. A plot is produced using matplotlib

    freq_resp(self,mode = 'dB',Npts = 1024)

    A method for displaying the filter frequency response magnitude,
    phase, and group delay. A plot is produced using matplotlib

    freqz_resp(b,a=[1],mode = 'dB',Npts = 1024,fsize=(6,4))

        b = ndarray of numerator coefficients
        a = ndarray of denominator coefficents
     mode = display mode: 'dB' magnitude, 'phase' in radians, or 
            'groupdelay_s' in samples and 'groupdelay_t' in sec, 
            all versus frequency in Hz
     Npts = number of points to plot; default is 1024
    fsize = figure size; defult is (6,4) inches

    Mark Wickert, January 2015
A method for displaying digital filter frequency response magnitude phase and group delay . A plot is produced using matplotlib freq_resp ( self mode = dB Npts = 1024 ) A method for displaying the filter frequency response magnitude phase and group delay . A plot is produced using matplotlib freqz_resp ( b a = [ 1 ] mode = dB Npts = 1024 fsize = ( 6 4 )) b = ndarray of numerator coefficients a = ndarray of denominator coefficents mode = display mode : dB magnitude phase in radians or groupdelay_s in samples and groupdelay_t in sec all versus frequency in Hz Npts = number of points to plot ; default is 1024 fsize = figure size ; defult is ( 6 4 ) inches Mark Wickert January 2015