For the impatient, here are the interesting bits:
def resample_poly_filter(up, down, beta=5.0, L=16001): # *** this block STOLEN FROM scipy.signal.resample_poly *** # Determine our up and down factors # Use a rational approximation to save computation time on really long # signals g_ = gcd(up, down) up //= g_ down //= g_ max_rate = max(up, down) sfact = np.sqrt(1+(beta/np.pi)**2) # generate first filter attempt: with 6dB attenuation at f_c. filt = firwin(L, 1/max_rate, window=('kaiser', beta)) N_FFT = 2**19 NBINS = N_FFT/2+1 paddedfilt = np.zeros(N_FFT) paddedfilt[:L] = filt ffilt = np.fft.rfft(paddedfilt) # now find the minimum between f_c and f_c+sqrt(1+(beta/pi)^2)/L bot = int(np.floor(NBINS/max_rate)) top = int(np.ceil(NBINS*(1/max_rate + 2*sfact/L))) firstnull = (np.argmin(np.abs(ffilt[bot:top])) + bot)/NBINS # generate the proper shifted filter filt2 = firwin(L, -firstnull+2/max_rate, window=('kaiser', beta)) return filt2 plt.figure(figsize=(15,3)) wfilt = resample_poly_filter(P, Q, L=2**16+1) plt.specgram(scipy_signal.resample_poly(sig, P, Q, window=wfilt)*30, scale='dB', Fs=P, NFFT=256) plt.colorbar() plt.axis((0,2,0,Q/2))
Recycling my test sweep from the previous post, I get:
Sweep resampled using my own filter |
No comments:
Post a Comment