[SciPy-User] problem with signal.group_delay

classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

[SciPy-User] problem with signal.group_delay

Neal Becker
Here is my test code (opt.loop_bw = 100e3, sps=2)

    from scipy import signal
    bw = opt.loop_bw/(sps*250e6)
    sos = signal.ellip (8, 0.5, 80, 0.5 * bw, output='sos')
    bas = [signal.sos2tf(s[np.newaxis,:]) for s in sos]
## my group delay code
    from group_delay import group_delay2  << my group_delay code
    delays = [group_delay2 (ba[0], ba[1])[0] for ba in bas]
## scipy group delay code
    delays2 = [signal.group_delay ((ba[0], ba[1]), w=1e-3) for ba in bas]
##
    print ('delays:', delays)
    print ('delay:', sum(delays))
    delay = sum(delays)

scipy group delay code gives me:
/home/nbecker/anaconda3/envs/py35/lib/python3.5/site-
packages/scipy/signal/filter_design.py:343: UserWarning: The group delay is
singular at frequencies [0.001], setting to 0
  format(", ".join("{0:.3f}".format(ws) for ws in w[singular]))

I have my own group delay code, that I wrote based on
https://ccrma.stanford.edu/~jos/fp/Numerical_Computation_Group_Delay.html

This seems to give a correct result.  My code is attached.  It is using my
own wrapper for fftw for fft, but of course you can use whatever you like.
,----[ /home/nbecker/sigproc/group_delay.py ]
| ## delay = real[dft(n*h(n))/dft(h(n))]
| '''https://ccrma-www.stanford.edu/~jos/filters/Numerical_Computation_Group_Delay.html'''
|
| import numpy as np
| from numpy_fncs import rotate
|
| def group_delay (a, fftsize=256):
|     from fft3 import fft, FORWARD
|     the_fft = fft (fftsize, FORWARD)
|     if (len(a) < the_fft.size):
|         padded = np.concatenate ((a, np.zeros (the_fft.size-len(a),
|         dtype=complex)))
|     else:
|         padded = a
|     ## num = rotate (the_fft (np.arange (the_fft.size) * padded),
|     ## the_fft.size/2) den = rotate (the_fft (padded), the_fft.size/2)
|     num = the_fft (np.arange (the_fft.size) * padded)
|     den = the_fft (padded)
|     poles = abs (den) < 1e-50
|     num[poles] = 0
|     den[poles] = 1
|     delay = (num / den).real
|     return delay
|
| def group_delay2 (a, b, fftsize=256):
|     c = np.convolve (a, np.conjugate (b)[::-1])
|     return group_delay (c) - len (b)
|
| if __name__ == "__main__":
|     ## a = np.array ((0,0,1,0,0),dtype=complex)
|     ## d = group_delay (a)
|     ## print d
|     ## b = np.array ((1,0),dtype=complex)
|     ## print group_delay2 (a, b)
|     from scipy.signal import butter
|     a,b = butter (7, .4)
|     ##a,b = ((0,1), (1,))
|     d = group_delay2 (a, b)
|     from pylab import plot, show
|     plot (np.linspace (0, 1, len(d)), d)
|     show()
|    
`----

_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user
Loading...